MTSv: rapid alignment-based taxonomic classification and high-confidence metagenomic analysis

Building the metagenomic FM-index

MTSv implements a custom metagenomic index (MG-index) based on the FM-index data structure (Ferragina & Manzini, 2000). To build the index, MTSv takes a multi-FASTA file with headers indicating the Taxonomic ID (TaxID) for each reference sequence (Fig. 1). The reference sequences can be sourced from any DNA sequence collection (i.e., GenBank, RefSeq, etc.) and customized to fit the project. Because MTSv was designed to be highly parallelizable, we recommend building multiple indices from smaller chunks of the reference sequences. This helps reduce the memory requirements and allows for faster processing for both index building and assignment. A tool for splitting up the reference files is included with MTSv.

During the index build, MTSv concatenates the sequences in the FASTA file while recording the location of sequence boundaries and the TaxID associated with each sequence. MTSv then builds a suffix array (Manber & Myers, 1993), Burrows-Wheeler Transform (Burrows & Wheeler, 1994) and FM-index (Ferragina & Manzini, 2000) from the concatenated sequences using the Rust-Bio v0.39.1 package (Köster, 2016). The FM-index with the associated sequence metadata constitutes the MG-index. One MG-index is created per FASTA file, and new indices can be added as the reference collection grows without needing to rebuild any of the existing indices.

Taxonomic assignment

When classifying metagenomic sequences, MTSv uses the MG-Index to narrow down reference sequence locations that are most likely to contain a good alignment. The binning step begins by extracting overlapping substrings of the same size (seed-size) using a user-provided offset (seed-interval) from each read and its reverse complement (Fig. 2). It then uses the MG-Index to search for exact, un-gapped matches for each seed. The seed matches are sorted by location and grouped into candidate regions using specified windows. The number of hits per candidate is tallied and any candidate that does not meet the user-specified minimum number of seed hits is filtered out.

The filtering step is based on a q-gram algorithm which defines the minimum number of exact k-mer matches (from all $n-k+1$ overlapping k-mers) that can be expected between an n-length read and a reference sequence with at most e mismatches. In the worst case, where all mismatches are evenly spaced across the alignment, the minimum number of matching k-mers is: $m=\left(n+1\right)-k\left(e+1\right)$ ( $m>0$ when $\lfloor \frac{n}{e+1}\rfloor >k$) (Burkhardt et al., 1999; Rasmussen, Stoye & Myers, 2006; Reinert et al., 2015; Ukkonen, 1992). If only every l^th overlapping k-mer is used, the minimum number of matching k-mers is expected to be ${\displaystyle \frac{m}{l}}$. Within these parameters, MTSv is highly likely to find all alignments with up to e mismatches but will still find many alignments with more than e mismatches up to the user provided maximum edit distance cutoff.

After the filtering step, the remaining candidate positions are sorted by the number of seed hits in descending order so that the most promising regions are evaluated first. For each candidate region, MTSv extracts the corresponding range from the reference sequence and looks up the TaxID associated with the reference in the MG-index. If the current read has already been successfully aligned to the TaxID associated with the candidate region, no additional alignment is attempted, and the next candidate region is checked. Otherwise, a Single Instruction Multiple Data (SIMD)-accelerated Smith-Waterman alignment (Zhao et al., 2013) is performed between the extracted reference sequence and the query sequence (using a scoring of 1 for matches and −1 for mismatches, gap opening, and gap extension). If the alignment score is sufficiently high, there is one final check to determine if the edit distance is less than or equal to the user-specified edit distance cutoff. If the alignment is considered successful, then no further alignments are attempted for that read against the same TaxID. Skipping all additional alignments to a TaxID avoids many expensive operations and reduces computation time. As the number of reference genomes (and thus the size of a taxonomic bin) continues to rapidly expand, the advantage of this computational speed-up will only become more critical.

Reads are binned against each MG-Index producing separate output files. Each line of the output file contains the read name and a comma separated list of the TaxIDs it aligned to with the edit distance (e.g., “R1:1392=2,1396=3”). The edit distance provided for a TaxID is recorded from the first successful alignment (because no more are attempted afterward) and therefore it is not guaranteed to be the minimum. However, because the alignments start with the most promising candidates first, the edit distance is likely to be near the minimum for each TaxID. Once all output files from the binning step are completed, the files are merged to combine read assignments across multiple MG-Indices.

Comparison to existing software

We compared the performance of MTSv to two traditional read alignment tools (Bowtie2 v2.4.5 (Langmead & Salzberg, 2012) and Minimap2 v2.24 (Li, 2018)) and two commonly used taxonomic assignment tools designed for metagenomic data (Centrifuge v1.0.4 (Kim et al., 2016), and Kraken2 v2.1.2 (Wood, Lu & Langmead, 2019)). Bowtie2 is a well-established alignment tool that uses an FM-Index assisted seed-and-extend approach similar to MTSv. Minimap2 was introduced more recently and appears to be a replacement for the BWA suite of tools (Li, 2007; Li & Durbin, 2009). It uses a hash table index and a seed-chain-align approach that is faster than Bowtie2 and BWA-MEM for Illumina short read alignments (Li, 2018). Centrifuge was designed for taxonomic classification of metagenomic sequences and uses an FM-Index data structure. However, unlike MTSv which performs a full alignment, Centrifuge only extends seed hits to find the longest exact matching segments and does not account for gaps or mismatches. Centrifuge assigns scores that favor taxa with longer exact matching segments and then assigns the read to the taxon/taxa with the highest score. Kraken2 is a widely used ultrafast metagenomic classification tool that uses an exact k-mer matching approach. Kraken2 builds a database that maps fixed size k-mers from reference sequences to the lowest common ancestor of all organisms that share the k-mer sequence. The k-mers extracted from sequencing reads can then be looked up in a hash table and assigned to the corresponding TaxID. The final assignment is chosen based on a score from a weighted path classification tree (Wood & Salzberg, 2014). Although Centrifuge and Kraken2 (with the help of Bracken (Lu et al., 2017)) can use their taxonomic assignments to estimate species abundance profiles, for comparison purposes, we focused only on their performance on taxonomic assignment.

Building the reference indices

To compare the taxonomic assignment tools, reference indices were built from bacterial/archaeal genome assemblies from RefSeq (O’Leary et al., 2016) at the Complete Genome assembly level (121 GB) downloaded on October 28^th, 2019 (Release 96). TaxIDs were assigned to each sequence using the accession to TaxID mapping file provided with the download (all TaxIDs were labeled at the species level or higher) and there was a total of 4,785 unique TaxIDs (4,530 bacteria and 255 archaea). To reduce the memory requirements for building and loading the MG-index, we split the reference file into 13 roughly 10 GB chunks and we built individual MG-indices for each chunk. Centrifuge and Kraken2 do not have an option to split up their data structures, so we built them with a single sequence file using default parameters (Kraken2 k-mer length was 35 bp). For comparison purposes, the Centrifuge index was built without reference sequence compression and to avoid variation due to masking methods, low complexity regions were not masked for any of the tools.

Two additional sets of toy reference indices were built to compare MTSv and the traditional alignment tools. The collection of references was chosen to highlight biases that can be introduced when using traditional alignment tools for taxonomic assignment of metagenomic reads. The first set contained reference sequences of closely related species in the Bacillus genus and was dominated by Bacillus anthracis genomes (Table 1). The second was dominated by Staphylococcus aureus but also contained a Staphylococcus epidermidis reference sequence. The accession numbers for the sequences are available in Supplemental Table S1. MTSv, Bowtie2, and Minimap2 indices were built for both reference sequence collections using default parameters.

Parameter settings

For Bowtie2, we used parameter settings that were optimized for taxonomic binning of metagenomic reads (Jaillard et al., 2016); namely, a seed length of 20 (-L 20), seed error of 0 (-N 0), effort level of 15 consecutive seed extensions (-D 25), and end-to-end mode (—end-to-end–min-score L,-0.6,-0.4). Default values were used for Minimap2. In some cases, the aligners were also run in a way that forced them to find and report all alignments for each read. For Bowtie2 this was achieved using the “all” mode (-a) and for Minimap2 we used the “-N” option to retain up to 5,000 alignments. Our goal was to ensure that the search effort extended enough to cover alignments to all references in the index and the limit of 5,000 was chosen because it was much larger than the number of sequences. Bowtie2 also has a -k method to set an upper bound on the number of alignments but in our example, they were effectively the same. When comparing the performance of these tools they were run on the 2.60 GHz Intel Skylake Xeon nodes on Northern Arizona University’s Monsoon high-performance computing cluster with 120 GB of memory using eight threads each.

Default parameters were used for Kraken2 and Centrifuge. Unless otherwise noted, the default parameters were also used for MTSv: the seed size was set to 18 bp, the space between subsequent overlapping seeds was 15 bp, the minimum number of seed matches to perform a full alignment was set to 1, and the maximum edit distance cutoff was 13% or 20 mismatches for 150 bp reads. For MTSv, only TaxIDs with the lowest reported edit distance were considered for the final assignment. For performance comparison these tools were also run using 2.60 GHz Intel Skylake Xeon nodes using 80 GB of memory and eight threads each.

Performance evaluation using CAMI datasets

To evaluate the performance of MTSv compared to Centrifuge and Kraken2, we used the human microbiome datasets produced for the second Critical Assessment of Metagenome Interpretation (CAMI2) challenge (Meyer et al., 2022). The simulated datasets were modelled on Human Microbiome Project Consortium (Huttenhower et al., 2012) profiles from five different body sites using NCBI RefSeq (O’Leary et al., 2016) complete genome sequences as genome sources (Downloaded on July 8, 2017). We analyzed the short-read data (Illumina HiSeq2000 150 bp paired-end reads) from the gastrointestinal tract, oral cavity, airway, skin, and urogenital tract. The ground truth source of each read was provided, and the performance metrics were calculated using Amber v.2.0.2 (Meyer et al., 2018). For both MTSv and Centrifuge, reads with assignments to multiple taxa were rolled up to the lowest common ancestor using ETE3 v3.1.2 (Huerta-Cepas, Serra & Bork, 2016) whereas Kraken2 does this by default.

Amber is a tool that helps standardize comparative assessment of metagenomic binning tools by providing performance metrics evaluated against the ground truth assignments for each read. Amber defines true positives (TP) as sequences that are assigned to the correct taxon at a given taxonomic level and false positives (FP) as those assigned to an incorrect taxon. False negatives (FN) include unassigned sequences and those that are not assigned at the current taxonomic level (e.g., when multiple assignments are classified at a higher common taxonomic level than the one being evaluated). Sample purity (or precision) is the sum of all true positive sequences in the sample divided by the total number of assigned sequences.

$$p={\displaystyle \frac{TP}{TP+FP}}$$

Average purity is the precision calculated for each taxonomic bin averaged across all bins. For average purity small bins have the same weight as larger bins whereas larger bins contribute more to the sample purity value. Sample accuracy was calculated as the sum of true positive sequences divided by the total number of sequences (including those that were unassigned).

$$a={\displaystyle \frac{TP}{FN+TP+FP}}$$

The average completeness is the sensitivity of each taxonomic bin averaged across all bins. Completeness is calculated as the number of true positive sequences in the bin divided by the true size of the bin. Again, the average completeness weighs all bins equally, regardless of their size.

$$c={\displaystyle \frac{TP}{TP+FN}}$$

Simulating reads

We used Mason v2.0.9 (Holtgrewe, 2010) to simulate 150 bp Illumina reads using default parameters. To compare the traditional alignment tools using the toy indices, we simulated and combined 1 million reads from a single reference sequence from each species in the index (Table S1). When analyzing how divergent genomes impacted taxonomic assignment, we used Mutation-Simulator v2.0.3 (Kühl, 2019) to simulate SNPs, insertions, and deletions in a reference genome for Faecalibacterium prausnitzii (GCA_003312465.1) at rates ranging from 0–15% with indels occurring at 1/10 the rate of the SNPs. From these genome sequences we simulated 1 million paired reads using Mason. These reads were then added to the CAMI2 gastrointestinal tract reads which we confirmed did not contain any sequences from F. prausnitzii. For benchmarking the taxonomic assignment tools, 10 million reads were simulated from the reference sequences in the bacterial RefSeq database.

Pathogen spike-in to soil background

Using real sequencing data, we tested the performance of MTSv, Centrifuge, and Kraken2 at detecting pathogen DNA spiked into a complex soil background. We extracted and quantified, via qPCR (Liu et al., 2012), DNA from strains of Bacillus anthracis, Francisella tularensis, and Yersinia pestis from the Pathogen and Microbiome Institute’s culture collection. The pathogen DNA was then spiked into DNA extracted from pooled environmental soil samples collected in Puerto Rico in 2016. We used qPCR to verify that none of the pathogens were present in the soil background prior to the spike-in (for B. anthracis, we used primers and protocols targeting a SNP in the plcR gene (Easterday et al., 2005), for F. tularensis we used the Ft-sp.FTS_0772 assay primers and protocols (Öhrman et al., 2021), and for Y. pestis, we targeted the pla gene using primers and protocols from Hinnebusch & Schwan (1993), see Supplemental Methods). DNA from each pathogen was spiked into separate aliquots of the soil background at a 2:98 ratio (roughly 2% of the final mixture). The mixtures and the soil background were sequenced on an Illumina MiSeq platform (V3, 2 × 300 kit). The reads were quality-filtered using Fastp v0.20.122 and trimmed to 150 bp. The raw reads are available through the NCBI Sequence Read Archive with the following accessions: SAMN31055686 (B. anthracis mixture), SAMN31055687 (Y. pestis mixture), SAMN31055688 (F. tularensis mixture), and SAMN31055689 (soil background).

Traditional read alignment tools were not designed or intended for taxonomic assignment against multi-species reference databases

The standard search effort and reporting modes of traditional read alignment tools result in misclassification when used for taxonomic assignment. Bowtie2 and Minimap2 missed thousands of alignments particularly for species that were underrepresented in the reference index (Table 1) or when genomes from closely related species were included. By default, Bowtie2 reports only the best alignment found after a defined level of search effort (choosing one at random if there is a tie) whereas Minimap2 reports up to five alternative alignments per read. When the standard search effort and reporting modes were used, many reads aligned to the dominant and/or closely related taxa rather than the correct taxon. One index contained closely related species (>99% identity in some cases) from the Bacillus genus (Helgason et al., 2000; Rasko et al., 2005) and was dominated by Bacillus anthracis references. Due to the similarity between the genomes, valid alignments for reads of one species were often found within other references. In some cases, the search effort was expended before the alignment to the correct genome was found. The other index included references for S. aureus and S. epidermidis which share many genes in common but the average nucleotide identity among homologous sequences is low around 75% (Méric et al., 2015). Despite this, many alignments to S. epidermidis were missed because S. aureus sequences were highly overrepresented (961:1) resulting in biased alignments. Minimap2 performed better than Bowtie2 because it reported some alternative alignments by default. But in either case, when using a multi-species reference index, the level of search effort did not extend to cover all taxa evenly resulting in biases and missed assignments.

When Bowtie2 and Minimap2 were forced to find and report all alignments, the number of missed assignments was reduced but the runtime and memory requirements increased (Fig. 3). Bowtie2 took over 7 h to align 2 million reads against the Staphylococcus index compared to just a couple of minutes in standard mode (~300 times longer). Minimap2 was faster than Bowtie2 although it did not report as many alignments (855 M vs 986 M), and it failed to find any alignments for over 5,000 reads. Minimap2 took over an hour compared to ~18 min in standard mode (~4 times longer). Minimap2 also required 111 GB of memory when searching for all alignments compared to 7 GB in standard mode. MTSv took on average 12 min to run (~35 times faster than Bowtie2 and 5 times faster than Minimap2) and required 9.5 GB of memory (only 8% of what Minimap2 required) to achieve the same or better results (in terms of correct taxonomic assignments). MTSv does an exhaustive search of all candidate alignments but it is taxon-aware so it will stop attempting alignments against reference sequences for a TaxID that already has a valid match. Consequently, MTSv was able to skip over 1 billion (93%) redundant and computationally expensive alignments while still correctly assigning the Staphylococcus reads.

MTSv has high precision and sensitivity and reduces false positive assignments

MTSv had nearly perfect sample purity and much higher average purity compared to Centrifuge and Kraken2 (Fig. 4, rows 1 and 2, respectively). This indicates that most of the assignments made by MTSv at each taxonomic level were correct whereas, in every case Centrifuge and Kraken2 had incorrect assignments all the way up to the kingdom level (assigning reads to a false archaea bin as indicated by the average purity of 0.5). High precision is the advantage of using a full alignment approach and this is particularly clear when looking at the number of false positive bins (Fig. 5, middle row). At the species level, Centrifuge and Kraken2 assigned reads to thousands of additional taxa that were not in the sample whereas the highest number of false bins for MTSv was 605 for the airways sample. Additionally, the sizes of the false bins were generally smaller for MTSv (Fig. 5, bottom row).

While MTSv has a more conservative approach to making assignments (here the alignment of a 150 bp read could only contain up to 20 mismatches compared to Centrifuge which will report a hit if it finds an exact match >22 bp long) it was still able to assign a large portion of the reads at the species level (Fig. 5, first row). As a result, MTSv maintained high sample accuracy and average completeness (Fig. 4, rows 3 and 4, respectively) by keeping the number of false negatives low. MTSv and Kraken2 had similarly high performance for these metrics while Centrifuge was penalized for having more unassigned reads. The oral sample contained viral sequences that were not included in our bacteria-only reference database, so all the tools had lower average completeness for that sample.

MTSv assigns more reads from genomes that diverge from the reference

The limited number of reference genomes available for taxonomic classification means that sequences in metagenomic samples will often share low similarity with existing database sequences. Divergent sequences contain more gaps and mismatches which will reduce the number of exact k-mer matches. As a result, exact matching approaches tend to miss these assignments completely or have a reduced number of reads assigned to these taxa which can influence downstream abundance estimations.

The number of reads assigned accurately by Centrifuge and Kraken2 to Faecalibacterium prausnitzii dropped more rapidly as the genome diverged from the reference sequence (Fig. 6, left). Even at 95% genome similarity, the number of reads assigned dropped to 91% and 94%, respectively, while MTSv was at 98% using default parameters. At 90% similarity, assigned reads dropped to ~50% for Centrifuge and Kraken2 while MTSv still assigned 74%. The reads were spiked into the CAMI2 gastrointestinal sample to show how the missed assignments affected the proportion of reads assigned to F. prausnitzii compared to other species in the sample (Fig. 6, middle). For Centrifuge and Kraken2, the unassigned reads caused F. prausnitzii to drop in rank below three other species in relative sequence abundance when genome similarity was 90% and below five species at 85% (Fig. 6, right). Using MTSv default parameters, F. prausnitzii dropped below only three species at 85% genome similarity.

With default parameters, MTSv outperformed Centrifuge and Kraken2 in assigning reads from genomes that diverge from reference sequences. By using more sensitive search parameters (seed size = 16, seed interval = 2) and relaxing the edit distance cutoff to 20%, MTSv performed even better. With these settings, MTSv assigned 84% of F. prausnitzii reads and read abundance dropped below only one other species when genome similarity was at 85%.

MTSv provides higher confidence in the detection of pathogens

In microbial forensics and clinical applications, it is important to have a high level of confidence that a detected organism is truly present in the sample. The presence of many false positive assignments can erode this confidence and may require additional time consuming and expensive validation tests. Ideally, adding a single organism to a complex background should not result in many additional false positive assignments. However, although each of the tools was able to detect the correct species in our pathogen spike-in experiments, Centrifuge and Kraken2 introduced many additional false positive assignments to near neighbor species.

The total number of reads assigned by each tool was low which is expected when performing taxonomic assignments of a complex soil mixture against a relatively limited collection of reference sequences (Table 2). Both Centrifuge and Kraken2 assigned more reads than MTSv but many of the assignments had very low support. The median number of k-mers mapped per read for Kraken2 was 19 (16.4%), 14 (12.1%), and 16 (13.8%) for B. anthracis, F. tularensis, and Y. pestis spike-ins, respectively and the median hit length for Centrifuge was 24 bp for all three samples. In contrast, all the assignments made by MTSv had 20 or fewer mismatches from the reference sequences. Relaxing the mismatch cutoff to 30 increased the number of assigned reads to 706,182, 680,015, and 604,819, respectively, while still requiring an 80% match with the reference sequence.

Despite having fewer total reads assigned, MTSv made the most assignments to the spike-in taxa at the species level. For B. anthracis and Y. pestis, fewer calls could be made at the species level, likely due to high genetic similarity between these and near neighbor species. For example, Bacillus anthracis shares >99% average nucleotide identity with its closest relatives in the Bacillus cereus group (Helgason et al., 2000; Rasko et al., 2005). When reads were assigned to multiple species in a genus, the assignment was correctly made at the genus level. Unlike the other tools, MTSv was less likely to assign reads to near neighbor species within the same genus and introduced fewer new species assignments.

A comparison to the soil background sample (Fig. 7) shows that the additional near neighbor assignments were likely false positives rather than a true signal in the soil background. Because Kraken2 and Centrifuge had many low confidence assignments, we opted to include filtered results for comparison. We used a confidence filter of 0.1 for Kraken2 (Wood, Lu & Langmead, 2019; Ye et al., 2019) and a score filter of 150 for Centrifuge (Ammer-Herrmenau et al., 2021; Sanderson et al., 2018); neither filter had a large impact on the number of assignments to the spike-in species. The unfiltered Kraken2 results showed large increases (>2-fold) in 13 additional taxa in the Bacillus genus with a maximum of a 16-fold increase in B. wiedmannii. Filtering the low confidence assignments reduced this to eight taxa with more than a 2-fold increase and maximum 9-fold increase for B. thuringiensis. Prior to filtering, Centrifuge showed >2-fold increase in six Bacillus species but after filtering this increased to seven species with a maximum of a 6-fold increase. Centrifuge appeared to perform better prior to filtering but this was because it had many more assignments to Bacillus species in the soil background compared to the other tools and most of them were of very low quality (score < 150, hit length < 32 bp). The higher number of (albeit low quality) assignments for Bacillus species in the background sample masked the increase due to the spike-in prior to filtering them out. MTSv had only three Bacillus species with greater than a 2-fold increase with a maximum of 5-fold.

MTSv is reasonably fast for a full-alignment algorithm and is flexible enough to run on lower memory or limited CPU environments

MTSv offers a lot of flexibility when it comes to managing the resources required for taxonomic assignment against large databases. This begins with the index building step. Unlike Centrifuge, Kraken2, and other tools, MTSv allows the indices to be split up into different sized chunks to accommodate resource limitations. For example, because the Centrifuge index must be built all at once and requires a large amount of RAM that scales with the size of the reference genome, we were restricted to using a smaller reference database in this article to allow direct comparison with Centrifuge. MTSv indices can be built using any reference database or even multiple reference databases because it breaks them into smaller chunks that require less memory to build and run. MTSv can easily combine the assignment results across multiple indices and additional indices can be added without needing to rebuild the entire database.

In this article, we chose to split the reference sequences among 13 different MG-indices for MTSv. Each index required 251 GB of memory and an average of 1 h 23 min to build. Alternatively for more memory limited situations, the reference sequences could be split into 66 ~2 GB files which would each require 49 GB of memory and about 18 min each to build (Table 3). Both approaches required less memory than Centrifuge’s build process which required 317 GB of RAM. The storage size of the MTSv indices (using the default suffix array and occurrence array sampling intervals of 32 and 64, respectively) was ~3.6 times the size of the original reference sequences. Although this is larger than Kraken2 (~0.3×) and Centrifuge (~0.5×), storage space is arguably the least expensive and most readily available resource. The alternative example in Table 3 (MTSv alt) was built with a higher occurrence array sampling rate of 512 so the files took up less space on disk (~2.3× the original FASTA) but they required slightly longer binning times.

Because MTSv performs full alignments, it cannot compete with the speeds of exact matching tools, but it is reasonably fast compared to traditional alignment tools. For taxonomic assignment, the average time per index was about 30 min for MTSv and each required about 34 GB of memory (Table 4). If each of the 13 indices were run simultaneously then the run time would be equal to the time for the longest index (~40 min). If each index was run in series, it would take about 6 and a half hours or about 13 times longer than Centrifuge.