Criteria for assessment of the quality of Pseudomonas aeruginosa genome sequences

Full Text

Introduction

With the development of high-throughput sequencing technologies and the decrease in their cost, the volume of genomic data produced is growing exponentially. Projects using large datasets of whole-genome sequencing (WGS) have many advantages: statistical power is increased, and it becomes possible to test various hypotheses about the micro- and macroevolution of genomes.

The continuous improvement of sequencing technologies and bioinformatics analysis has increased the significance of WGS in biology, medicine, pharmaceuticals, and agriculture, stimulating comparative genomic research. However, the growth in the number of sequencing projects and laboratories has led to an increase in the number of genome assemblies that are not always suitable for analysis. This highlighted the need to evaluate the quality of whole genome assembly data for researchers, who use it. This in turn created a necessity to develop standard metrics for comparing the quality of genome assemblies and annotations, as well as for evaluating the effectiveness of different methods used to obtain them.

Recent studies on genome assembly quality assessment have focused either on pre-assembly quality control or on the assembly evaluation in terms of contiguity and correctness. However, the assessment of correctness depends on the reference and is not applicable to de novo assembly projects. Therefore, it is worth studying methods that allow for quality assessment reports to be obtained both after and before assembly, to check the quality/correctness of de novo assembly and input data [1].

For genome assemblies, metrics such as the number of contigs, the number of scaffolds, and N50 (the maximum contig’s length at which the total length of all contigs no shorter than this value accounts for at least 50% of the total length of all contigs in the assembly) provide only a brief overview of genome quality, not always reflecting its analytical suitability.

In turn, there are currently a sufficient number of resources and methods for the post-analysis stage of work, as well as for assessing genome quality: Picard [2], SQUAT [1], Plantagora [3], QUAST [4], CheckM1 [5], CheckM2 [6], GenomeQC [7], BUSCO [8]. However, they are unified and represent algorithms with different orientations, sometimes suitable for analyzing only eukaryotic organisms, while not taking into account the specific genome organization of a particular species. One of the most versatile and widely used instruments that utilizes genes to assess WGS data is BUSCO. Unlike the solutions mentioned above, BUSCO focuses on genome analysis using evolutionarily conserved orthologous genes, which are considered universal for certain taxonomic groups (bacteria, fungi, plants or animals). However, BUSCO does not provide an answer about the quality of the analyzed genome, only indicating the percentage of found/not found orthologous genes, and the final conclusion must be drawn by the specialist themselves. However, orthologous genes can be lost without affecting bacterial viability, unlike housekeeping genes, which can lead to an underestimation of genome quality.

Currently, WGS of infectious disease pathogens is widely used to study them, determine their origin, and track their spread. To assess the quality of such a large amount of data, domestic software tools are necessary. The latter is particularly important given that import substitution is becoming one of the strategic objectives in conditions where access to foreign technologies and foreign databases is difficult [9].

There are currently no evaluation criteria for WGS data for Pseudomonas aeruginosa. There are software services that perform assessments based on general (non-specific) criteria (N50, number of contigs, etc.) and do not take into account the characteristics of a specific microorganism.

The aim of the study is to create an algorithm and criteria for assessing the quality of WGS data from P. aeruginosa representatives, as well as to develop a domestic software capable of evaluating the quality of WGS data.

Materials and methods

The study used 108 genomes of P. aeruginosa strains: 24 strains were obtained from the Collection of Pathogenic Microorganisms of the Rostov-on-Don Antiplague Institute of Rospotrebnadzor (isolated in Rostov-on-Don, Khabarovsk, and Mariupol in 2022–2024), and 84 strains were obtained from the international NCBI database. WGS was conducted as part of the implementation of the federal project for the socio-economic development of the Russian Federation until 2030, “Sanitary Shield of the Country — Health Security (Prevention, Detection, Response)”. Sequencing was performed on the MiSeq platform (Illumina) using the MiSeq Reagent Kit v2 (500-cycles) (Illumina). This method allows for reads 2 × 251 nucleotides long, with genome coverage ranging from 8 to 20.

The assessment of the primary sequencing data was performed using the FastQC program. The collected WGS data was analyzed using the QUAST program [4, 10]. The Trimmomatic [11] and Lighter [12] algorithms were used for trimming and for reads correction. Genome assembly from reads was performed using the Spades program [13]. All genomes have passed an initial assessment using the Kraken 2 program, which allows for the identification of DNA fragment belonging to various prokaryotic species [14]. The WGS data of strain PAO1 from the international NCBI database [15] were used as a reference genome.

The proprietary software was developed in the Java and Python programming languages. The algorithm for searching for gene sequences in the assembly was performed with the use of Smith-Waterman local alignment with a minimum similarity threshold of 80%.

The confidence interval was calculated, differences were considered significant at p < 0.05.

Results

It is known that the genome of the causative agent of Pseudomonas aeruginosa infection contains a number of genes that are critical for its viability. These genes are called housekeeping genes. It is evident that if any of these genes are missing from the WGS data, it is a sequencing and/or genome assembly error. This very feature was the basis of our proposed algorithm – all essential genes should be detected in a good-quality sequence. Of course, the selection of genes to be used for quality control is of great importance in this process.

One of the criteria we devised for the algorithm to assess genome quality is the selection of genes based on the following criteria:

• nucleotide sequences must be within 1000 bp;
• the gene must be a single-copy;
• the gene must be directly involved in the microorganism's physiological activity, performing essential functions for its life processes;
• the gene must present in all strains of aeruginosa.

The oprI gene was chosen for rapid species identification of P. aeruginosa. The main task is to assess the quality of sequencing data not only based on the identification of essential genes but also on translating their sequences. Taking in account that these genes are critical for the existence of a microbial cell, their absence from the genome or critical translation errors (stop codons) are considered as sequencing errors.

Housekeeping genes used for validating the selected whole-genome sequences of P. aeruginosa: fur, algU, dinB, dnaQ, holA, holB, PA0472, fpvI, tonB1, cntL, sigX, capB, cspD, groES, rpoH. The selected genes are essential for functioning and survival in the environment and in a macroorganism. The following parameters were chosen as criteria for evaluating the whole-genome sequences: the GC content of the P. aeruginosa genomes, the size of the P. aeruginosa whole-genome sequence, and the N50 scaffold value.

After conducting the research and selecting the quality assessment criteria for genomes, the Genomes Validator software was developed, which, for convenience, operates in "batch mode," analyzes an unlimited number of genomes, and presents the results in tabular form. For each genome, the original file name, species, quality (poor, average, good), length, N50 value, GC content, as well as the reason for the invalidity of the genome are indicated (Fig. 1).

 

Fig. 1. Practical demonstration of the program Genomes Validator; the results of genome analysis are presented in table format.

 

The developed program Genomes Validator is a cross-platform, which has a graphical interface, does not require installation, allows to analyze multiple genomes at the same time, and is available for downloading at https://github.com/alexeyvod/GenomesValidator. It has an intuitive interface and is user-friendly for those without programming skills.

The program was validated on a sample of 108 whole-genome sequences of P. aeruginosa strains. Following validation, genomes of good (63%), medium (29%), and poor (8%) quality were identified. Further analysis identified 37% of the genomes analyzed (of medium and poor quality), which can help avoid errors in subsequent calculations using bioinformatics methods. The average N50 value among the sample was 1,250,527.

The parameters N50, genome length, and GC content were identical to the values obtained from the programs used for comparison: CheckM2 and QUAST. However, these programs do not provide genome quality assessment metrics.

When assessing the quality of bacterial genomes using CheckM2 software, we found out that genomes with a Completeness value of 100 showed significant variability in the Contamination index. At the same time, the use of the Genomes Validator made it possible to estimate additionally the size of the whole-genome sequence, which may be more informative for practical analysis of WGS data. (Fig. 2). This parameter allows for a preliminary assessment of the presence of extrachromosomal elements in the genome of the strain under study. It should be noted that contamination with foreign DNA usually affects the overall GC content and causes a significant change in genome size, whereas the presence of plasmids or other mobile genetic elements does not lead to significant changes in this parameter.

 

Fig. 2. Fragment of a table comparing the performance characteristics of the CheckM2 and Genomes Validator programs.

 

Based on the statistical analysis of the Completeness and Contamination parameters (Table), it was found that Contamination values in the range of 2 to 8 may indicate a possible low reliability of the obtained WGS data. However, such results could also be due to specific characteristics of the clinical isolate's genome. Thus, the genome of strain Ps-agn-2889, analyzed in the CheckM2 program, has a Completeness score of 100 with a Contamination score of –35.65, but the reason for the contamination is not clear from the data obtained. Analysis in the program Genomes Validator revealed a 1.5-fold increase in genome size and GC content, indicating clear contamination with foreign bacterial DNA. The genome of clinical strain 44269, analyzed using the CheckM2 program, has a Completeness score of 100 with a Contamination score of –12.04, which casts doubt on its quality. Nevertheless, when using the Genomes Validator program, the genome size and GC content indicate the clear presence of extrachromosomal elements that affect the Contamination score, rather than contamination with foreign DNA, as evidenced by the research of the strain authors [16].

 

Comparison of the quality metrics of the CheckM2 and Genomes Validator programs

«Genomes validator»	«CheckM2»
	completeness	contamination
Good/High	99.99 ± 0.001	0.98 ± 0.203
Average	83.89 ± 1.865	2.23 ± 0.222
Poor/Low	81.01 ± 6.667	8.72 ± 3.708

Note. The confidence interval is indicated at p < 0.05

 

Discussion

The oprI gene was chosen as the species-defining gene for several reasons: its nucleotide sequence is 253 bp long, which allows for species identification even with very poor WGS data quality; the OprI protein plays an important role in binding to peptidoglycan, participates in immunological reactions, and is responsible for susceptibility to antimicrobial peptides [17–19]. This gene was chosen because one meta-analysis showed that it is successfully used to identify the P. aeruginosa species with high accuracy [20].

Housekeeping genes (fur, algU, dinB, dnaQ, holA, holB, PA0472, fpvI, tonB1, cntL, sigX, capB, cspD, groES, rpoH) were selected for validation of the chosen whole-genome sequences of P. aeruginosa based on their functional significance, as determined by a literature data analysis.

The fur gene is the main regulator of iron uptake in prokaryotic organisms, is essential for P. aeruginosa to cause pathogenesis, and for survival under iron-deficient conditions [21].

The sigma factor algU is a key stress response regulator that controls the expression of over 300 genes, plays a crucial role in virulence factor synthesis and pathogenesis thru quorum sensing, and enhances alginate production by increasing the expression of the algD operon [22].

SOS-mediated mutagenesis involves the products of the dinB gene, which perform translesion DNA synthesis, TLS (through damage), exhibiting low accuracy but helping to rapidly replicate DNA in response to various damaging agents. However, mutations accumulate, which in turn help acquire adaptive mechanisms in response to antibacterial drugs [23].

The DNA polymerase III ε subunit, encoded by the dnaQ gene, is very important and provides 3'-5' exonuclease activity, correcting mismatches encountered during DNA repair, which allows for the removal and correction of mismatched base pairs. Mutations in the dnaQ gene can disrupt these processes, leading to more than 1000-fold increase in the mutation rate in the genome [24]. The DNA polymerase III holoenzyme consists of δ and δ' subunits, which are encoded by the holA and holB genes, forming a complex with the ε subunit of the dnaQ gene and jointly participating in DNA repair [25]. The PA0472 gene encodes the RNA polymerase σ factor. It's difficult to judge what role a specific σ factor plays in the P. aeruginosa genome, but it is known that RNA polymerase σ factors perform a huge range of vital functions: promoter recognition, double-stranded DNA unwinding, binding to RNA polymerase, and transcription control. They are also involved in the transcription of specific regulons associated with the response to environmental changes and are included in iron transport [26].

One of the RNA polymerase σ-factors involved in iron assimilation processes is the FpvI protein, encoded by the fpvI gene, which is involved in regulating the uptake of the high-affinity siderophore pyoverdine, an important virulence factor as it can displace iron from the iron–transferrin complex [27].

P. aeruginosa has 3 genes in its genome that encode TonB proteins (tonB1, tonB2, and tonB3), and only the TonB1 protein, encoded by tonB1, interacts with TonB-dependent transporters involved in iron or heme uptake [28].

In addition to the main siderophores, P. aeruginosa produces another metallophore encoded by the cntL gene, called pseudopalin, which is essential for the uptake and utilization of zinc, cobalt, and nickel in its pathogenesis. Urease, which is a nickel-dependent enzyme, is produced by P. aeruginosa, while cobalt is essential for the cobalamin-dependent ribonucleotide reductase (NrdJab), which functions in biofilm formation under oxygen-limited conditions [29].

It is known that in P. aeruginosa, sigX is involved only in the transcription of its own gene and is largely responsible for the transcription of oprF, which encodes the major outer membrane protein OprF, which in turn is involved in several crucial functions: maintaining cell structure, outer membrane permeability, and recognition by the host immune system [30]. Deletion or knockout of the algU and sigX genes in the PAO1 genome disrupts biofilm formation [31].

The capB and cspD genes are responsible for encoding cold shock proteins involved in adaptation to cold in the environment [32].

The groES gene encodes a heat shock protein that helps the microorganism survive at 42°C [33]. As is known, these resistance mechanisms are an integral part of the physiology of P. aeruginosa cells [33].

The σ³² factor, encoded by the rpoH gene, is the main regulator of the heat shock response, controlling the function of groES, among others [34].

The selected housekeeping genes confirmed their relevance in terms of their key role in the viability of P. aeruginosa, demonstrating the importance of their functioning for the biological processes of this microorganism. Furthermore, gene identification is not only based on the nucleotide sequence but is also translated into an amino acid sequence. This method was chosen to detect the stop codon in the gene and demonstrate not only its location in the genome but also its functionality. Thus, when assessing data quality, the absence of one or more of the selected genes will be considered a criterion for poor genome quality. If the N50 values for the genome selected for analysis are 10,000 or higher, but one of the candidate genes has a stop codon, it can be classified as a medium-quality genome. However, in our opinion, using a genome of this quality for phylogenetic analysis, SNP typing, or MLST analysis is not recommended. At the same time, searching for certain genes in the genome is possible, but without using them for typing the analyzed strain.

Despite the high N50 score and other evaluation parameters, the absence of two or more genes indicates poor WGS data quality. The N50 parameter is used to assess and compare the quality of genome assembly, allowing for the selection of the best among good/high-quality options.

The next criterion used to assess the quality of the WGS data was the GC content of the P. aeruginosa genomes. Analyzing the genomes of the strains using the CheckM2 program, we observed that genomes with Completeness scores > 97% and Contamination < 3% have a GC content ranging from 63.8% to 66.6%, which led us to establish threshold values of 65.2 ± 2.5%. The range was chosen wider to account for possible changes in the genomic composition. This criterion was supported by a literature review, which did not contradict our results and allowed us to include this parameter in a comprehensive quality assessment criterion for genome assemblies [35, 36]. This criterion demonstrates whether there is contamination of foreign DNA or reads from related species in the selected genome(s) for subsequent analysis.

WGS data validation using this criterion works as follows: if the genome selected for analysis falls within the established GC content values, it is considered a high-quality genome. If the selected genome does not fall within the established GC content values, it is considered a low-quality genome.

The quality of the genome is also assessed by the size of the P. aeruginosa whole-genome sequence. Thus, to evaluate the WGS data, the criteria for the minimum and maximum permissible genome sizes were used. The minimum genome size was 5.84 Mb, and the maximum was 8.26 Mb. The decision to use these values was based on literature data: for example, studies have reported that the auxiliary genome can vary within the range of 6.9–18.0% [38, 39]. The standard value for the length of the P. aeruginosa whole-genome sequence was taken as 6–7 Mb [37, 39].

Based on the above, the criterion for assessing the good quality of the P. aeruginosa genome will be a genome ranging in size from 5.84 to 8.26 Mb. If the analyzed genome falls outside the specified values, its quality will be assessed as poor or average, or the option of a more detailed and thorough analysis of this genome should be considered to exclude its structural features.

Genomes with an average level can be used limitedly for phylogenetic analysis, SNP typing, or MLST analysis, but they can be used to search for specific genes (without typing them) or for INDEL analysis.

Genomes with a low quality level are recommended not to be included in bioinformatics analysis and should be corrected by re-sequencing.

In addition to the CheckM2 and QUAST programs, which were selected as comparison tools, there are the SQUAT and Plantagora programs, but they do not meet the criteria of our research objects, as they are primarily developed for eukaryotic organisms. At the same time, CheckM2 is a tool developed for assessing the quality of prokaryotic genomes, while QUAST is a universal program. In developing our evaluation criteria, we tried to move away from complex tables with mathematical parameters assessing the quality of the WGS data provided by QUAST after the analysis. This involves the participation of bioinformatics specialists in the analysis and, in our opinion, does not fully reflect the quality of the WGS data, but rather assesses how well the genome assembly was performed [4]. At the same time, CheckM2 provides digital data on the parameters of the analyzed genome across various metrics, without drawing clear conclusions about the quality of the genome or whether it can be used for further research. The Contamination index does not always reflect the quality of the genomes of clinical isolates containing extrachromosomal elements.

Thus, we have tried, on the one hand, to select clear and concise parameters for evaluating WGS data, and on the other hand, to simplify the process for the user to obtain a specific result without resorting to in-depth bioinformatics analysis or using command linesк.

Conclusion

A comprehensive study was conducted in which we selected housekeeping genes that allow us to assess the quality of the P. aeruginosa WGS data. Quality assessment criteria have been defined: genome length and GC content, which allow for the evaluation of the P. aeruginosa genome assembly.

Based on validated assessment criteria tested on a sample of genomes, the assembly of the P. aeruginosa genome can be classified into three categories based on the quality level of the source material: good, medium, and low. Good quality — the genome length is within the average genome size for the species ± 18%, the GC content is ± 2.5% of the average for P. aeruginosa, all essential genes have been found, and their protein product translation is not disrupted by a stop codon. Average quality – all essential genes found, but errors in their translation were detected due to the formation of a stop codon as a result of a sequencing error. Low quality — at least one gene in the life support system is missing, or the genome size or GC content does not match the value characteristic of the species.

An algorithm and a publicly available program for rapid analysis based on WGS data of P. aeruginosa, Genomes Validator, have been developed.

About the authors

Alexey A. Kovalevich

Rostov-on-Don Antiplague Scientific Researsh Institute

Author for correspondence.
Email: kovalevich_aa@antiplague.ru
ORCID iD: 0000-0001-6926-0239

junior researcher, Laboratory of molecular biology of natural focal and zoonotic infections

Russian Federation, Rostov-on-Don

Alexey S. Vodopianov

Rostov-on-Don Antiplague Scientific Researsh Institute

Email: vodopyanov_as@antiplague.ru
ORCID iD: 0000-0002-9056-3231

Cand. Sci. (Med.), leading researcher, Laboratory of molecular biology of natural focal and zoonotic infections

Russian Federation, Rostov-on-Don

Ruslan V. Pisanov

Rostov-on-Don Antiplague Scientific Researsh Institute

Email: pisanov.ruslan@yandex.ru
ORCID iD: 0000-0002-7178-8021

Cand. Sci. (Biol.), leading researcher, Laboratory of molecular biology of natural focal and zoonotic infections

Russian Federation, Rostov-on-Don

References

Supplementary files