Improved metagenomic analysis with Kraken 2 — TEI-tagged view

Derrick E. WoodpersName, Jennifer LupersName, Ben LangmeadpersName (2019date). “Improved metagenomic analysis with Kraken 2.” Genome Biologytitle·j.

Object of study — Kraken 2term'sterm·var data structures, classification algorithms, resource use, accuracy, and translated-search behavior

Research question — whether minimizerterm·def indexing and a compact probabilistic table can reduce metagenomic-classifier memory and runtime without materially sacrificing accuracy

Contribution — Kraken 2term, a minimizer-based classifier with compact hashing, internal taxonomy compression, spaced seedsterm, translated search, and Brackenterm compatibility

Full text — tagged

Improved metagenomicterm analysis with Kraken 2term Abstract Although Krakenterm’sterm·var kterm·var -mer-based approach provides a fast taxonomic classification of metagenomicterm sequence data, its large memory requirementsterm can be limiting for some applications. Kraken 2term improves upon Kraken 1term by reducing memory usageterm·var by 85%, allowing greater amounts of reference genomic data to be used, while maintaining high accuracy and increasing speed fivefold. Kraken 2term also introduces a translated search modeterm, providing increased sensitivity in viral metagenomicsterm analysis.

Keywords: Metagenomicsterm, Metagenomicsterm classification, Microbiome, Probabilistic data structures, Alignment-free methods, Minimizersterm Assigning taxonomic labelsterm to sequencing readsterm is an important part of many computational genomics pipelinesterm for metagenomicsterm projects. Recent years have seen several approaches to accomplish this task in a time-efficient manner [ 1 – 3 ]. One such tool, Krakenterm [ 4 ], uses a memory-intensive algorithm that associates short genomic substrings ( k -mersterm) with the lowest common ancestorterm (LCAterm) taxa. Krakenterm and related tools like KrakenUniqterm [ 5 ] have proven highly efficient and accurate in independent tool comparisons [ 6 , 7 ]. But Krakenterm’sterm·var high memory requirementsterm force many researchers to either use a reduced-sensitivity MiniKraken databaseterm [ 8 , 9 ] or to build and use many indexes over subsets of the reference sequencesterm [ 10 , 11 ]. Its memory requirementsterm can easily exceed 100 GB [ 7 ], especially when the reference data includes large eukaryotic genomes [ 12 , 13 ]. Here, we introduce Kraken 2term, which provides a major reduction in memory usageterm·var as well as faster classification, a spaced seedterm searching scheme, a translated search modeterm for matching in amino acid space, and continued compatibility with the Brackenterm [ 14 ] species-level sequence abundanceterm estimation algorithm.

Kraken 2term addresses the issue of large memory requirementsterm through two changes to Kraken 1term’sterm·var data structures and algorithms. While Kraken 1term used a sorted list of k -merterm/LCAterm pairs indexed by minimizersterm [ 15 ], Kraken 2term introduces a probabilistic, compact hash tableterm to map minimizersterm to LCAsterm. This table uses one third of the memory of a standard hash table, at the cost of some specificity and accuracy. Additionally, Kraken 2term only stores minimizersterm (of length ℓ , ℓ ≤ kterm·var ) from the reference sequenceterm library in its data structure, whereas Kraken 1term stored all k -mersterm. This change means that, during classification, the minimizerterm·def ( ℓ -mer) is the substring compared against a reference set in Kraken 2term, while Kraken 1term compared k -mersterm (Fig. 1 a, b). Kraken 2term’sterm·var index for a specific reference database with 9.1 Gbp of genomic sequences uses 10.6 GB of memory when classifying. Kraken 1term’sterm·var index for the same reference uses 72.4 GB of memory for classification (Fig. 2 a, Additional file 1 : Table S1). In general, a Kraken 2 databaseterm is about 85% smaller than a Kraken 1term database over the same references (Additional file 2 : Figure S1).

Differences in operation between the two versions of Krakenterm. a Both versions of Krakenterm begin classifying a k -merterm by computing its ℓ bp minimizerterm·def (highlighted in magenta). The default values of kterm·var and ℓ for each version are shown in the figure. b Kraken 2term applies a spaced seedterm mask of sterm·var spaces to the minimizerterm·def and calculates a compact hash code, which is then used as a search query in its compact hash tableterm; the lowest common ancestorterm (LCAterm) taxon associated with the compact hash code is then assigned to the k -merterm (see the “ Methods ” section for full details). In Kraken 1term, the minimizerterm·def is used to accelerate the search for the k -merterm, through the use of an offset index and a limited-range binary search; the association between k -merterm and LCAterm is directly stored in the sorted list. c Kraken 2term also achieves lower memory usageterm·var than Kraken 1term by using fewer bits to store the LCAterm and storing a compact hash code of the minimizerterm·def rather than the full k -merterm. d Impact on speed, memory usageterm·var, and prokaryotic genus F1-measureterm in Kraken 2term when changing kterm·var with respect to ℓ ( ℓ = 31, sterm·var = 7 for all three graphs). e Impact on prokaryotic genus sensitivity and positive predictive valueterm·def (PPV) when changing the number of minimizerterm·def spaces sterm·var ( kterm·var = 35, ℓ = 31 for both graphs). In d and e , the data are from our parameter sweep results in Additional file 1 : Table S2, and the default values of the independent variables for Kraken 2term are marked with a circle.

Comparison between Kraken 2term and other sequence classification tools. a Processing speedterm (in millions of reads per minute) and memory usageterm·var (measured by maximum resident set size, in gigabytes) are shown for each classifier, as evaluated on 50 million paired-end simulated reads with 16 threads. Accuracy results are shown for b 40 prokaryotic genomes and c 10 viral genomes. The results here are shown for sensitivity, positive predictive valueterm·def (PPV), and F1-measureterm as evaluated on a per-fragment basis at the genus rank, with 1000 reads simulated from each genome. The strains from which reads were simulated were excluded from the reference libraries for each classification tool. “Kraken 2Xterm” is Kraken 2term using translated search against a protein database. Full results for these strain-exclusion experiments are available in Additional file 1 : Table S1 Kraken 2term’sterm·var approach is faster than Kraken 1term’sterm·var because only distinct minimizersterm from the query (read) trigger accesses to the hash table. A similar minimizer-based approach has proven useful in accelerating read alignment [ 16 ]. Kraken 2term additionally provides a hash-based subsamplingterm approach that reduces the set of minimizerterm·def/LCAterm pairs included in the table, allowing the user to specify a target hash table sizeterm·var; smaller hash tablesterm yield lower memory usageterm·var and higher classification throughputterm at the expense of lower classification accuracyterm (Fig. 1 d, Additional file 1 : Table S2).

Kraken 2term also features other improvements to accuracy and runtime. A new translated search modeterm (Kraken 2Xterm) uses a reduced amino acid alphabetterm and increases sensitivity on viral datasets compared to nucleotide-based search. Block- and batch-based parsing within the critical section is used to improve thread scaling, in a manner similar to that used in recent versions of Bowtie 2 [ 17 ]. We also added a form of spaced seedterm search and automated masking of low-complexity reference sequencesterm to improve accuracy.

To assess the accuracy and performance of Kraken 2term, we selected 40 prokaryotic and 10 viral genomes for which we had reference genomes for at least 2 sister subspecies and at least 2 sister species (Additional file 1 : Table S3). We then created a reference genome (or protein) set that excluded the 50 taxa for the genomes we selected. This reference set and taxonomy were held constant between the various classifiers we examined, avoiding any confounding due to the differences in the reference database. A similar approach has been recently used for this same purpose in another study [ 7 ].

We simulated 1 million Illumina 100 × 100 nt paired-end reads from each of the 50 selected genomes, for a total of 50 million reads (25 million fragments). We processed these data with 4 nucleotide search-based sequence classification programs (Centrifuge [ 1 ], CLARK [ 2 ], Kraken 1term [ 4 ], and KrakenUniqterm [ 5 ]) and a translated search classifier (Kaijuterm [ 3 ]). We additionally processed these data with Kraken 2term, using several different databases created with different parameters (the “ Methods ” section).

This strain-exclusion approachterm mimics the real-world scenario where reads likely originate from strains that are genetically distinct from those in the database. The addition of simulated sequencing errorsterm also provides further genetic distance between the test data and the reference sequencesterm. Through this approach, we sought to avoid overly optimistic estimates of a classifier’sterm·var performance.

We found that Kraken 2term exhibited similar, and often superior, per-sequence accuracyterm to the other nucleotide classifiers and that Kraken 2Xterm provided similar (though slightly lower) accuracy compared to Kaijuterm (Fig. 2 b, Additional file 1 : Table S1). The nucleotide-based classifiers exhibited lower accuracy on the viral read dataterm than did the translated search classifiers, demonstrating the advantage of translated search in scenarios marked by high genetic variability and sparsity of available reference genomes [ 3 ].

In some cases, we found that Kraken 2term would not classify a large proportion of reads correctly at the species level, despite the presence of at least two sister strains in the reference database (Additional file 2 : Figure. S2). This was often the result of classifications that were either incorrect at the species level or correct but only made at the genus level (or higher). Such classifications can occur when genomes from different species or genera share a high genomic identity, which is the case in multiple places of the taxonomy, including the Shigella [ 18 ], Bacillus [ 19 ], and Pseudomonas [ 20 ] genera. A redefinition of the taxonomy based on the phylogenyterm as recently proposed [ 21 ] would likely improve sensitivity at the species levelterm.

Following our evaluation of the classifiers’ accuracy, we then examined the runtime and memory requirementsterm of each program. Kraken 2term provided substantial increases in processing speedterm, classifying paired-end data at over 93 million reads per minute while using 16 threads, a speed over 5 times faster than Kraken 1term, the next-fastest classifier (Fig. 2 a, Additional file 1 : Table S1). Additionally, Kraken 2term exhibited superior thread scaling to Kraken 1term (Additional file 1 : Table S4). Kraken 2term’sterm·var memory requirementterm is also 15% of Kraken 1term’sterm·var, and only 2.5 times as much as that of the least memory-intensive classifier we examined, Centrifuge. With respect to the translated search programs, Kraken 2Xterm is over 3 times faster and uses 47% less memory than Kaijuterm.

To determine if Kraken 2term exhibited similar analytical performance on real sequencing data, we classified read data from the FDA-ARGOS project [ 22 ]. We compared the fragment classifications obtained by the various classification programs to the taxonomic labelsterm attached to the corresponding ARGOS experiment. Kraken 2term exhibits similar genus-level concordance and discordance statistics to the other nucleotide search classifiers, while Kraken 2Xterm exhibits similar but less agreement with the ARGOS labels than does Kaijuterm (Additional file 1 : Table S5). These results agree with those obtained in the strain-exclusion experiment on simulated data.

As a continuation of the strain-exclusion experiments, we applied Brackenterm [ 14 ] to the Kraken 1term and Kraken 2term results, estimating species- and genus-level sequence abundanceterm for prokaryotic species. Brackenterm uses a Bayesian algorithmterm to integrate reads Krakenterm classified at higher taxonomic levels into the abundance estimates. Although the true strain-level taxa are excluded from the database, Brackenterm recaptured most of the true genus-level and species-level sequence abundancesterm using both Kraken 2term and Kraken 1term classification results. Comparing the results, the Brackenterm estimates were more accurate with Kraken 2term than with Kraken 1term at both the genus and species levels, likely owing to Kraken 2term’sterm·var higher sensitivityterm (Additional file 2 : Figure S2). Brackenterm ran in less than 1 sterm·var, a minute fraction of the runtime of any of the classification programs we examined.

As databases of assembled genomes continue to grow, databases of reference sequencesterm used for metagenomicsterm studies will also grow [ 21 , 23 ]. We presented Kraken 2term, an extremely memory-efficient metagenomicsterm classification tool that replaces Kraken 1term’sterm·var k -merterm database with a probabilistic data structure that is substantially smaller, allowing six to seven times more reference data compared to Kraken 1term. The algorithms introduced in Kraken 2term to subsample the set of genomic substrings also provide Kraken 2term with the ability to further reduce the size of its database and accelerate the processing of sequencing data. We showed Kraken 2term’sterm·var accuracy is comparable to that of Kraken 1term and other competing tools, consistent with other studies [ 6 , 7 ]. We also showed that its new translated search modeterm has accuracy approaching that of the protein-focused Kaijuterm tool, while using less memory and runtime. Also, Kraken 2term is compatible with the Brackenterm software for species-level quantification, making Kraken 2term straightforwardly usable for that application.

In the future, it will be important to consider additional use cases for Kraken 2term. For example, other data structures similar to our compact hash tableterm, such as the counting quotient filter [ 24 ], could be implemented and used in computing environments and applications that may benefit from a particular data structure’sterm·var design and properties. Additionally, the KrakenUniqterm [ 5 ] tool uses the HyperLogLog sketch [ 25 ] to estimate the number of distinct k -mersterm matched at each node of the taxonomy, a statistic that is used in turn to better determine the presence or absence of individual genomes. We plan to add this functionality in the future, as it enables applications in the diagnosis of infections where the infectious agent is present at low abundance.

Methods Compact hash tableterm The hash table used by Kraken 2term to store minimizerterm·def/LCAterm key-value pairsterm is very similar to a traditional hash table that would use linear probingterm for collision resolution, with some modifications. Kraken 2term’sterm·var compact hash table (CHT)term uses a fixed-size array of 32-bit hash cellsterm to store key-value pairsterm. Within a cell, the number of bits used to store the value of the key-value pairterm will vary depending on the number of bits needed to represent all unique taxonomy ID numbers found in the reference sequenceterm library; this was 17 bits with the standard Kraken 2 databaseterm in September 2018. The value is stored in the least significant bits of the hash cellterm and must be a positive integer. Values of 0 represent empty cells. Within the remaining bits of the hash cellterm, the most significant bits of the key’sterm·var hash code (a compact hash code) are stored. Searching for a key Kterm·var in the CHTterm is done by computing the hash code of the key h ( Kterm·var ) then linearly scanning the table array starting at position h ( Kterm·var ) mod | T | (where | T | is the number of cells in the array) for a matching key. Examples of this search process—including both key/value insertion and querying—are shown in Additional file 2 : Figure S3. In Kraken 2term, the hash function h used is the finalization function from MurmurHash3 [ 26 ].

Compacting hash codes in this way allows Kraken 2term to use 32 bits for a key-value pairterm, a reduction compared to the 96 bits used by Kraken 1term (64 bits for key, 32 for value) (Fig. 1 c). But it also creates a new way in which keys can “collide,” in turn impacting the accuracy of queries to the CHTterm. Two distinct keys can be treated as identical by the CHTterm if they share the same compact hash code and their starting search positions are close enough to cause a linear probe to encounter a stored matching compact hash code before an empty cell is found. This property gives the CHTterm its probabilistic nature, in that two types of false-positive query resultsterm are possible: either (a) a key that was not inserted can be reported as present in the table or (b) the values of two keys can be confused with each other. In Kraken 2term, the former error is indeed a false positiveterm, whereas the latter results in a less specific LCAterm being assigned to the minimizerterm·def (Additional file 2 : Figure S3). The probability of either of these errors is < 1% with Kraken 2term’sterm·var default load factorterm·var of 70% (Additional file 2 : Figure S4). The adverse effect on read-level classification is further mitigated by the algorithm Kraken 2term uses to combine information from across the read, which is unchanged from Kraken 1term and utilizes information from all k -mersterm in a sequence to counteract low-frequency erroneous LCAterm values that could be returned by a key-value store.

The probabilistic nature and comparisons involving parts of a key’sterm·var hash code make the CHTterm similar to the counting quotient filter (CQF) described by Pandey et al. [ 24 ] Like the CQF, Kraken 2term’sterm·var CHTterm features high locality of memory access during an individual query due to the linear probingterm that the CHTterm employs. Unlike the CQF, however, our CHTterm does not allow the full hash code to be recovered from a stored value (the CQF’sterm·var remainder ), and so we are unable to resize a CHTterm once it is instantiated. Additionally, our CHTterm has an additional possibility of error compared to the CQF, where two keys that do not have the same full hash code but share a truncated hash code will be treated as identical. The CQF can avoid such “soft” hash collisions.

Internal taxonomy of a Kraken 2 databaseterm While Kraken 1term used the taxonomy provided by the user without modification, Kraken 2term makes some modifications to its internal representation of the taxonomy that causes that representation to differ from the user-provided taxonomy. First, Kraken 2term finds a minimal set of nodesterm in the user-provided taxonomy. This minimal set consists of all nodes to which a reference sequenceterm is assigned, as well as all of those nodes’ ancestors; vertices between nodes in this set remain as they were in the user-provided taxonomy, maintaining the tree structure in the internal representation. Kraken 2term then assigns nodes in the minimal set sequentially increasing internal taxonomy ID numbers using a breadth-first searchterm (BFS) beginning at the root, with the root having an internal ID number of 1. This BFS provides a guarantee that ancestor nodes will have smaller internal ID numbers than their descendants; an example of this numbering is shown in Additional file 2 : Figure S3. Kraken 2term stores a mapping of its internal taxonomy numbers to the external taxonomy ID numbersterm to make its results more easily interpretable, and performs all output using the external taxonomy ID numbersterm.

Kraken 2term’sterm·var use of this internal taxonomy representationterm allows for the easier computation of the LCAterm of two nodes because the ID numbers themselves give information as to their relative depths in the tree, while the National Center for Biotechnology Information (NCBI) taxonomy IDs lack this property. The internal taxonomy representationterm also allows Kraken 2term to use the minimal number of bits for storage of taxonomy ID numbers, giving maximal space for the compact hash codes and reducing the probability of CHTterm errors (or “hash table collisionsterm,” as we describe elsewhere in this paper).

A Kraken 2 databaseterm consists of a CHTterm and this internal taxonomy representationterm. Typical databases will be built using the NCBI taxonomy [ 27 ], but users can override this default to create custom databases for atypical use cases.

Minimizer-based subsamplingterm In contrast to Kraken 1term’sterm·var use of all k -mersterm in the standard use case, Kraken 2term subsamples the set of genomic substrings and inserts only the distinct minimizersterm into its database (Fig. 1 b). We define the ℓ bp minimizerterm·def of a k -merterm ( ℓ ≤ kterm·var ) to be the lexicographically smallest canonical ℓ -mer found within the k -merterm. An ℓ -mer is called canonical if it is lexicographically less than or equal to its reverse complement. Note that if kterm·var = ℓ , no subsampling occurs and Kraken 2term inserts the same substrings into its data structure that Kraken 1term would. Additionally, as the difference between kterm·var and ℓ grows, fewer substrings are inserted into the CHTterm, reducing its size along with Kraken 2term’sterm·var memory usageterm·var and runtime (Fig. 1 d, Additional file 1 : Table S2). The default values for Kraken 2term, kterm·var = 35 and ℓ = 31, were determined after the analysis of the parameter sweep results we show in Additional file 1 : Table S2.

Kraken 2term determines which ℓ -mers are minimizersterm by the use of a sliding window minimum algorithmterm, in contrast to Kraken 1term’sterm·var implementation which examined each k- mer anew. This allows for a faster determination of minimizersterm, as less work is required when moving from one k- mer to the next overlapping k- mer (in terms of computational complexity, the new approach uses an average of O (1) time to calculate a new minimizerterm·def vs. Θ ( kterm·var ) time with the older algorithm). The sliding window minimum calculation uses a double-ended queue (or “deque”) in which canonicalized candidate ℓ -mers are inserted in the back, along with the candidates’ position in the original sequence. As a new candidate is encountered, enqueued candidates are removed from the back of the deque until the candidate at the back has a greater value than the new candidate (as determined by lexicographical ordering). The new candidate is then pushed onto the back of the deque.

Once a k- mer’sterm·var worth of ℓ -mers has been processed in this way, the front of the deque contains the minimizerterm·def of that k- mer. This property is then maintained during scanning subsequent bases by removing the front element in the deque if it is from a position in the original sequence that is not in the current k- mer. In this way, the front element of the deque holds the minimizerterm·def of the k- mer currently being examined.

We further augmented the sliding window algorithm to include the exclusive or (XOR) shuffling operation from Kraken 1term. This operation serves to permute the ordering of the ℓ -mers when calculating minimizersterm and helps to avoid a bias toward low-complexity ℓ -mers when selecting the minimizerterm·def of a k -merterm [ 4 , 15 ]. To shuffle, we calculate the XOR value of the ℓ -mer and a predefined constant and use this value as the “candidate” that is put in the deque. When the original ℓ -mer value is needed again, the operation is reversed by XORing a second time with the same constant.

Spaced seedterm usage Spaced k -mersterm, a similar concept to spaced seedsterm, have been shown to improve the ability to classify reads within the Krakenterm framework [ 28 ]. Kraken 2term uses a simple spaced seedterm approach where a user specifies an integer sterm·var when building a database that indicates how many positions in the minimizerterm·def will be masked (i.e., not considered when searching). Beginning with the next-to-rightmost position in the minimizerterm·def, every other position is masked until sterm·var positions have been masked. For example, if sterm·var = 3 and ℓ = 12, the positions in the bit string 1111 1101 0101 with a “0” would be masked. When using Kraken 2term, Kraken 1term’sterm·var classification results can be most closely approximated by setting kterm·var = ℓ = 31 and sterm·var = 0, as these settings will avoid any minimizer-based subsamplingterm and spaced seedterm usage. Kraken 2term’sterm·var default value for sterm·var is 7 and was determined after the analysis of the parameter sweep results we show in Additional file 1 : Table S2.

The canonical ℓ -mers that are minimizerterm·def candidates are masked with the spaced seedterm mask prior to their insertion into the deque for the sliding window calculation. By performing canonicalization of the minimizerterm·def candidates prior to applying the spaced seedterm mask, we ensure the result is the same whether applied to the ℓ -mer or its reverse complement.

Kraken 1term’sterm·var sensitivity performance was governed by the value of kterm·var (the length of the searched substring). By comparison, the use of spaced seedsterm and minimizer-based subsamplingterm means that Kraken 2term’sterm·var sensitivity performance will be largely governed by ℓ - sterm·var (the number of compared bases in Kraken 2term’sterm·var searched substring). Thus, increasing sterm·var will generally increase sensitivity while decreasing positive predictive valueterm·def (Fig. 1 e, Additional file 1 : Table S2).

Hash-based subsamplingterm Kraken 2term estimates the required capacity of the hash table given the kterm·var , ℓ , and sterm·var values chosen along with the sequence data in a database’sterm·var reference genomic library. Some users will not have access to large memory computers, and therefore, this estimate may be greater than the maximum possible hash table size that they can work with. To aid such users, Kraken 2term allows them to specify a maximum size when building a database. If the estimated required capacity is larger than the maximum requested size, then the minimizersterm will be subsampled further using a hash function. Given an estimated required capacity Sterm·var ′ and a maximum user-specified capacity of Sterm·var ( Sterm·var < Sterm·var ′), we can calculate the value f = Sterm·var / Sterm·var ′, which is the fraction of available minimizersterm that the user will be able to hold in their database. A minimum allowable hash value of v = (1 − f )∙ M can also be calculated, where M is the maximum value output by hash function h . Any minimizerterm·def in the reference library with a hash code less than v will not be inserted into the hash table. This value v is also provided to the classifier so that only minimizersterm with hash codes greater than or equal to v will be used to probe the hash table, saving the search failure runtime penalty that would be incurred by searching for minimizersterm guaranteed not to be in the hash table.

Evaluation of k -merterm level discordance rates At a k -merterm level, there are two main types of discordance between Kraken 1term and Kraken 2term’sterm·var results: those caused by two distinct k -mersterm sharing the same minimizerterm·def (a “minimizer collisionterm”) and those caused by two distinct minimizersterm being indistinguishable by the CHTterm (a “hash table collisionterm”). Minimizer collisionsterm are not always damaging. When it occurs between k -mersterm from very closely related genomes, such a collision might detect true homology even in the face of single nucleotide polymorphisms and/or sequencing error. That said, minimizer collisionsterm between k -mersterm from distantly related genomes could produce either elevated LCAterm values (if both genomes are in the reference library) or incorrectly classified k -mersterm (if one of the genomes is not in the reference library). Hash table collisionsterm are a consequence of the probabilistic nature of the CHTterm and can also cause either elevated LCAterm values or incorrectly classified k -mersterm (Additional file 2 : Figure S3). We note that these different discordant results are all at a k -merterm level and may not always affect a query sequence’sterm·var classification due to the many k -mersterm’ worth of data that are used to classify a query sequence; aside from slight modifications to handle the subsampling methods we use in Kraken 2term, the classification method of Kraken 2term is identical to Kraken 1term.

We wished to estimate the rate at which these collisions would cause discordance at a k -merterm level between the Kraken 1term and Kraken 2term results. To do so, we selected a specific bacterial genome for which we had neighboring genomes at each taxonomic rankterm·var from species to phylum. The selected genome was our “reference sequenceterm,” and eight others were progressively more taxonomically distant from the reference sequenceterm. We list the nine genomes used in these experiments in Additional file 1 : Table S6. We additionally created a synthetic genome with 4 Mbp of uniformly random DNA. Together, these ten sequences formed a set of “query sequences” and were the basis for our evaluation of collision rates. For these experiments, we used the default Kraken 2term values of kterm·var = 35, ℓ = 31, and sterm·var = 7, unless otherwise noted.

To determine the rates of discordance caused by minimizer collisionsterm, we compared each of the ten query sequences’ k -mersterm to the set of reference sequenceterm k -mersterm. For each sequence, the minimizer collisionterm rate is the proportion of distinct k -mersterm in a query sequence that (a) are not in the set of reference sequenceterm k -mersterm and (b) share a minimizerterm·def with a reference sequenceterm k -merterm. The various sequences’ minimizer collisionterm rates are summarized in Additional file 1 : Table S7. We hypothesized that the minimizer collisionterm rate would be influenced by the length of the minimizerterm·def used, due to the length’sterm·var direct relationship to the number of possible minimizersterm. To test this, we repeated the minimizer collisionterm rate estimation experiment focusing on the reference genome and using the random synthetic genome as the sole query sequence. Setting kterm·var = 35 and sterm·var = 0, we varied the ℓ parameter from 8 to 31. Minimizerterm·def lengths greater than 15 had collision rates under 1%. Minimizerterm·def lengths greater than 22 had 0 collisions. The full results are shown in Additional file 2 : Figure S5.

To determine the rates of discordance caused by hash table collisionsterm, we compared each of the ten query sequences’ minimizersterm to a CHTterm populated with the reference sequenceterm minimizersterm. The CHTterm was created with a load factorterm·var of 70% and 15 bits reserved for the truncated hash code (the same parameters used in Kraken 2term’sterm·var standard database in September 2018). For each sequence, the hash table collisionterm rate is the proportion of distinct minimizersterm in a query sequence that (a) are not minimizersterm in the set of reference sequenceterm minimizersterm and (b) are reported by the CHTterm as being inserted in the hash table. The various sequences’ hash table collisionterm rates are summarized in Additional file 1 : Table S8. To investigate the impact of load factorterm·var and truncated hash code size on hash table collisionterm rates, we repeated the hash table collisionterm rate experiment, but focused only on the reference genome and used the random synthetic genome as the sole query sequence. We used the same default values of kterm·var , ℓ , and sterm·var as before (35, 31, and 7, respectively) and calculated hash table collisionterm rates while varying both the load factorterm·var and truncated hash code size. The impact of these two parameters on hash table collisionterm rates is shown in Additional file 2 : Figure S4. The parameters adopted for Kraken 2term’sterm·var default mode had an error rate of 0.016%, consistent with the results seen when comparing genomes of different species (Additional file 1 : Table S8).

Processing of a standard genomic reference library The CHTterm’sterm·var modest memory requirementsterm, and the additional savings yielded by minimizer-based subsamplingterm, allow more reference genomic data to be included in Kraken 2term’sterm·var standard reference library. Whereas Kraken 1term’sterm·var default database had data from archeal, bacterial, and viral genomes, Kraken 2term’sterm·var default database additionally includes the GRCh38 assembly of the human genome [ 29 ] and the “UniVec_Core” subset of the UniVec database [ 30 ]. We include these in Kraken 2term’sterm·var default database to allow for easier classification of human microbiome reads and more accurate classification of reads containing vector sequences.

Additionally, we have implemented masking of low-complexity sequencesterm from reference sequencesterm in Kraken 2term, by using the “dustmasker” [ 31 ] (for nucleotide sequences) and “segmasker” [ 32 ] (for protein sequences) tools from NCBI. Using the tools’ default settings, nucleotide and protein sequences are checked for low-complexity regions, and those regions identified are masked and not processed further by the Kraken 2 databaseterm building process. In this manner, we seek to reduce false positivesterm resulting from these low-complexity sequencesterm, similar to the build process for Centrifuge [ 1 ].

Populating the Kraken 2term hash table Kraken 2term begins building a CHTterm by first estimating the number of distinct minimizersterm present in the reference library for the selected values of kterm·var , ℓ , and sterm·var . This is done through a form of zeroth frequency moment estimation [ 33 ] where Kraken 2term creates a small set structure implemented with a traditional hash table. In this set Q , we insert only the distinct minimizersterm that satisfy the criterion h ( m ) mod F < E , where h ( m ) is the hash code of the minimizerterm·def m and E ≪ F (in practice, Kraken 2term uses E = 4 and F = 1024). We then find the estimate of the total number of distinct minimizersterm by multiplying the number of satisfactory distinct minimizersterm (| Q |) by F / E . This form of estimation requires storing in memory only a fraction of all distinct minimizersterm (approximately E / F ) and allows us to quickly set the capacity of our CHTterm properly without needing to first store all elements in it.

After estimating the number of distinct minimizersterm D = | Q |( F / E ) present in the reference library, Kraken 2term then allocates memory for a CHTterm containing D /0.7 hash table cells. We selected the divisor of 0.7 so that the resultant hash table will have approximately 30% of its cells remain empty after the population of the CHTterm (i.e., the CHTterm will have a load factorterm·var of 70%). As stated earlier, the cells of this table are 32 bits each, and so the total memory required for Kraken 2term’sterm·var CHTterm is 32 D /0.7 bits or 4 D /0.7 bytes.

Kraken 2term then proceeds to scan each genome in the reference library. Each genome must be associated with a taxonomic ID number so that Kraken 2term can calculate LCAterm values; genomes without associated taxonomy IDs are therefore not processed by Kraken 2term. For a minimizerterm·def M in a genome G, Kraken 2term attempts to insert a key-value pairterm containing M (key) and the taxonomic ID T (value) associated with G into the CHTterm. If the CHTterm does not report that M was previously inserted, then the < M , T > key-value pairterm will be inserted, indicating that the LCAterm of M is currently T . If M was previously inserted into the CHTterm, with LCAterm value T* , then its associated LCAterm value is updated to equal the LCAterm of T and T* . All minimizersterm are processed in this way; once the reference library’sterm·var minimizersterm are all processed, the LCAterm values are properly set for each of the minimizersterm and the database build is complete. The LCAterm operation is both commutative and associative, facilitating parallel index construction.

Classification of a sequence fragment with Kraken 2term Kraken 2term classifies sequence fragments similarly to Kraken 1term, with modifications to facilitate minimizer- and hash-based subsamplingterm. For each k -merterm in an input sequence, Kraken 2term finds its minimizerterm·def and, if it is distinct from the previous k -merterm’sterm·var minimizerterm·def, uses it as a key to probe the CHTterm. If the minimizerterm·def matches a key in the CHTterm, Kraken 2term considers the associated LCAterm value to be the k -merterm’sterm·var LCAterm (Fig. 1 b). Classification then proceeds in the same manner as Kraken 1term, taking note of how many k -merterm hits mapped to each taxon, constructing a pruned classification treeterm, and using the leaf of the maximally scoring root-to-leaf path of that tree to classify the sequence [ 4 ]. If hash-based subsamplingterm was used to build the CHTterm, each minimizerterm·def has its hash code compared against the table’sterm·var maximum allowable hash code, and minimizersterm with higher-than-allowed hash codes are not searched against the CHTterm. Any k -merterm containing an ambiguous nucleotide code is also not searched against the CHTterm.

We note that although Kraken 2term only uses the minimizerterm·def to query the CHTterm, the LCAterm found via this query is assigned by Kraken 2term to the k -merterm rather than only the minimizerterm·def. This means that a stretch of n overlapping k -mersterm that share a minimizerterm·def will all be assigned the same LCAterm value by Kraken 2term and that n hits to that LCAterm will be part of the classification treeterm, even though only one distinct minimizerterm·def was present among the k -mersterm.

Parsing of input files Previous work by LangmeadpersName et al. [ 17 ] has shown the importance of removing parsing work from critical sections, i.e., portions of the program that can be executed by only 1 thread at a time. Kraken 2term uses 2 different methods to defer a majority of parsing work from the critical section to thread-local execution. The first method (referred to as “batch deferred” parsing by LangmeadpersName et al.) reads a set number of lines (40,000 in Kraken 2term) of input in a thread-local buffer within the critical section and then parses the input within a single thread’sterm·var execution. This method is used to perform reading of paired-end FASTQ input, where the lengths of a fragment’sterm·var mates can be different and reading a consistent number of lines from both input files is necessary to ensure a thread is working with complete mate pairs. For FASTA or single-end FASTQ input, Kraken 2term instead uses a more efficient method that reads in a set number of bytes (3 MB in Kraken 2term) of input into a thread-local buffer within the critical section and continues reading input into that buffer until a record boundary is found, at which point a thread leaves the critical section and parses its input. These modifications allow Kraken 2term to more efficiently use multiple threads than did Kraken 1term (Additional file 1 : Table S4).

Translated search To perform a translated search, Kraken 2Xterm first builds a database from a set of reference proteins in the same manner that Kraken 2term does for nucleotide sequences. The usual alphabet of 20 amino acids is reduced to 15 using the 15-character alphabet of Solis [ 34 ]; we add a single additional value representing selenocysteine, pyrrolysine, and translation termination (stop codons). This gives us 16 characters in our reduced alphabet, allowing us to represent a character with 4 bits. Minimizersterm of reference proteins are calculated using the same methods for nucleotide sequences (i.e., using spaced seedsterm if requested and a sliding window minimum algorithmterm), but reverse complements are not calculated and by default kterm·var = 15, ℓ = 12, and sterm·var = 0.

When searching against a protein minimizerterm·def database, Kraken 2Xterm translates all six reading framesterm of the input query DNA sequences into the reduced amino acid alphabetterm. Minimizersterm from all six frames are pooled and used to query the CHTterm, and therefore, all contribute to the Kraken 2Xterm classification of a query sequence.

Generation of data for strain exclusion experiments We downloaded the reference genome and protein data used for the clade exclusion experiments from NCBI in January 2018 from the archaeal, bacterial, and viral domains. We also downloaded the taxonomy from NCBI at this same time. Using the taxonomy ID information for each sequence, we obtained a set of all taxonomy IDs represented by the reference genomes. From this set, we selected a subset of “eligible strains” that had both two sister sub-species taxa present and two sister species taxa present in the set of reference genomes. We selected this subset by examining only those nucleotide sequences with the phrase “complete genome” in their FASTA record header but excluding those that were plasmids or second or third chromosomes. In this manner, we sought to ensure we did not count a genome multiple times due to multiple sequences being associated with that genome. From the eligible strain subset, 40 prokaryotic taxonomy IDs and 10 viral taxonomy IDs were selected arbitrarily to be the strains of origin for our experiments. The strains selected are listed in Additional file 1 : Table S3.

After selecting the taxonomy IDs that represented the strains of origin, we gathered all of the nucleotide sequences we had downloaded—including chromosome and plasmid sequences excluded from our examination when creating the eligible strain subset—into a single file and did the same for the protein sequences. For both the nucleotide and protein files, we placed sequences with taxonomy IDs that were outside the strains of origin into a strain exclusion reference file. Then, for each taxonomy ID in our strain of origin set, we created a single “strain reference” file containing all nucleotide sequences that were associated with that taxonomy ID.

We used Mason 2 [ 35 ] to simulate 100-bp paired-end Illumina sequence data from our strains of origin, with 500,000 fragments being simulated from each strain. When simulating the reads, we used the default options for simulating sequencing errors with Mason 2’sterm·var “mason_simulator” command. These defaults caused the simulator to simulate sequencing errors at rates of 0.4% for mismatches, 0.005% for insertions, and 0.005% for deletions. We combined simulated reads from the strains of origin into a single set of read data. We also shuffled the order of the fragments in this set to control for ordering effects that might affect runtime.

Execution of strain exclusion experiments To evaluate the accuracy and computational performance of Kraken 2term, we compared it to Kraken 1term and several other programs. In selecting these programs, we concentrated on three main properties. First, because Krakenterm’sterm·var principal aim is to provide high-speed taxonomic sequence classification, we looked for taxonomic sequence classification tools that were high in classification speed (within approximately an order of magnitude of Kraken 1term). Secondly, because our experiments rely on holding fixed the reference data between programs, we selected tools which had the ability to customize the underlying reference sequenceterm set and taxonomy using whole-genome reference data. These two requirements led to our selection of KrakenUniqterm, CLARK, Centrifuge, and Kaijuterm as comparator programs. We note that these requirements exclude an accuracy evaluation against programs that are not taxonomic sequence classifiers (programs that output a mapping of sequences to taxa). Sequence abundanceterm estimation programs (which map taxa to sequence counts or frequencies), such as Brackenterm, and population abundance estimation programs (which map taxa to organism counts or frequencies), such as MetaPhlAn [ 36 ], are answering related but different problems than those in our comparator set. For example, Brackenterm does not actually change any of the taxonomic labelsterm associated with the sequenced fragments but rather adjusts the fragment counts associated with low-rank taxa. We also note that although MetaPhlAn does, as part of its operation, classify a small proportion of reads that map to marker genes, this proportion can be less than 10% of reads [ 6 ] in whole-genome shotgun metagenomicterm experiments (such as ours), and thus, MetaPhlAn would yield far lower per-sequence sensitivity relative to the tools in our comparison.

In brief, we used the nucleotide search-based classification programs (Kraken 1term, KrakenUniqterm, Kraken 2term, CLARK, and Centrifuge) to build a strain-exclusion database from reference genomes, and we used the translated search-based classification programs (Kraken 2Xterm and Kaijuterm) to build a strain-exclusion database from reference protein sequences. We compared Kraken 2term and Kraken 2Xterm (both using the code base from Kraken 2term.0.8) against Kraken 1term.1.1, KrakenUniqterm 0.5.6, CLARK 1.2.4, Centrifuge 1.0.3-beta, and Kaijuterm 1.5.0. Because CLARK requires a rank to be specified at the time of building a database, and our evaluations center on genus-rank accuracy, we built a CLARK database for the genus rank for our evaluation work in this paper.

Classifiers received the simulated read data as paired-end FASTQ input. To evaluate runtime and memory usageterm·var, we sought to eliminate the performance impact of reading or writing from disk or from a network storage location. To accomplish this, we copied simulated read data and classifier databases onto a random access memory (RAM) filesystem and directed the classifiers to read input from and write output to that RAM filesystem.

Accuracy was evaluated on a smaller subset of the simulated data containing 1000 fragments per genome of origin or 50,000 fragments in total. To obtain processing speedterm and memory usageterm·var information, we ran each classifier using 16 threads on 25 million sequences’ worth of simulated read data. We used the taskset command to restrict each classifier to the appropriate number of processors (e.g., “taskset -c 0-15” was used with our 16 thread experiments); this ensures that a classifier that uses an external process to aid in its execution has that process’ runtime properly counted against its runtime here. The “/usr/bin/time -v” command provided us with elapsed wall clock time and maximum resident set size data (memory usageterm·var) for each experiment and allowed us to verify that no major page faults were incurred by a classifier during its execution (the absence of which indicates minimal disk- or network-related input/output effects on the runtime). Classifiers were run on a computer with 32 Xeon 2.3 GHz CPUs (16 hyperthreaded cores) and 244 GB of RAM.

Evaluation of accuracy in strain exclusion experiments We evaluated the accuracy of each classifier at a per-fragment level, with respect to a particular taxonomic rankterm·var. Each fragment had a known true subspecies taxon of origin, which implied a true taxon of origin at both the species and genus ranks, which is where we measured accuracy. We now describe how we counted true-positive (TP), false-negative (FN), vague positive (VP), and false-positive (FP) results at the genus and species levels. We describe this at the genus level specifically, but the analogous procedure was also used at the species level. For a given true genus of origin, a TP classification is a classification at that genus or at a descendant of that genus. Because we excluded the strains of origin from our reference databases, we expected all classifiers to make incorrect strain-level classifications and so allow classifications of descendants of the true genus to be judged as TP. We define an FN classification as a failure of a classifier to assign any classification to a sequence and a VP classification as a classification at an ancestor of the true genus of origin. Finally, we define an FP classification as a classification that is incorrect, that is, not at the true genus of origin nor an ancestor or descendant of that true genus. These four categories are mutually exclusive, and all fragments run through a classifier will have their classification (or lack thereof) categorized by one of these categories.

These categories are different from those typically used for binary classification problems; they are used here because these methods can make classifications that are not at leaves of the taxonomic tree but are still correct. For example, a classification of an Escherichia coli fragment as Escherichia would be evaluated as TP for genus-rank accuracy, but as VP for species-rank accuracy. Classification of that same fragment as Vibrio would be evaluated as FP at any rank below class (because the LCAterm of Vibrio and Escherichia is the class taxon Gammaproteobacteria) and would be evaluated as TP for the class rank and above.

Using these categories, we define rank-level sensitivityterm as the proportion of input fragments that were true-positive classificationsterm·def, or TP/(TP + VP + FN + FP)term. We define rank-level positive predictive valueterm (PPV) as the proportion of classifications that were true positives (excluding vague positives), or TP/(TP + FP)term. Along with these definitions of rank-level sensitivityterm and PPV, we also define an F1-measureterm as the harmonic meanterm of those two values.

Evaluation of thread scaling efficiency To evaluate Kraken 1term’sterm·var and Kraken 2term’sterm·var ability to efficiently use multiple threads, we performed an experiment using the strain exclusion databases and simulated read data we describe previously in this section. We ran both Kraken 1term and Kraken 2term on the same data using 1, 4, and 16 threads. The 2 programs were run once on the data as paired-end read data and once as single-end read data. Read data and Krakenterm database files were all placed on a RAM filesystem, and the “taskset” command was used to limit the classifier programs to only as many cores as the number of threads being used. These conditions mirror those of our main strain exclusion experiments, only varying the number of threads between the various runs of the classifiers. The results for this experiment are shown in Additional file 1 : Table S4. In short, Kraken 2term exhibits superior speedup with respect to the number of threads allocated compared to Kraken 1term. This is especially true for paired-end reads.

FDA-ARGOS experimental concordance evaluation The FDA-ARGOS (dAtabase for Reference Grade micrObial Sequences) project provides sequencing experiments for many microbial isolates [ 22 ]. We used the NCBI’sterm·var Sequence Read Archive [ 37 ] to find all 1392 experiments related to the FDA-ARGOS project (accession PRJNA231221). Because some tools are unable to properly process reads of differing lengths, we selected only those 263 experiments that were run on an Illumina HiSeq 4000 instrument and produced 151-bp reads. We then randomly selected 1 experiment from each genus to download and used reservoir sampling to select a subset of 10,000 paired-end fragments from each selected experiment. We also removed experiments for which our strain-exclusion reference genome set did not have a reference genome of the same species as the sequenced isolate. These steps yielded 25 experiments’ worth of data, for 250,000 paired-end fragments in total. Using the strain-exclusion databases created earlier, we then used each classifier to classify the data and examined the percentage of each experiment’sterm·var fragments that were classified.

Because the FDA-ARGOS data are from real sequencing experiments, several factors could explain discordance between a classifier’sterm·var results and the experiments’ assigned taxa, including the evolutionary distance between sequences and reference data, low-quality sequencing runs, and contamination. The true causes of such discordance may not be discernable, and even when they are, they often require an in-depth examination of the sequencing and reference data. For these reasons, we do not report sensitivity and PPV for these data because we cannot be certain of the true taxonomic origin of each individual fragment of real sequencing data. Rather, we evaluated the concordance of the SRA-assigned taxa with the fragments’ classifications at the genus rank and report for each classifier the following quantities: (a) the percentage of fragments with a concordant classification at the genus rank, (b) the percentage of fragments with a discordant classification at the genus rank, (c) the percentage of fragments with a classification of an ancestor of the SRA-assigned genus taxon, and (d) the percentage of fragments that were not classified. The results of this concordance evaluation are provided in full in Additional file 1 : Table S5.

Parameter sweeps We examined various values for parameters to ensure Kraken 2term’sterm·var default parameters would provide an advantageous balance of accuracy, classification speed, and memory usageterm·var. Specifically, we looked at parameters relating to minimizer-based subsamplingterm ( kterm·var and ℓ ), hash-based subsamplingterm ( f = Sterm·var / Sterm·var′ ), and spaced seedterm usage ( sterm·var ). For Kraken 2term, we performed two parameter sweeps, with one focused on minimizer-based subsamplingterm and one focused on hash-based subsamplingterm. The first parameter sweep looked at values for ℓ in the interval [ 25 , 31 ], values for kterm·var in the interval [ ℓ , ℓ + 10], and values for sterm·var in the interval [0, 7]; the second parameter sweep looked at values of ℓ in the interval [ 25 , 31 ], fixed kterm·var = ℓ , values for f in the set {0.125, 0.25, 0.5}, and values for sterm·var in the interval [0, 7]. We also performed a third parameter sweep, focused on translated search (Kraken 2Xterm), where we looked at values for ℓ in the interval [ 11 , 15 ], values for kterm·var in the interval [ ℓ , ℓ + 3], and values for sterm·var in the interval [0, 3].

Each parameter sweep used the strain exclusion data that we previously created to build databases, and we used the same accuracy and timing methods for these databases that we did in the cross-classifier comparison. The results of the first two parameter sweeps, run on nucleotide databases, are provided in Additional file 1 : Table S2, while the results of the third parameter sweep, run on protein databases, are provided in Additional file 1 : Table S9. We note that the parameter sweeps yielded a large number of parameter combinations giving approximately the same, near-optimal levels of accuracy. This suggests performance is not overly sensitive to particular parameter settings.

Evaluation of database sizes of Kraken 1term and Kraken 2term We began by shuffling the reference DNA sequences in our strain exclusion set and recorded the total number of bases in each sequence. We modified Kraken 2term’sterm·var capacity estimator to report an estimate of the number of distinct minimizersterm after each sequence processed, rather than only after all sequences are processed. Finally, we ran the capacity estimator twice on the shuffled genomic data, once with kterm·var = 31, ℓ = 31, sterm·var = 0 (corresponding to Kraken 1term’sterm·var defaults—effectively counting the number of distinct k -mersterm) and again with kterm·var = 35, ℓ = 31, sterm·var = 7 (Kraken 2term’sterm·var defaults).

The size of a Kraken 1term database is a function of the number of distinct k -mersterm in the reference data. If there are X distinct k -mersterm, the size of Kraken 1term’sterm·var database.kdb (sorted list of k -merterm/LCAterm pairs) file will be 1072 + 12 X bytes; the 1072-byte term is the size of the Jellyfish/Krakenterm header data, and 12 bytes are used for each k -merterm/LCAterm pair. The database.idx (minimizerterm·def offset index) file is 8,589,934,608 bytes, a function of Kraken 1term’sterm·var default minimizerterm·def length of 15. The full database size is the sum of the sizes of those two files.

Similarly, the size of a Kraken 2term hash table is a function of the estimate of the number of distinct minimizersterm in the reference data. If there are an estimated Y distinct minimizersterm, Kraken 2term’sterm·var hash table will be 32 + ⌊4 Y /0.7⌋ bytes in size (representing 32 bytes of metadata and using 4 bytes per cell and a load factorterm·var of 0.7).

We used the estimates of the numbers of distinct k -mersterm and distinct minimizersterm to calculate the database sizes of Kraken 1term and Kraken 2term for successively larger subsets of the strain exclusion set. The results of this evaluation are shown in Additional file 2 : Figure S1, with raw data available in Additional file 1 : Table S10.

Reviewing the results when all genomic sequences were added, our results indicate that the number of distinct k -mersterm is approximately 3.1 times the number of distinct minimizersterm for the settings we have selected for Kraken 1term and Kraken 2term. It is not possible to draw a direct relationship between the number of distinct k -mersterm or minimizersterm and the number of sequence bases processed. For example, homology between similar strains and species will cause the number of distinct k -mersterm/minimizersterm to grow slower than the total number of bases. Examining the linear-term coefficients from the database-size expressions (12 X and 4 Y /0.7) indicates a Kraken 2 databaseterm will be approximately 15% of the size of a Kraken 1term database of the same reference data; this is because X ≈ 3.1 Y , and (4/0.7)/(12 × 3.1) = 0.15. When we examine the full reference set, the 15% estimate is consistent with the ratio of Kraken 2term’sterm·var hash table size (10.456 GB) to Kraken 1term’sterm·var database.kdb file size (77.490 − 8.589 = 68.901 GB), which is 10.456/68.901 = 0.152.

Brackenterm experiments on strain exclusion data We first generated Brackenterm metadata from each of the Kraken 1term and Kraken 2term reference libraries used in the strain exclusion experiments. We then used Brackenterm to estimate genus- and species-level abundance from the Kraken 1term and Kraken 2term classification results on the prokaryotic strain exclusion read data. Due to the low sequence similarity between our simulated viral reads and the strain-exclusion reference data, none of the nucleotide search programs exhibited high sensitivity on these reads, including Kraken 1term and Kraken 2term. Such low classification rates prevent Brackenterm from inferring taxonomy for a large proportion of the viral reads. Additionally, the taxonomy for viruses has several examples where species are not grouped by ancestry and lack similarity in both gene organization and genomic sequence [ 38 ]. For these reasons, we chose to exclude the simulated viral reads from our analysis of Brackenterm.

For overall evaluation of the accuracy of Brackenterm in these strain exclusion experiments, we calculated the mean absolute percentage error (MAPE): where Sterm·var x is the estimated number of reads and T x is the true number of reads for taxon x . In this strain exclusion experiment, n = 40, the total number of distinct prokaryotic species and genera in the sample and T x = 1000 for each taxon.

Supplementary information Additional file 1: Table S1. Comparison of accuracy and computational performance. Table S2. Comparison of Kraken 2term with other classifiers, using various parameter values. Table S3. Genomes excluded in strain-exclusion simulation. Table S4. Thread scaling evaluation results. Table S5. Evaluation of FDA-ARGOS sequencing data. Table S6. Sequences used for evaluation of collision rates. Table S7. Minimizer collisionterm evaluation results. Table S8. Hash table collisionterm evaluation results. Table S9. Comparison of Kraken 2Xterm with other classifiers, using various parameter values. Table S10. Database size evaluation results.

Additional file 2: Figure. S1. Estimation of database sizes for Kraken 1term and Kraken 2term as sequences are added to the reference set. Figure S2. Brackenterm performance on strain exclusion simulated prokaryotic data. Figure S3. Examples of compact hash tableterm usage with Kraken 2term. Figure S4. Evaluation of compact hash tableterm error rates as a function of two variables. Figure S5. Evaluation of minimizer collisionterm rates as a function of minimizerterm·def length.

Additional file 3: Review history.

Acknowledgements The authors would like to thank James R. White and Steven Salzberg for the helpful discussions about the manuscript.

Peer review information Barbara Cheifet was the primary editor of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.

Review history The review history is available as Additional file 3 .

Authors’ contributions DEW and BL designed the algorithms for Kraken 2term. DEW developed the Kraken 2term software. DEW, JL, and BL designed the experiments. DEW and JL performed the experiments. DEW, JL, and BL prepared and reviewed the manuscript. All authors read and approved the final manuscript.

Funding BL and DEW were supported by NSF grant IIS-1349906. BL was additionally supported by NIH grant R01-GM118568. JL was supported by NIH grant R35-GM130151.

Availability of data and materials We have made the data for our strain exclusion experiments publicly available for download, including all reference sequencesterm, taxonomy, and simulated read data [ 39 ]. Code to generate all databases from these reference sequencesterm, to generate simulated read data, and to run the comparison of classifiers’ accuracy and performance is also available for public download in a GitHub repository [ 40 ] and via permanent storage at 10.5281/zenodo.3520278 .

Kraken 2term’sterm·var source code is open-source, licensed under the MIT License, and available in a GitHub repository [ 41 ]. The specific version of Kraken 2term evaluated here, version 2.0.8, is also permanently available at 10.5281/zenodo.3520272.

Ethics approval and consent to participate Not applicable.

Consent for publication Not applicable.

Competing interests The authors declare that they have no competing interests.

Footnotes Supplementary information Supplementary information accompanies this paper at 10.1186/s13059-019-1891-0.

References Associated Data Supplementary Materials Additional file 1: Table S1. Comparison of accuracy and computational performance. Table S2. Comparison of Kraken 2term with other classifiers, using various parameter values. Table S3. Genomes excluded in strain-exclusion simulation. Table S4. Thread scaling evaluation results. Table S5. Evaluation of FDA-ARGOS sequencing data. Table S6. Sequences used for evaluation of collision rates. Table S7. Minimizer collisionterm evaluation results. Table S8. Hash table collisionterm evaluation results. Table S9. Comparison of Kraken 2Xterm with other classifiers, using various parameter values. Table S10. Database size evaluation results.

Additional file 2: Figure. S1. Estimation of database sizes for Kraken 1term and Kraken 2term as sequences are added to the reference set. Figure S2. Brackenterm performance on strain exclusion simulated prokaryotic data. Figure S3. Examples of compact hash tableterm usage with Kraken 2term. Figure S4. Evaluation of compact hash tableterm error rates as a function of two variables. Figure S5. Evaluation of minimizer collisionterm rates as a function of minimizerterm·def length.

Additional file 3: Review history.

Data Availability Statement We have made the data for our strain exclusion experiments publicly available for download, including all reference sequencesterm, taxonomy, and simulated read data [ 39 ]. Code to generate all databases from these reference sequencesterm, to generate simulated read data, and to run the comparison of classifiers’ accuracy and performance is also available for public download in a GitHub repository [ 40 ] and via permanent storage at 10.5281/zenodo.3520278 .

Kraken 2term’sterm·var source code is open-source, licensed under the MIT License, and available in a GitHub repository [ 41 ]. The specific version of Kraken 2term evaluated here, version 2.0.8, is also permanently available at 10.5281/zenodo.3520272.