MetaBIDx: a new computational approach to bacteria identification in microbiomes

Phase 1: Building the index of a microbiome

MetaBIDx employs a modified Bloom filter^[32] as its indexing mechanism. A Bloom filter is a space-efficient probabilistic data structure used for efficient membership queries. Although a Bloom filter may mistakenly identify items not stored in it (false positive mistakes), it can accurately recall stored items.

Index construction for a microbiome relies on the reference genomes of species present. Without prior knowledge of the microbiome’s composition, one can use comprehensive bacterial and viral genomes from public databases for index creation. However, if specific species within the microbiome are of interest or there is partial knowledge about its composition, the index can be built using only those species’ reference genomes.

In constructing the index F for referenced bacterial genomes, each genome’s k-mers are processed. This includes both the main and reverse complement strands of a genome. A set of n randomly generated hash functions maps each k-mer to n entries in the index. These entries can hold three types of values: 0 for an empty entry, -1 for a “dirty” entry, and a positive genome id. An entry in F has value -1 if two k-mers from different genomes hash to it, indicating the k-mer is not unique. If a k-mer from one genome hashes to an entry already holding a different genome’s id, that k-mer is not unique, and the entry is marked as “dirty”. Otherwise, the entry is updated with the genome id of the currently processed genome if it is empty or already contains the same genome id.

This approach allows F to function similarly to a Bloom filter, aiding in the detection of genomes present in a metagenomic sample. The algorithm for building this index processes each genome sequentially, updating the index entries accordingly.

Phase 2: Determining species in a microbiome

The process consists of two steps. First, the index F is used to assign reads in the sample to the species stored in the index. Second, these reads are further clustered into groups of similar coverages to determine which species are present in the sample.

Data collection and preparation

To evaluate the proposed method, we employed metagenomic shotgun sequencing data without imposing data quality constraints or specific requirements. Three widely recognized datasets were utilized for this assessment, comprising two mock community datasets and one derived from a human sample.

Mende Dataset: this dataset (available at https://swifter.embl.de/~mende/simulated_data), comprises three metagenomic samples. These samples are distinguished by their species complexity, featuring 10, 100, and 400 species, respectively. Each sample contains 75 bp long reads. The number of reads varies from 26,665,674 to 26,667,004 pairs. This dataset, originally used in a study on metagenomic assembly^[35], was constructed using simulated Illumina sequencing errors and quality values, reflecting the characteristics of actual metagenomic data.

CAMI Challenge Dataset: this dataset was obtained from the CAMI challenge^[36] (accessible at https://data.cami-challenge.org). It was also used in another benchmark^[37]. It includes eight metagenomic samples representing a gradient of complexity: low (RL_S001), medium (RM_S001, RM_S002), and high (RH_S001, RH_S002, RH_S003, RH_S004, RH_S005). These samples are characterized by experimental conditions and features akin to real datasets, such as the inclusion of multiple, closely related strains, the presence of plasmid and viral sequences, and realistic abundance distributions. The reads in this dataset are 150 bp in length. The number of reads varies from 49,898,179 to 49,905,935 pairs, from low complexity to high complexity.

PT-8 (S2): this sample was used in a study^[38]. It was derived from brain tissue biopsies of a 67-year-old patient with osteomyelitis, lung disease, and multifocal brain and spinal lesions, and was diagnosed with Mycobacterium tuberculosis. The PT-8 (S2) sample, with human reads excluded, can be accessed at the NCBI SRA repository (https://www.ncbi.nlm.nih.gov/sra/SRX1621515). This sample consists of bacterial, viral, and fungal species.

Reference Genomes: methods that use reference genomes to build indices for species identification. For the Mende dataset, we collected reference genomes that contain reads from all four metagenomic samples. This resulted in a set of 457 bacterial genomes, averaging approximately 3.5 Mbp each, with a cumulative size of around 1.6 GB. In contrast, the reference genomes for the CAMI dataset were more extensive, encompassing 2,850 bacterial genomes. This collection included all bacterial genomes present in the reads from the eight samples, along with genomes of closely related strains within the same species or subspecies. The average genome size in this collection was about 5.7 Mbp, culminating in a total size of approximately 16 GB. The genome collection for the CAMI dataset was utilized in the experiment with real dataset as it also contains Mycobacterium tuberculosis. All reference genomes were downloaded from NCBI.

Comparative analysis of different methods for species prediction

We conducted a comprehensive comparative analysis of our tool, MetaBIDx, and various metagenomic tools, which can be used for species detection, including CLARK^[24], Kraken2^[39], KrakenUniq^[40], Centrifuge^[41], and Sourmash^[42]. These tools were selected based on their robustness, documented accuracy, processing speed, and the ability to create custom indexes, a vital feature for our analysis. Tools that lacked comprehensive documentation or presented installation and experimental challenges were not considered for this study.

After experimenting with several k-mer lengths, MetaBIDx was built with k-mers of length 31. A previous study^[43] experimented with different k-mer sizes and observed that k-mer similarity between genomes approximated various degrees of taxonomic similarity, and that a k-mer length of 31 appeared to correspond to species-level similarity. Most tools used in the evaluation also have the default k-mer size of 31. An index with a size of 8 GB, utilizing 3 hash functions, was created for the Mende dataset, which included 457 reference genomes. For the CAMI dataset, a more extensive index of 16 GB was built using 2 hash functions, accommodating 2,850 reference genomes. The other tools in the study also employed k-mers of length 31 and were run with their default settings for a fair comparison. To ensure consistency across all methods, the same collections of reference genomes were used to construct the genome libraries or indexes. The experiments were conducted on a standardized computational setup, using a machine with 32 cores and 330 GB of RAM, and all tools were run in multi-threaded mode to utilize the full computational capacity. The script for building index for all tools is shared in the Supplementary Materials.

The evaluation of prediction performance utilized three widely recognized metrics: precision, recall, and F1-score. Precision quantifies the proportion of true positive predictions out of all positive predictions made (sum of true positives and false positives), essentially reflecting the accuracy in predicting species as a fraction of all species predictions. Recall measures the proportion of true positive predictions relative to the total actual species present in the sample, indicating the method’s ability to identify all relevant species. The F1-score, derived as the harmonic mean of precision and recall, offers a composite metric that equally weights precision and recall, providing a single measure to assess the balance between them.

The comparative assessment of our method, MetaBIDx, alongside other tools, was structured into two key experiments. The first experiment focused on evaluating the ability of MetaBIDx and other tools to predict bacterial species based on identified reads alone, which is the standard approach adopted by these tools for species prediction in metagenomic samples. This experiment’s main objective was to show that the default behavior of these tools could be enhanced for more accurate species identification.

In the second experiment, we aimed to facilitate a more equitable comparison by augmenting the other methods with our strategy for reducing false positives. This was achieved by applying our technique of clustering “approximate” coverages to each tool. The inclusion of this approach in the assessment was intended to improve the precision of the other tools, thereby enabling a fairer and more balanced comparison with MetaBIDx. Through this two-pronged experimental design, we sought to comprehensively evaluate and demonstrate the effectiveness of our method in the context of metagenomic species prediction.

Additionally, we also explored the effect of using high-quality k-mers at different thresholds to enhance species prediction accuracy. This dual-phase approach was designed to provide a thorough understanding of the capabilities and limitations of each tool in metagenomic species identification.

Comparative analysis of species identification

First, we compared MetaBIDx against other methods that predict species solely based on read classification. This approach predicates species prediction on the detection of reads originating from the species in question. If a read from a particular species is detected in a sample, that species is predicted to be present. Such methods rely on the accuracy of reads identification to determine the presence or absence of species in a sample.

Mende Dataset. Table 1 reports the performance of our method versus the others on the Mende dataset. The comparison clearly demonstrates that MetaBIDx significantly outperforms other tools in predicting species across different sample complexities.

In the 10-species sample, MetaBIDx achieved perfect precision (1.000) and a high recall (0.800), resulting in an F1-score of 0.889. Sourmash also achieved a competitive precision (0.900) and a high recall (0.900), as well as the highest F1-score of 0.900. This is markedly superior to the other classification tools, which, despite high recall values (0.800 for KrakenUniq, Kraken2, and Centrifuge, and 1.000 for CLARK), had very low precision, leading to considerably lower F1-scores (ranging from 0.041 to 0.075 for the two Kraken and Centrifuge tools, and 0.053 for CLARK).

In the 100-species sample, MetaBIDx again maintained perfect precision (1.000) and an increased recall (0.976), resulting in an F1-score of 0.988. Sourmash still had competitive precision (0.852), recall (0.750), and an F1-score of 0.800. The other tools showed improvement in precision compared to the 10-species sample but were still significantly lower than MetaBIDx. Their recall values ranged from 0.713 to 1.000, and F1-scores varied between 0.345 and 0.382, still considerably lower than MetaBIDx.

For the most complex sample with 400 species, MetaBIDx maintained its high performance with perfect precision (1.000) and a recall of 0.970, leading to an F1-score of 0.985. Other tools showed a notable improvement in this category, with CLARK reaching an F1-score of 0.933. However, MetaBIDx still outperformed them, as the F1-scores for KrakenUniq, Kraken2, Centrifuge, and Sourmash were 0.773, 0.807, 0.879, and 0.843, respectively.

Overall, MetaBIDx consistently exhibits superior performance in species prediction across various sample complexities, particularly in maintaining high precision without sacrificing recall, leading to significantly higher F1-scores compared to other tools.

CAMI dataset. Table 1 reports the performance of our method versus the others on the Mende dataset. The results highlight the comparative efficacy of these tools across a range of samples with varying complexities.

High complexity samples (RH_S001 to RH_S005): In these samples, MetaBIDx consistently demonstrated superior precision, ranging from 0.839 to 0.885, and recall values varied from 0.449 to 0.778. This resulted in F1-scores between 0.591 and 0.807, significantly higher than the other tools. In contrast, tools like CLARK, KrakenUniq, Kraken2, and Centrifuge exhibited much lower precision (consistently below 0.1) and F1-scores, despite having high recall values. This indicates that while these tools were able to identify a broad range of species (high recall), they also misidentified many species (low precision), reducing their overall accuracy. Sourmash’s F1-scores varied from 0.365 to 0.385, slightly higher than other tools.

Medium complexity samples (RM_S001 and RM_S002): MetaBIDx again outperformed the other tools with higher precision and F1-scores. Particularly in RM_S001, MetaBIDx achieved a high precision of 0.821, albeit with a lower recall of 0.397, resulting in an F1-score of 0.535. In RM_S002, Sourmash achieved the highest F1-score (0.569), followed by MetaBIDx (0.427). The other tools, while maintaining perfect or near-perfect recall, had very low precision and F1-scores, indicating a high rate of false positives in species identification.

Low complexity sample (RL_S001): MetaBIDx achieved perfect precision (1.000) and a recall of 0.421, leading to an F1-score of 0.593. This is substantially better than the other tools, which, despite having high recalls, had extremely low precision and F1-scores.

In summary, MetaBIDx consistently outperformed the other tools across the CAMI dataset, particularly in terms of precision and F1-score. This suggests that MetaBIDx is more effective in accurately identifying the species present in a sample, with fewer false positives compared to other methods. Its performance is particularly notable in complex samples, where accurate species identification is more challenging.

Enhancing precision with clustering of “approximate” coverages

Here, we aim to present a more equitable comparison between our tool and the others. The inherent limitation of these read-classification methods lies in their tendency to generate a high number of false positives, leading to lower precision in species prediction, as observed in the previous section.

We applied the false-positive-reduction strategy to all methods. This strategy groups species with similar coverages into the same clusters. The rationale is that true positives (actual species present in the sample) will typically show higher coverage compared to false positives (species incorrectly identified due to random read matches). By filtering out species with low coverage, we aim to reduce the number of false positives, thereby increasing the precision of species identification across all methods.

In this experiment, we excluded Sourmash as this tool is not a read classification tool. The result is summarized in Table 2.

Mende Dataset: Applying the clustering strategy significantly improved the precision of all methods, particularly in the 10-species sample, where MetaBIDx, CLARK, KrakenUniq, and Kraken2 all achieved perfect precision (1.000). MetaBIDx maintained its high performance across all samples, consistently showing high precision and recall, resulting in F1-scores ranging from 0.889 to 0.985. CLARK, KrakenUniq, and Kraken2 also showed notable improvements in precision and F1-scores compared to their performance without clustering, particularly in samples with 100 and 400 species.

CAMI Dataset: In the CAMI dataset, the application of clustering also enhanced precision for all methods. MetaBIDx consistently demonstrated high precision and recall, with F1-scores ranging from 0.535 to 0.807 across different samples. Other tools, including CLARK and KrakenUniq, exhibited considerable improvements in precision, leading to higher F1-scores compared to their initial performance without clustering. However, MetaBIDx maintained an edge in terms of overall accuracy.

In conclusion, by adopting the clustering of “approximate” coverage, all methods showed an increase in precision, thereby reducing false positives. This approach demonstrates that integrating coverage-based clustering can significantly enhance the accuracy of species prediction in metagenomic analysis. MetaBIDx, with its inherent design to utilize this technique, consistently outperformed or matched the performance of other tools under this enhanced comparison framework.

Identification of pathogens in human samples

We evaluated the performance of all tools in identifying the pathogen in the human sample PT-8 (S2) dataset, at the species level using an index built from 2,850 reference genomes. This sample was diagnosed with a disease organism, which we assumed as the ground truth.

We found that MetaBIDx had the highest rate of identified reads at 83%, followed by Kraken2, KrakenUniq, and Centrifuge with similar rates. CLARK had the lowest rate of identified reads, only reaching 42%. All tools assigned approximately 70% of identified reads to Mycobacterium tuberculosis and the remaining 30% to other species.

When clustering of species based on coverage derived from identified reads was used, all tools identified Mycobacterium tuberculosis as the predicted species. It is important to note that PT-8 (S2) was used in a prior study^[38] and was derived from brain tissue biopsies of a 67-year-old patient with osteomyelitis, lung disease, and multifocal brain and spinal lesions. The patient was diagnosed with Mycobacterium tuberculosis and responded promptly to anti-tuberculous treatment. This suggested that our approach to reducing false positives via clustering based on genome coverage was effective and could be clinically beneficial.

The impact of using high-quality k-mers

Sequencing errors can lead to false positives, reducing the precision of species prediction. We evaluated the impact of k-mer quality on the accuracy of bacterial prediction using MetaBIDx. K-mer quality was determined by averaging the quality scores of its constituent bases. For the Mende dataset, thresholds of 33 (the lowest quality) and 49 (the median quality) were used; for the CAMI dataset, thresholds of 18 (the lowest quality) and 30 (the median quality) were used. MetaBIDx lets users adjust the k-mer quality parameter to a desirable level.

The findings, summarized in Table 3, indicate that higher quality thresholds for k-mers lead to improved precision in bacterial species identification. This improvement was particularly notable in the CAMI dataset, where a significant increase in precision was observed, although it was accompanied by a slight reduction in recall.

Mende Dataset: For the 10 and 100 species samples, there was no change in precision and recall when the k-mer quality threshold was increased from 33 to 49. In the 400 species sample, both thresholds (33 and 49) resulted in high precision and recall, with an F1-score of 0.972.

CAMI Dataset: Across all samples, increasing the k-mer quality threshold from 18 to 30 led to a notable improvement in precision and F1-scores. The increase in k-mer quality, however, resulted in a decrease in recall, though the overall F1-score improvement suggests a favorable balance between precision and recall. The trade-off between precision and recall highlights the importance of selecting an optimal k-mer quality threshold that balances the need for accurate species identification and comprehensive species detection.

Running time analysis

We reported the running time of MetaBIDx on building the Mende index and querying reads from 400 species sample using different numbers of CPU(s) in Supplementary Tables 1 and 2. The running time of building the Mende index decreases significantly from 165 to 65 min as the CPU increases from 1 to 32. Similarly, the running time of querying reads from 400 species sample also decreases significantly from 145 min with 1 CPU to 16 min with 32 CPUs. The results in Supplementary Tables 1 and 2 indicate the effectiveness of parallelization in reducing computational time.

The comparison of running times between our method, MetaBIDx, and other tools such as CLARK, KrakenUniq, Kraken2, and Centrifuge indicates that MetaBIDx generally has longer running times across different samples in both the Mende and CAMI datasets. For the Mende dataset, on average, for each sample, MetaBIDx took 16 min, CLARK took 6 min, Kraken took 2 min, KrakenUniq and Centrifuge took 7 min. For the CAMI dataset, the average running times are 71, 17, 17, 8, and 16 min for MetaBIDx, CLARK, KrakenUniq, Kraken2, and Centrifuge, respectively. Supplementary Table 3 reports detailed information on this comparison.