MetaPalette: a k-mer Painting Approach for Metagenomic Taxonomic Profiling and Quantification of Novel Strain Variation

Taxonomic profiling is a challenging first step when analyzing a metagenomic sample. This work presents a method that facilitates fine-scale characterization of the presence, abundance, and evolutionary relatedness of organisms present in a given sample but absent from the training database. We calculate a “k-mer palette” which summarizes the information from all reads, not just those in conserved genes or containing taxon-specific markers. The compositions of palettes are easy to model, allowing rapid inference of community composition. In addition to providing strain-level information where applicable, our approach provides taxonomic profiles that are more accurate than those of competing methods.

M etagenomics is a developing field used to characterize the organismal composition of microbial communities in environmental or clinical samples (1). A key step in most metagenomic analyses is to identify the organisms in the sample and their relative frequencies. A wide variety of different algorithms have been developed for this purpose.
Most approaches, including the one described here, are based on relating sequenced reads to reference organism genome sequences. Conceptually, the aim of these approaches is to place the organisms in the sample on a "tree of life" that has been defined in advance. In practice, the available reference organisms are extremely unevenly scattered through the true tree of life. Many medically important branches, such as enterobacteria, are relatively well sampled, with many strains from the same species, while there are entire phyla of unculturable organisms that are unrepresented (2)(3)(4)(5).
A further difficulty, both in theory and in practice, is that a fully resolved tree of life cannot be established, even from complete reference genomes. At the scale of individual species, homologous recombination scrambles variation so that a tree is not necessarily an appropriate representation of organismal relationships, while more distant phylogenetic relationships can be difficult to estimate due to the various technical challenges of reconstructing ancient evolutionary events (6)(7)(8)(9)(10).
Based on these practical considerations, an effective metagenomic method should both identify the closest organism or set of organisms in the reference set and also estimate the genetic difference between the closest reference(s) and the organism present in the sample. The method should work both if the closest neighbor is a distant member of the same phylum and if there are multiple strains within the species in question. Fine-scale classification is important because the detailed knowledge we have, of E. coli, for example, shows that organisms from the same species can have entirely different ecologies and phenotypic effects on their host (11).
Given these difficulties, a number of different approaches are taken to characterize metagenomic samples. A commonly used approach is to first place individual reads onto a tree constructed for a particular set of genes and then attempt to sum the phylogenetic information across the reads (12)(13)(14)(15). Phylogenetic analysis of each read can be computationally challenging for large datasets, and individual reads can often only be placed inaccurately. It is challenging to appropriately represent this uncertainty in later stages of the analysis. These approaches also break down if a tree is not a good representation of relationships among organisms, e.g., within species. Furthermore, while utilizing specific genes (so-called marker genes) can increase computational efficiency, this approach throws away a considerable amount of information from sequences that do not align to the marker genes. As a result of these issues, these methods are typically accurate for genus level or higher classification but not for fine-scale classification.
Another approach identifies features that are characteristic of particular organisms, such as the frequency of k-mers (16)(17)(18). These features are used either for taxonomic binning of individual reads or in order to compute the overall composition. Depending on the k-mer size utilized, these methods either are suitable only for higher-level phylogenetic analysis (for small k-mers) or are highly dependent on the training database utilized (for larger k-mers). In either case, no existing method using this approach can accurately detect and classify organisms that are highly diverged from ones in the training database, and the existing methods still struggle with quantifying strain-level variation. Using longer k-mers allows for higher specificity, but using k-mers that are unique to specific taxa in the reference data set (as in references 16 and 17) ignores a great deal of information about evolutionary relatedness provided by other k-mers. It also makes the approach highly dependent on the specific composition of the reference data set. We argue that utilizing all k-mers in a reference database and multiple k-mer sizes allows the modeling of the k-mer signatures of organisms absent from a given training database.
In this paper, we present an approach based on defining a "palette" for each reference organism. Specifically, we count the number of k-mers found in the sample DNA that are present in each reference organism. Our approach thus uses all k-mers of a particular length in the reference data set, while discarding the specific information provided by matches of individual k-mers. This is similar in spirit to the so-called pseudoalignment approach in reference 19, except that here, we use k-mer counts of the entire sample, not of individual reads whose origins may be ambiguous. We model these palettes using a simple linear mixture model which includes both the reference organisms and "hypothetical organisms" of different degrees of genetic relatedness to the reference organisms. The algorithm is called MetaPalette, and the outputs of the algorithm are demonstrated in Fig. 1.
We first introduce the concept of a "common k-mer matrix" and demonstrate how utilizing multiple k-mer sizes allows accurate quantification of evolutionary relatedness. We then develop a mixture modeling procedure that utilizes this information to taxonomically profile a metagenomic sample and indicate the evolutionary relatedness of novel organisms. Evidence on simulated and real data is given that demonstrates that this approach can accurately capture strain-level variation, and we then benchmark this approach against other, commonly utilized metagenomic-profiling techniques.

MetaPalette design.
We begin by introducing the basic analytic approach. For a more detailed description, see Materials and Methods.
(i) Common k-mer training matrix. To quantify the similarity of two genomes, we count (with multiplicity) the fraction of each genome's k-mers that are in common with the other. Rigorous mathematical definitions of this and other quantities are contained in Materials and Methods. This quantity, denoted for pckm k ͑·,·͒ "percentage of common k-mers," is similar to the well-known Jaccard index (20) except that, among other differences, pckm k ͑·,·͒ is not symmetric but does incorporate the counts of k-mers, not just their occurrence.
When given a set of genomes (i.e., a training database), a pairwise similarity matrix can be formed using the equation A i,j ͑k͒ ϭ pckm k ͑g i ,g j ͒ for g i and g j training genomes. The column vector pckm k ͑·,g j ͒ can be thought of as a palette, representing the particular k-mer profile of g i in relation to those of other genomes. We call each of these matrices a "common k-mer matrix." These matrices reflect the relatedness of the training genomes based on k-mer similarity. For larger k-mer sizes, one can clearly extract taxonomic information from these matrices, as shown in Fig. 2.
Beyond genus-level variation, strain-level variation can be captured through these common k-mer matrices. For example, using all the strains of the species Burkholderia multivorans accessible via NCBI, we formed a neighbor-joining tree using the average of the 30-mer and 50-mer common k-mer matrices. This tree, shown in Fig. 3, demonstrates how the common k-mer matrices can capture variations among these strains.  The entries of A ͑k͒ can be calculated in a computationally efficient manner. We take the approach of forming bloom count filters (using Jellyfish [21]) for each of the training genomes and then counting the common k-mers using a simple Cϩϩ program based on heap data structures.
(ii) Modeling related organisms. To model the k-mer counts for organisms with various degrees of relation from the training database, we take advantage of the differing behavior of pckm k ͑·,·͒ as a function of k for closely related organisms and distantly related organisms. In particular, the percentage of common k-mers, pckm k ͑·,·͒, decays more slowly as a function of k for closely related organisms than for distantly related ones. This is consistent with the intuitive idea that, for example, two organisms from different phyla will have similar percentages of shared 1-mers but very few common 50-mers. Conversely, two closely related strains will have both high percentages of shared 1-mers and high percentages of shared 50-mers. This is demonstrated in Fig. 4a. This property means that using more than one k-mer length should in principle allow us to distinguish between having an organism that is identical to a training organism at a low frequency and having an organism that is distantly related to all training organisms but present in the sample at a higher frequency.
We focus on two particular k-mer sizes, k ϭ 30 and k ϭ 50, due to the predictability of pckm k for these k-mer sizes. Indeed, using 6,914 whole bacterial genomes downloaded from a variety of publicly accessible repositories (via RepoPhlAn [https://bitbucket.org/nsegata/repophlan]), we observed that the percentage of shared 30-mers can be predicted from the percentage of shared 50-mers (Fig. 4b). A degree 3 polynomial was used (as it resulted in the lowest root mean square error [RMSE] and R 2 values, which did not improve for higher-degree polynomials). Namely, we observed that for the polynomial p͑x͒ ϭ Ϫ .5141x. 3 For k-mer lengths substantially shorter than 30, the behavior of pckm k is more variable-for example, because of convergence of sequence composition between distantly related organisms. On the other hand, k-mers much larger than 50 are increasingly time consuming to compute and are likely to be more sensitive to sequencing error and other technical artifacts.
We can augment the matrices A ͑k͒ with columns that represent hypothetical organisms which are related by different degrees to the reference organism. For a given organism with genome g i , if we wish to include a hypothetical organism h that is 90% similar to genome g i in its 30-mers, we can round down each entry of the column vector pckm 30 ͑·,g i ͒ to be no more than 0.90. Call this vector pckm 30 ͑·,h͒. The entries below 90% do not need to be changed, since we assume that the hypothetical organism has the same patterns of k-mer sharing to more distantly related "outgroup" taxa as to the reference organism. We model the 50-mer similarity by setting pckm 50 ͑·,h͒ ϭ p͑pckm 30 ͑···,g i ͒͒ for the previously defined polynomial p. Adding these vectors to A ͑30͒ and A ͑50͒ effectively adds a hypothetical organism that has a common k-mer signature 90% similar to that of genome g i . We then repeat this procedure for all training genomes g i and for similarities ranging from 90%, 80%,. . .,10% and append these columns to A ͑30͒ and A ͑50͒ .
(iii) Sample k-mer signature. Given a metagenomic sample, we form two vectors, y ͑30͒ and y ͑50͒ , consisting of the total counts in the sample of the 30-mers and 50-mers shared with the training organisms. In "Mathematical formulation" in Materials and Methods, we show that these vectors are linearly related to the organism abundances via the common k-mer matrices A ͑30͒ and A ͑50͒ .
Note that in forming y ͑k͒ , we count the k-mers in the entire sample, not of the individual reads. This allows for a very computationally efficient approach: as the training genomes typically have low error, their k-mers can be efficiently stored in de Bruijn graphs (formed using Bcalm [22]). We can then query the bloom count filter formed from the sample in a highly parallel fashion.
(iv) Sparsity-promoting optimization procedure. After forming y ͑k͒ , we note that some of the entries y i ͑k͒ may be nonzero due not to the presence of organism i in the sample but to the fact that there exists an organism j that shares portions of its genome with organism i. Since A i,j ͑k͒ represents the "overlap" of these two organisms, we can deconvolute this linear mixture relationship by solving the equation A ͑k͒ x ϭ y ͑k͒ for x, the vector of organism abundances. However, after having augmented A ͑k͒ with the hypothetical organisms, this system of equations is underdetermined (10 times more columns than rows). We can employ a sparsity-promoting optimization procedure to infer the most parsimonious x consistent with the equation A ͑k͒ x ϭ y ͑k͒ for k ϭ 30, 50. This procedure, first introduced in reference 23 and proven correct in reference 24, is detailed in "Optimization procedure" in Materials and Methods.
(v) Inferring taxonomy. The abundances of the hypothetical organisms are then mapped back onto the taxonomy (for the output taxonomic profile) or the neighborjoining tree formed from A ͑k͒ (for the output strain variation figures), utilizing a least-common-ancestor approach detailed in "Inferring taxonomy" in Materials and Methods.
Quantification of strain-level variation. We demonstrate in two ways that the inclusion of the hypothetical organisms allows the inference of strain-level variation. First, we spike novel organisms into a mock metagenomic community and show that MetaPalette can accurately predict their presence. Second, we utilize a real metagenomic soil sample to give evidence for a novel strain that MetaPalette predicts.
(i) HMP mock community. We first formed the common k-mer matrices A ͑k͒ using 31 strains of Lysinibacillus sphaericus. We then used Grinder (25) to simulate a data set consisting of two novel strains (not included in the training database). These reads were then spiked into the Human Microbiome Project (HMP) mock even community (añ 6.6 million-read metagenome consisting of 22 select organisms sampled using an Illumina GA-II sequencer; NCBI accession number SRR172902). The output of Meta-Palette is shown in Fig. 5, demonstrating the ability of the method to correctly infer the presence of organisms absent from the training data.
Decreasing the number or changing the identity of the training organisms does not impede the method. As shown in Fig. 6, 50,000 simulated reads from the species Providencia alcalifaciens were again spiked into the HMP mock even community, and the inferred abundance was again placed optimally on the neighbor-joining tree. See Fig. 11 to 14 for a variety of such figures spanning all domains of life. These results provide evidence that MetaPalette can correctly infer the presence of organisms related to but absent from the training database.
(ii) Metagenomic soil sample. To assess MetaPalette on a real metagenomic sample, we utilized the Iowa prairie metagenomic sample from reference 26 (corresponding to MG-RAST project accession number 6377). After running MetaPalette on a subset of these data (metagenome accession number 4539594.3), the taxonomic profile that was returned predicted the presence of the genus Bradyrhizobium. Generating the tree plot on a subset of this genus resulted in, among others, a prediction of a novel organism in the clade defined by strains of Bradyrhizobium valentinum (Fig. 7a). To verify this, we aligned the entire soil metagenome to the reference genome of the strain B. valentinum LmjM3 using Bowtie2 with -very-sensitive-local settings (27) and extracted the aligned reads. Interestingly, 0.29% of the reads aligned, while the MetaPalette-predicted abundance for this putative novel organism of interest was 0.33%. The depths of coverage of the extracted reads are pictured in Fig. 7b; the mean depth was 74.3ϫ.
To assess the evolutionary relatedness of this predicted organism, we utilized the B. valentinum LmjM3 nifH gene sequence (NCBI accession number KF806461), which was used in the work reported in reference 28, along with other genes, to determine the taxonomy of B. valentinum LmjM3. Aligning the extracted reads to nifH resulted in a mean depth of coverage of 22ϫ. We collapsed the aligned reads (via a majority vote) in regions of coverage of at least 22ϫ and called this the maximum-likelihood sequence. We then performed a multiple-sequence alignment of this sequence along with the nifH sequences of 20 other organisms closely related to B. valentinum. The topology of the bootstrap consensus neighbor-joining tree is pictured in Fig. 7c and shows that the maximum-likelihood sequence is placed at the same location as was predicted by MetaPalette. While this is not enough evidence to unequivocally claim the existence of a novel strain in this sample, this gives support that MetaPalette correctly inferred the abundance and placement, as shown in Fig. 7a, of a potentially novel strain in the clade defined by strains of B. valentinum.
Comparison to other metagenomic profiling methods. To facilitate an objective comparison with other methods with minimal "author bias," we utilized the same data and metrics used by other authors in a recent paper evaluating metagenomics methods (29). This allowed comparison to the following algorithms: CLARK (i) Training data. Each of the methods was trained using the default recommended databases. We trained our method using 6,914 whole-genome sequences and assemblies obtained from various public repositories via RepoPhlAn (https://bitbucket. org/nsegata/repophlan). The training procedure for MetaPalette on these 6,914 organisms took a total of approximately 7 h on a 48-core server.
(ii) Testing data. The testing data consisted of 6 samples and are fully explained in reference 29 (in Methods), but we briefly summarize them here. Three replicates were formed from two different distributions of over 900 different genomes spanning the tree of life (including Eukaryote genomes). Included in each test sample were shuffled/ randomized genomes (not meant to be assigned to any known taxon), as well as sequences from the genome of Leptospira interrogans that were evolved using Rose (40) to simulate novelty. Error profiles were based on those of 6 real soil metagenomic samples sequenced using an Illumina HiSeq 2000. Each of the resulting test samples contains between 27 and 37 million read pairs.
(iii) Error metrics. We utilized the same divergence error metric as was used in reference 29; that is, for x i , representing the true frequency of taxon i in the sample, and x i ‫ء‬ , representing the predicted frequency of taxon i for a given method, where the summation is over those indices such that X i Ͼ 0 and x i ‫ء‬ Ͼ 0. Since this error metric does not take into consideration the number of spurious assignments (that is, taxa predicted by a method to be in a sample but not actually present), we also use the number of false positives (FP) at a given taxonomic rank, as follows: FP ϭ |{i : x i ‫ء‬ Ͼ 0 and x i ϭ 0}|· (iv) Comparison results. Each method was run using the default parameters. For each method, we averaged the divergence error metric over all the test samples at the genus level (Fig. 8a). Furthermore, we selected a number of the more accurate methods and averaged the number of false positives over all the test samples at the phylum level (Fig. 8b). The results in Fig. 8a and b clearly show the competitive nature of MetaPalette, as it has the lowest error using both metrics. However, when comparing to other methods, one should be careful of their intended use. For example, taxator-tk is intended to be used on an assembled metagenome (and here unassembled reads were used), and QIIME only uses the 16S rRNA sequences in a sample. Furthermore, most of these methods assign individual reads and then summarize them to obtain a taxonomic profile, while our method only profiles the entire sample and returns relative proportions of organisms. Figure 9 shows the execution time of each of the methods (on a log scale, obtained from reference 29), further showing the competitive nature of MetaPalette.
Software and pretrained data. (i) Software. The source code for MetaPalette, along with installation instructions and directions, is accessible at https://github.com/ dkoslicki/MetaPalette. MetaPalette is written primarily in python and accepts input reads in uncompressed fasta or fastq format, as well as compressed fasta/fastq using bzip2 and gzip. For fastq input, optional parameters can be given to specify counting k-mers only above a certain quality score (Phred), thereby attenuating the negative impact of sequencing error in the correct inference of relative abundances. The output taxonomic profile is compliant with the Bioboxes profiling format, version 0.9, found at https://github.com/bioboxes/rfc/tree/master/data-format. Python scripts are also in-cluded to aid in downloading data, forming custom databases, and creating the appropriate taxonomy files.
A preliminary version of this software was submitted to Critical Assessment of Metagenomic Interpretation (CAMI: http://www.cami-challenge.org/) under the name CommonKmers. However, since significant changes have been made since that point, we strongly recommend using the current MetaPalette software instead.

DISCUSSION
We have described a fast, flexible, and accurate method for estimating the taxonomic composition of organisms which is based on reconstructing a k-mer-based profile of a sample. Each reference organism has a k-mer "palette," and we fit the sample as a mixture of different palettes, both of the reference organisms and organisms absent from the training data at various degrees of relatedness to the training database. Our approach is in part inspired by the chromosome-painting method used to deduce fine-scale population structure in human genetics (45,46), which is also based on mixture modeling of palettes. A particular advantage of MetaPalette over other metagenomic profiling methods is that MetaPalette provides an indication of how related the organisms in a given sample are to the closest matching organisms of the training database, whether they are within the same species or distantly related organisms from the same phyla.
Furthermore, the standard approach to summarizing composition information has been to place organisms at different taxonomic levels. We produce a standard taxonomic profile which we have shown to be more accurate than that produced by other methods. This fixed-rank approach is sensible at the genus level and above but omits fine-scale information. Hence, for branches of the tree of life that are well represented in the training database, we can also output a phylogenetic tree giving detailed information on how the sampled taxa relate to the organisms in the training database ( For many applications, it is of interest to understand which individual reads belong to which organisms (1). A principled approach to this problem is to first estimate the overall composition of the sample, using MetaPalette or an equivalent, and then to assign individual reads conditional on the overall assignment. This represents a promising avenue for future methodological development.

Mathematical formulation.
We include here all rigorous mathematical definitions of the quantities discussed in the main text.
Given the alphabet A ϭ ͕A,C,T,G͖, let A n denote the set of all words v of length |v| ϭ n on A, and let A ‫ء‬ ϭ ഫ nՆ0 A n be the set of all finite words on A. Hence, words containing non-ACTG characters are ignored. Let D ϭ ͕g 1 ,···,g M ͖ be a database of genomic sequences g j ʦ A ‫ء‬ and let S ϭ ͕s 1 ,···,s N ͖ be a set of sample sequences (the reads to be classified). For notational simplicity, assume that the read length is fixed; i.e., for all t ϭ 1,ѧ,N |s t | ϭ r. Fix a k-mer size and endow A k ϭ ͕v 1 ,···,v 4 k͖ with the lexicographic order. Let occ v ͑w͒ represent the number of occurrences (with overlap) of the subword v in the word w. That is, for w,v ʦ A n , let For a fixed k-mer size and two genomes, g i and g j , we calculate the number of k-mers in genome j common to both g i and g j . That is, the ͑i,j͒ th entry of the common k-mer training matrix A ͑k͒ is as follows: Refer to the entries of the common k-mer matrix as pckm k ͑g i ,g j ͒ ϭ A i,j ͑k͒ . Let s i ʚ g j denote the relationship that read s i was derived from genome g j . We represent the taxonomic profile of the sample S by the probability vector x as follows: where is the indicator function. Now let the measurement vector y be given by the probability vector We assume that the reads s t are uniformly randomly selected from the genome g j . Then, for w ʦ A k , let ‫͑ސ‬w|g j ͒ be the probability that k-mer w is found in genome g j . Then, we have that the proportion of k-mers w in the sample is similar to the proportion of the appearance of w in the genome g j when weighted by the relative abundance of the genome g i in the sample, as follows: We then calculate Equation 10 is justified since if w SW g j ͑k͒, then occ g i ͑w͒ ϭ 0. For computational reasons, we make the assumption in equation 11 that SW g i ͑k͒പSW g j ͑k͒പSW S ͑k͒ ϭ SW g i ͑k͒പSW g j ͑k͒. However, this assumption can be mitigated by adding hypothetical organisms (see "Hypothetical organisms" below). Our assumptions imply that We will try to recover the vector x satisfying x j Ն 0 for all j ϭ 1,···,M from equation 13. Further improvements. A few further improvements are possible, but not pursued here. Namely, we could use just the k-mers that are actually in the sample to form the training matrix. That is, use SW g i ͑k͒പSW g j ͑k͒പSW S ͑k͒ in the formation of A ͑k͒ , as follows: The disadvantage of this is that the (slow) training step would need to be rerun for each sample. For a second improvement, we could make the approximation in equation 5 more delicate by incorporating the coverage Koslicki and Falush Volume 1 Issue 3 e00020-16 msystems.asm.org 12 So A ͑k͒ would have the form: This effectively multiplies column j of A ͑k͒ by the percent coverage of genome j. Finally, in equation 16, we could put a weighting factor that represents how unique a k-mer is to the genome in question. This would down-weight k-mers shared among many diverse genomes and up-weight those unique to certain strains/species/genera/etc. Hypothetical organisms. To simulate an organism that is, say, 90% related to a database genome g i , we augment the common k-mer matrix A ͑k͒ with a column derived by rounding down the entries of the column vector ͕A i,j ͑k͒ ͖ iϭ1,···,M that are above 90%. Two k-mer sizes are needed to form the hypothetical organism's common k-mer matrices. For the first k-mer size, k 1 , we define A ͑k 1 ͒,h for a fixed number of hypothetical bins h ʦ ͕0.9,0.8,ѧ,0.1͖, where For the second k-mer size, k 2 , using the polynomial p͑x͒ ϭ Ϫ .5141x. 3 ϩ 1.0932x. 2 ϩ 0.3824x, we define Instead of thresholding, as we did here, one can imagine other scalings obtained from studying the relationship between a given taxonomy and the common k-mer matrix A ͑k͒ . In particular, to deal with differing rates of evolution in the tree of life, a fruitful area of future investigation would be to modify the polynomial p depending on the taxonomy of the organisms under consideration. Optimization procedure. We choose two k-mer sizes to be k 1 ϭ 30 and k 1 ϭ 50, as this seems to give a good trade-off between reconstruction fidelity and computational performance. We then collect the common k-mer matrix and hypothetical matrices blockwise into the 2|D| ϫ 10|D| size matrix as follows: A ϭ ͫ A (30) , A (30),0.9 , ··· , A (30),0.1 A (50) , A (50),0.9 , ··· , A (50),0.1 ͬ Collect also the k-mer sample vectors y ͑k͒ : y ϭ ͫ y (30) y (50) ͬ (20) The problem at hand is then to reconstruct the phylogenetic profile x by solving the linear system For the output taxonomic profile, a hybrid of the fixed rank and LCA approaches can increase sensitivity or specificity. We thus include three options: the default option is the LCA approach, while the sensitive and specific options are various hybrids of the two methods.
Sequence analysis details. For the tree shown in Fig. 7c, the evolutionary history was inferred using the neighbor-joining method (49). The bootstrap consensus tree, inferred from 500 replicates, is taken to represent the evolutionary history of the taxa analyzed (50). Branches corresponding to partitions reproduced in less than 50% of the bootstrap replicates are collapsed. The percentages of replicate trees in which the associated taxa clustered together in the bootstrap test (500 replicates) are shown next to the branches. The analysis involved 21 nucleotide sequences. All positions with less than 95% site coverage were eliminated. That is, less than 5% alignment gaps, missing data, and ambiguous bases were allowed at any position. There were a total of 652 positions in the final data set. Evolutionary analyses were conducted in MEGA6 (51).   10 depicts a tree constructed using the same method as just described but with evolutionary distances computed using the maximum-composite-likelihood method (52). The unit is the number of base substitutions per site.
Additional figures. We provide here a number of additional output figures from MetaPalette ( Fig. 11 to 14) to demonstrate that the ability to correctly infer the presence of organisms related to but absent from the training database is not dependent on the particular kingdom/phyla/etc. used. Unless otherwise noted, a total of 50,000 simulated reads from the novel organisms were spiked into the HMP mock even community. Figures are included for Bacteria, Archaea, Eukaryota, and viruses.