Fast and robust estimate of bacterial genus novelty using the percentage of conserved proteins with unique matches (POCPu)

Standardisation of protein sequences and taxonomy via GTDB

Portions of this text were previously published as part of a preprint (bioRxiv https://doi.org/10.1101/2025.03.17.643616). As GTDB provides curated taxonomy along with genomes and genome-derived protein sequences (Parks et al., 2022), we used it as a reliable source of high-quality data in our benchmark. We used inclusion criteria to facilitate the selection of a diverse range of taxonomic groups from GTDB (release 214) (N = 394,932 bacterial genomes), while maintaining achievable comparisons with the time, human, and computing resources available: (1) the bacteria had a valid name according to the list of prokaryotic names with standing in nomenclature (Parte et al., 2020) and a representative genome was available (N = 11,699), (2) they belonged to a family with at least two genera (N = 5,904), and (3) to a genus with at least ten genomes (N = 4,767). Based on these criteria, the protein sequence files for the shortlisted bacteria were obtained from GTDB (Table S1) which uses Prodigal v2.6.3 (Hyatt et al., 2010) for protein sequence prediction. A single representative genome from each species, as designated by GTDB, was selected for further analysis (Table S1).

Definition of percentage of conserved proteins (POCP)

The percentage of conserved proteins (Fig. 1A) between two genomes Q and S is defined as: (1)${\displaystyle POCP=\frac{C_{QS}+C_{SQ}}{T_Q+T_S}\times 100\text{\%}}$ where C_QS represents the conserved number of proteins from Q when aligned to S and conversely C_SQ represents the conserved number of proteins from S when aligned to Q; T_Q + T_S represents the total number of proteins in the two genomes being compared (adapted from Qin et al. 2014). The range of POCP is theoretically [0; 100%]. Conserved proteins are defined as protein sequences matches from the query with an e-value <10⁻⁵, a sequence identity >40%, and an aligned region >50% of the query protein sequence length (Qin et al., 2014).

Definition of Percentage of Conserved Proteins using only unique matches (POCPu)

During the alignment process, protein sequences from the query can match multiple subject sequences in the case of duplicated genes (Fig. 1A). Whilst briefly mentioned in the original article that “The number of conserved proteins in each genome of strains being compared was slightly different because of the existence of duplicate genes (paralogs)” (Qin et al., 2014, p. 2211), the impact of this on POCP values was not determined. Therefore, we defined the POCP with unique matches (POCPu) between two genomes Q and S as: (2)${\displaystyle POCPu=\frac{C_{uQS}+C_{uSQ}}{T_Q+T_S}\times 100\text{\%}}$ where C_uQS represents the conserved number of proteins from the unique matches of Q when aligned to S and, conversely, C_uSQ the conserved number of proteins from the unique matches of S when aligned to Q; T_Q + T_S represents the total number of proteins in the two genomes being compared. We hypothesized that the impact of paralogs on POCP will decrease with unique matches as only one pair of duplicated conserved genes will be counted instead of two (Fig. 1A).

Pairwise comparisons between a query sequence (Q) and a subject sequence (S) were defined by banning self-comparisons (Q ≠ S) and considering reciprocal comparisons (Q − S and S − Q) only within the same family to avoid unnecessary expansion of the comparison landscape.

Finding suitable protein sequence alignment methods to scale BLASTP-based POCP calculations

In Qin et al. (2014), guidance was provided on how to implement the computation of POCP, including the use of BLASTP. As ‘standard’ POCP method, we used BLASTP v2.14.0+ (Camacho et al., 2009) with parameters from Qin et al. (2014) (Table S2). We also considered a modified implementation of the BLASTP method, named BLASTPDB where BLAST databases are first built for the two genomes considered (Table S2). This allows parallel alignments on multiple CPU, which is not possible with BLASTP. We then included two tools that were designed as faster local-protein-alignment methods and used as alternatives to BLASTP: DIAMOND v2.1.6 (Buchfink, Reuter & Drost, 2021) and MMseqs2 v15.6f452 (Steinegger & Söding, 2017). The former is used in Hölzer (2024) while the latter is used in EzAAI (Kim, Park & Chun, 2021). Similar to BLASTPDB, these methods require that a protein database is built for each genome before performing the alignment (Table S2). DIAMOND and MMseqs2 were both used with four different sensitivity thresholds proposed recently (Table S2) (Buchfink, Reuter & Drost, 2021). This benchmark (Fig. 1B) was run on a subset of 1,235 GTDB genomes. We used genomes from a wide range of bacterial phylogenetic diversity across four phyla: Bacillota (6 families, 23 genera, 561 species), Pseudomonadota (6 families, 16 genera, 333 species), Bacteroidota (3 families, 7 genera, 217 species) and Actinomycetota (2 families, 6 genera, 124 species).

All protein matches were filtered to only keep matches with >40% identity to all the query sequences matches for POCP (C_QS and C_SQ in Eq. (1)) and only unique query sequence matches for POCPu (C_QS and C_SQ in Eq. (2)). The filtering was adapted to the method as the range of percentage of identity in MMseqs2 is [0 − 1] and [0 − 100] for BLAST and DIAMOND. The total number of proteins per genomes (T_Q and T_S in Eqs. (1) and (2)) was computed using seqkit stats v2.2.0 (Shen et al., 2016). Linear regressions were implemented using R version 4.3.1 (2023-06-16) to fit the expected POCP (or POCPu) values obtained via the BLASTP reference method against the other methods considered. The coefficient of determination R² of the linear regression is used as an interpretable and bounded goodness-of-fit measure between the expected and measured values instead of other measurement errors (Chicco, Warrens & Jurman, 2021). We did not rely on the adjusted coefficient of determination as the linear regressions had only one predictor, namely the POCP (or POCPu) values of the evaluated method.

Evaluate POCP and POCPu with known within-genus versus between-genera pairs

This analysis was run as described above except that only the scalable approach to compute POCP (or POCPu) was used on the full set of genomes (N = 4,767) to capture additional diversity (Fig. 1B): Pseudomonadota (15 families, 66 genera, 1,736 species), Actinomycetota (eight families, 27 genera, 1,584 species), Bacteroidota (six families, 27 genera, 886 species) and Bacillota (six families, 23 genera, 561 species). Next, we computed classification metrics with a positive event defined as “both genomes belong to the same genus”. Thus, for a pair of bacterial genomes with a POCP (or POCPu) >50%, the pair was a true positive (TP) if it belonged to the same genus, else it was considered as a false positive (FP). Conversely, for a pair of bacterial genomes with a POCP (or POCPu) ≤ 50%, the pair was a false negative (FN) if it belonged to the same genus, else it was considered as a true negative (TN). We then assessed the classification performance of both POCP and POCPu using Matthews correlation coefficient (MCC; Eq. (3)). (3)${\displaystyle MCC=\frac{TP\times TN-FP\times FN}{\sqrt{\left(TP+FP\right)\left(TP+FN\right)\left(TN+FP\right)\left(TN+FN\right)}}.}$ The MCC coefficient ranges from −1 to 1 and is high in the case of a perfect classification, whilst a value of 0 indicates random classification. Negative MCC values indicate perfect misclassification, as in a swap between positive and negative events. In addition, the MCC compensates for imbalanced datasets compared to other metrics such as accuracy or F1-score (Chicco & Jurman, 2020). Finally, we used one dimensional optimization (Brent, 1972) to find family-specific POCPu thresholds maximizing MCC values and separating between-genera from within-genera distributions. The optimization was run using the optimize () function from the stats package (R Core Team, 2023).

Benchmarking workflow implementation

Automatic protein sequences download, data pre-processing, many-versus-many protein alignments, POCP computation and delineation metrics calculations have been included in a workflow using nextflow v23.10.0 (Di Tommaso et al., 2017), based on components of nf-core (Ewels et al., 2020). The tools used are provided within Docker container (Merkel, 2014) or bioconda (The Bioconda Team et al., 2018) environments to ensure reproducibility and scalability and to facilitate future extension of the present benchmarking work. Nextflow natively keeps track of the time, CPU, memory, and disk usage of each process in an execution trace log file, which we used to evaluate the computing resources utilization. Process duration is available as wall-time and real-time, the CPU usage is reported as a percentage of usage of a unique CPU, meaning multi-threaded processes will have a value higher than 100%. Statistical analyses and visualization were conducted in R using targets v.1.7.0 (Landau, 2021).

Finding an alternative to BLASTP-based POCP calculations

First, we set out to identify a scalable alternative to BLASTP to compute the POCP to delineate genera (Fig. 1). We evaluated ten protein alignment methods (Table S2) based on three tools—BLASTP (Camacho et al., 2009), DIAMOND (Buchfink, Reuter & Drost, 2021) and MMseqs2 (Steinegger & Söding, 2017).

We processed 1,235 genomes and conducted 1,412,040 pairwise comparisons with a total of 32,316 CPU-Hours (3.7 in years). All the methods tested were faster than the reference BLASTP (Table S3). BLASTPDB, the database method of BLASTP (Table S2) that enables paralleled computations, was only half the time of BLASTP on average. In contrast, all DIAMOND and MMSEQS2-based methods ran at 20x and 11x the speed of BLASTP, respectively at the cost of using more memory, CPU, and disk usage (Table S3). Thus, as expected, more sensitive methods consumed more resources in general. An important criterion for a BLASTP alternative is to ensure that the increased speed does not compromise the accuracy of POCP calculation.

DIAMOND provides POCP values as accurate as with BLASTP

BLASTPDB produced the exact same POCP values as BLASTP (Fig. S2). The other methods did not perform as good. All methods of DIAMOND had a coefficient of determination (R²) above 0.99, except for DIAMOND_FAST that deviated from the expected values (Fig. 2).

All DIAMOND methods, especially DIAMOND_FAST, tended to underestimate POCP values (all dots were below the reference dashed line), meaning that they might assign genomes to different genera when they are from the same genus (top panels). Deviation from the BLASTP reference was aggravated when using the MMSEQS2_S1DOT0 method, mainly through underestimation of POCP values (bottom panels). The other MMSEQS2 methods performed better, but still less good than the DIAMOND methods. They also tended to overestimate more than underestimate. All in all, the DIAMOND methods, especially with increased sensitivity, generated POCP values nearly as accurate as BLASTP for a fraction of the time, but we refrained from using the MMSEQS2 methods due to being less accurate.

Proposal for clear and fast computation of POCP values

We observed that all methods generated POCP values exceeding the supposed upper limit of 100%. Hence, we investigated the underlying reasons and provide a clearer definition of POCP, termed POCPu (see Eq. (2)).

POCP values above 100% disappeared when using POCPu (Fig. 3 and Fig. S3). In general, the same patterns observed for POCP hold for POCPu, though with higher values of coefficient of determination (Fig. 3 and Table 1). The three different sensitive methods of DIAMOND produced POCPu values that matched perfectly the ones produced by the reference method BLASTP, with no underestimation as in the case of POCP. In contrast, the MMSEQS2 methods, whilst better with POCPu than POCP, still tended to underestimate POCPu values. All in all, DIAMOND-based POCPu is closer to its reference than POCP (Fig. 3), guiding our choice to create an accurate BLASTP alternative.

Table 1:

POCP and POCPu linear regressions results for each alternative method to BLASTP.

Method	POCP		POCPu
	R2	p-value	R2	p-value
BLASTPDB	1.0000000	<0.001	1.0000000	<0.001
DIAMOND_ULTRASENSITIVE	0.9920389	<0.001	0.9997605	<0.001
DIAMOND_VERYSENSITIVE	0.9919744	<0.001	0.9997574	<0.001
DIAMOND_SENSITIVE	0.9916398	<0.001	0.9997438	<0.001
DIAMOND_FAST	0.8967444	<0.001	0.9913959	<0.001
MMSEQS2_S7DOT5	0.9379505	<0.001	0.9722865	<0.001
MMSEQS2_S6DOT0	0.9369705	<0.001	0.9720813	<0.001
MMSEQS2_S2DOT5	0.9242783	<0.001	0.9713977	<0.001
MMSEQS2_S1DOT0	0.8659630	<0.001	0.9713360	<0.001

DOI: 10.7717/peerj.20259/table-1

Notes:

Coefficient of determination (R²) and associated p-value for linear regressions matching the POCP and POCPu values computed by each method against the respective POCP and POCPu values of the reference method BLASTP. Each linear regression are based on n = 70, 602 comparisons per method. Methods are sorted by decreasing POCPu values.

In summary, the BLASTPDB method performed as good as BLASTP (Figs. S2, S3 and Table 1) in half the time (Table S3), at the cost of using more resources. However, the DIAMOND sensitive methods were even faster with excellent adequacy with BLASTP, especially for POCPu (Table 1). Whilst DIAMOND_ULTRASENSITIVE had the highest R² value using POCPu (Table 1), it also had the highest memory consumption and disk usage (Table S3). A more sustainable alternative is DIAMOND_VERYSENSITIVE that performed 10 times faster than BLASTPDB, in less than 5% of the time of BLASTP, while still maintaining reasonable usage of the computing resources (Table S3). Importantly, POCPu values calculated using DIAMOND_VERYSENSITIVE delivered results extremely close to the reference BLASTP (Fig. 3 and Table 1) and were essentially identical to DIAMOND_ULTRASENSITIVE POCPu R² up to 5 digits (Table 1). Therefore, we consider DIAMOND_VERYSENSITIVE to be a valid and scalable alternative to BLASTP for POCP/POCPu computations.

Evaluate POCP and POCPu using within- and between-genera pairs

Unique matches enhance the accuracy of genus delineation

Next, we evaluated the 50%-threshold of POCP and POCPu to determine their reliability to delineate bacterial genera. We included all 4,767 genomes (Fig. 1) and calculated POCP and POCPu for 1,087,630 pairwise comparisons using DIAMOND_VERYSENSITIVE (Fig. 4).

Instead of two bell-shaped distributions, with the 50% threshold separating between-genera (left) from within-genus POCP values (right), we observed overlapping POCP distributions (Fig. 4A). This was associated with a high number of false positives (FP = 188,155), where between-genera values were >50%, especially compared with the number of true negatives (TN = 133,034), where between-genera values were <50%. Additionally, most of the within-genus values were >50% (TP = 220,307), with only few below the threshold (FN = 2,319).

In contrast, POCPu was much closer to the expected results given the taxonomic assignments of each genome (Fig. 4B). Between-genera POCPu values followed a bimodal distribution, with the highest peak and most of the distribution remaining below the 50%-threshold (TN = 253,860). Nonetheless, a fraction of between-genera values were higher than the threshold of 50%, representing false positives, i.e., different genera when they are not (FP = 67,329). As in the case of POCP, within-genus POCPu values were above the threshold of 50% (TP = 209,660), with a few below the threshold (FN = 12,966). All in all, considering only unique protein matches improves genus delineation.

To quantify these findings on the confusion matrix (true/false positives and negatives), we used the MCC (Eq. (3) in the Methods), which is a binary classification rate that gives a high score only when the classifier correctly predicts most of positive and negative cases. POCPu (MCC = 0.72) surpassed POCP (MCC = 0.46) to delineate bacterial genera, which quantitatively confirmed the visual findings (Fig. 4).

Family-specific POCPu thresholds enable clearer genus delineation

Analysing the bacterial families separately questioned the universal threshold of 50% conserved proteins (Fig. 5 and Fig. S4). The large family of Streptomycetaceae (Actinomycetota, Fig. 5A), was characterized by a very low MCC, and thus many false cases, including FN = 2,974 (2.8%) and FP = 7,110 (6.8%), despite even more true case TN = 1,231 (1.2%) and TP = 93,338 (89.2%). In 7 families out of 35, POCPu was clearly not adequate to delineate genus using a threshold of 50%, as indicated by low MCC (MCC ≤ 0.25; Fig. 5B). In contrast, POCPu delineated bacterial genera accurately for 18 families (MCC ≥ 0.7; Fig. 5B).

Due to these differences between families, we explored family-specific POCPu thresholds by maximizing MCC to improve genus delineation (Fig. S4 and Table 2). With this procedure, thresholds other than the default 50% would enhance classification for 19 families out of the 35 families. The genus delineation of eight families was improved, with at least 0.1-point increase in MCC (white squares in Table 2; Actinomycetota: Streptomycetaceae and Streptosporangiaceae; Bacillota: Amphibacillaceae, Lactobacillaceae and Metamycoplasmataceae; Pseudomonadota: Acetobacteraceae, Burkholderiaceae_B and Vibrionaceae). For 11 additional families maximum MCC above 0.7 were even obtained (black squares in Table 2; Actinomycetota: Micromonosporaceae and Pseudonocardiaceae; Bacteroidota: Bacteroidaceae and Weeksellaceae; Pseudomonadota: Burkholderiaceae, Enterobacteriaceae, Halomonadaceae, Pseudomonadaceae, Rhizobiaceae, Rhodobacteraceae and Xanthomonadaceae). Interestingly, in two cases, the optimal POCPu threshold was lower than the standard threshold: 43% for Bacteroidaceae (Bacteroidota) and 45.5% for Metamycoplasmataceae (Bacillota). In 17 cases, new thresholds higher than 50% conserved proteins better separated genomes from within genus and between genera (Table 2).

Table 2:

Proposal of family-specific POCPu thresholds for genus delineation.

	Family	ΔMCC	Threshold value (%)		MCC	Max. MCC
Actinomycetota
	Mycobacteriaceae	0.02	52.3		0.89	0.91
	Micrococcaceae	0.07	53.9		0.88	0.95
	Nocardioidaceae	0.00	48.8		0.74	0.74
□	Streptosporangiaceae	0.15	63.1	↑	0.72	0.87
	Microbacteriaceae	0.10	53.7		0.63	0.73
■	Pseudonocardiaceae	0.23	56.2	↑	0.63	0.86
■	Micromonosporaceae	0.41	57.9	↑	0.57	0.98
□	Streptomycetaceae	0.33	55.7	↑	0.16	0.49
Bacillota
	Planococcaceae	0.00	51.3		1.00	1.00
	Streptococcaceae	0.03	47.7		0.93	0.97
□	Metamycoplasmataceae	0.14	45.5	↓	0.81	0.95
	Paenibacillaceae	0.06	48.3		0.75	0.81
□	Amphibacillaceae	0.21	53.2	↑	0.74	0.95
□	Lactobacillaceae	0.11	57.9	↑	0.71	0.82
Bacteroidota
	Hymenobacteraceae	0.01	53.6		0.99	1.00
	Spirosomaceae	0.02	52.9		0.98	1.00
	Sphingobacteriaceae	0.01	50.4		0.92	0.92
■	Bacteroidaceae	0.17	43.0	↓	0.63	0.80
	Flavobacteriaceae	0.01	51.4		0.59	0.60
■	Weeksellaceae	0.63	60.1	↑	0.22	0.85
Pseudomonadota
	Xanthobacteraceae	0.00	45.8		1.00	1.00
	Beijerinckiaceae	0.00	50.1		0.97	0.97
	Moraxellaceae	0.02	53.2		0.92	0.94
□	Burkholderiaceae_B	0.12	54.1	↑	0.85	0.97
	Alteromonadaceae	0.01	49.3		0.78	0.79
	Sphingomonadaceae	0.10	52.9		0.63	0.73
■	Burkholderiaceae	0.22	64.3	↑	0.52	0.74
■	Rhizobiaceae	0.43	62.5	↑	0.48	0.91
■	Rhodobacteraceae	0.27	58.1	↑	0.43	0.70
■	Xanthomonadaceae	0.51	63.7	↑	0.33	0.85
□	Vibrionaceae	0.28	59.3	↑	0.23	0.51
■	Enterobacteriaceae	0.60	71.2	↑	0.19	0.79
■	Halomonadaceae	0.56	59.6	↑	0.19	0.76
■	Pseudomonadaceae	0.82	63.2	↑	0.07	0.89
□	Acetobacteraceae	0.63	64.3	↑	0.05	0.68

DOI: 10.7717/peerj.20259/table-2

Notes:

The thresholds were obtained after maximizing the MCC value to separate between-genera from within-genera distributions. Squares indicate that MCC value change was greater than 0.1 with the optimized threshold, filled squares denote rescued families from MCC < 0.7 to MCC > 0.7, whilst empty squares indicate improved genus delineation. Arrows highlight potential family-specific threshold worth considering to replace the default of 50% with the direction of change. Families without squares already delineate genera correctly with the default of 50%.

Because POCP was previously proposed to be influenced by genome size (Riesco & Trujillo, 2024), we used the large pairwise comparisons dataset to assess whether the changes in threshold were linked to differences in genome size. If POCPu is influenced, we reasoned that its genus delineation power should also be influenced, therefore we expected stronger genome size differences in the families for which an alternative POCPu threshold was found. We found no evidence that POCPu is affected by differences in genome size (Fig. S5A) nor proteins number (Fig. S5B).