{'fraction_non_alphanumeric': 0.03009512827904295, 'fraction_numerical': 0.020121072355145576, 'mean_word_length': 5.016649323621228, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pseudoalignment for metagenomic read assignment

We explore connections between metagenomic read assignment and the quantification of transcripts from RNA-Seq data. In particular, we show that the recent idea of pseudoalignment introduced in the RNA-Seq context is suitable in the metagenomics setting. When coupled with the Expectation-Maximization (EM) algorithm, reads can be assigned far more accurately and quickly than is currently possible with state of the art software.


Introduction
The analysis of microbial communities via whole-genome shotgun sequencing has led to exceptional bioinformatics challenges [Chen and Pachter, 2005] that remain largely unsolved [Scholz et al., 2012]. Most of these challenges can be characterized as "de novo" bioinformatics problems: they involve assembly of sequences, binning of reads, and annotation of genes directly from sequenced reads. The emphasis on de novo methods a decade ago was the result of a paucity of sequenced reference microbial and archaeal genomes at the time. However this has begun to change in recent years [Land et al., 2015]. As sequencing costs have plummeted, the number of fully sequenced genomes has increased dramatically, and while a large swath of the microbial world remains uncharacterized, there are now thousands of "reference quality" genomes suitable for the application of reference-based methods.
One of the fundamental metagenomics problems that is amenable to reference-based analysis is that of "sequence classification" or "read assignment". This is the problem of assigning sequenced reads to taxa. The MEGAN program [Huson et al., 2007] was one of the first reference-based read assignment programs and was published shortly after sequencing-by-synthesis methods started to become mainstream. It provided a phylogenetic context to mapped reads by assigning reads to the lowest taxonomic level at which they could be uniquely aligned, and became popular in part because of a powerful accompanying visualization toolkit. One of the drawbacks of MEGAN was that its approach to assigning ambigously mapping reads limited its application to quantification of individual strains, an issue which was addressed in a number of subsequent papers, culminating in GRAMMy [Xia et al., 2011] and GASiC [Lindner and Renard, 2013], which were the first to statistically assign ambigously mapped reads to individual strains. Unfortunately, these approaches all relied on read alignment, a computational problem that is particularly difficult in the metagenomic setting where reference genome databases are large and read sets gigantic.
In a breakthrough publication in 2014 [Wood and Salzberg, 2014] it was shown that it is possible to greatly accelerate read assignment utilizing fast k-mer hashing to circumvent the need for read alignment. An implementation called Kraken was used to show that analyses that previously took hours were tractable in minutes, and the removal of the read alignment step greatly simplified workflows and storage requirements. However the Kraken speed came at a cost. An examination of the Kraken algorithm and output reveals that the method takes a step back from GRAMMy and GASiC by discarding statistical assignment of reads at the strain level in favor of direct taxonomic assignment as in MEGAN (notably according to [Lindgreen et al., 2015], Kraken is more accurate than MEGAN although, as we'll see, it is not as accurate as GASiC). A natural question to ask is whether the strengths of Kraken and GASiC can be combined, i.e. whether it is possible to leverage fast k-mer based hashing to map reads not at the taxonomic but at the strain level, while assigning the resulting ambigously mapped reads using a statistical framework that allows for probabilistic assignment of reads.
To answer this question we turned to RNA-Seq [Cloonan et al., 2008, Lister et al., 2008, Nagalakshmi et al., 2008, Mortazavi et al., 2008, an experiment for which there has been extensive methods development that we hypothesized could be adapted and applied to metagenomics. Many of the challenges of metagenomic quantification translate to problems in RNA-Seq via a dictionary that replaces genome targets with transcript targets. For example, ambigously mapped genomic reads that are difficult to resolve at the strain level in the metagenomics setting are analogous to reads that are difficult to assign to specific isoforms in RNA-Seq. Statistical questions at the heart of "comparative metagenomics" [Huson et al., 2009, Rodriguez-Brito et al., 2006, Tringe et al., 2005 are analogous to the statistical problems in differential expression analysis. In fact, the only significant differences between metagenomics and RNA-Seq are that genome sequences are much larger than transcripts and reference databases are less complete. These differences have engineering implications, but statistically and computationally, metagenomics and transcriptomics are very much the same.
In this paper we show that technology transfer from RNA-Seq to metagenomics makes it possible to perform read assignment both rapidly and accurately. Specifically, we show that it is possible to accurately assign reads at the strain level using a fast k-mer based approach that goes beyond the hashing of Kraken and takes advantage of the principle of pseudoalignment [Bray et al., 2015]. The idea of pseudoalignment originates with RNA-Seq, where it was developed to take advantage of the fact that the sufficient statistics for RNA-Seq quantification are assignments of reads to transcripts rather than their alignments. The same applies in the metagenomics setting, and we show that just as in RNA-Seq, application of the EM algorithm to "equivalence classes" [Nicolae et al., 2011] allows for accurate statistical resolution of mapping ambiguities. Using a published simulated dataset [Mende et al., 2012] and an implementation of pseudoalignment coupled to the EM algorithm in kallisto [Bray et al., 2015], we demonstrate significant accuracy and performance improvements in comparison to state of the art programs.

Results
To test the hypothesis that RNA-Seq quantification methods can be applied in the metagenomics setting we began by examining the performance of eXpress, a program that implements a streaming EM algorithm for RNA-Seq read assignment from alignments, on simulated data [Roberts and Pachter, 2013]. We chose eXpress because it utilizes traditional read alignments directly to a transcriptome but is more memory efficient than other approaches (e.g. RSEM [Li and Dewey, 2011]) and therefore more suitable in the metagenomics setting. Other RNA-Seq quantification tools such as Cufflinks [Trapnell et al., 2010] were not suitable for our needs because of their dependence on read alignments to genomes and not transcriptomes, a requirement that does not translate easily to the metagenomics setting.
To test eXpress we aligned a simulated dataset of Illumina-like reads from 100 microbial genomes to a reference database containing only those genomes, allowing us to compare results to a ground truth (the Illumina100 data) [Mende et al., 2012]. We began by comparing eXpress to GASiC, which also utilizes read alignments for read assignment. The results are shown in Table 1. We found that eXpress outperforms GA-SiC at the exact genome, species, genus, and phylum levels, which we believe is because the statistical model of eXpress takes into account data-dependent read error profiles in assigning reads.
A major problem with GASiC and eXpress is that the alignments they require are slow to generate. The alignments, made with Bowtie2 [Langmead and Salzberg, 2012], took days. As reported in [Wood and Salzberg, 2014], Kraken was much faster on the data, taking only 22 minutes 38s. We also tested CLARK [Ounit et al., 2015], another recently published k-mer based assignment tool and, in agreement with the benchmarks in [Lindgreen et al., 2015], we found it to be slightly faster taking 20 minutes 30s. However, as seen in Table 1, both Kraken and CLARK have significantly worse performance than both GASiC and eXpress (in concordance with [Lindgreen et al., 2015] but in contradiction to [Ounit et al., 2015], we found that Kraken is more accurate than CLARK).
We next turned to a comparison of kallisto with Kraken and CLARK using the Illu-mina100 simulation (i100) but utilizing a larger and more realistic reference database of 1,958 genomes from [Martin et al., 2012]. The results, shown in Tables 1,2 and in Figure 1 (where the database is called "i100+Martin") show that kallisto is much more accurate than Kraken and CLARK at all four taxonomic levels tested. The performance of kallisto at the exact genome level is significantly better than Kraken at the genus level. Notably, Kraken's performance at the exact genome level (AVGRE 18.67 and RRMSE 38.26) is too poor to be of practical use.
The running time of kallisto is faster than Kraken with the i100 dataset (5m55s vs. 22m38s using a single core) but slower on the i100+Martin (50m55s vs. 27m3s). However kallisto index building is slightly faster than Kraken, and when utilizing multiple threads for pseudoalignment, kallisto's running time is negligible.
While the results of kallisto on the i100 and i100+Martin databases were convincing, we decided to examine its performance in the case of missing strains, a situation that is commonplace in metagenomic analyses. We examined two different scenarios: the effect  Table 3. Normalized count based classification accuracy after removing single genomes from index.
As expected, Table 3 shows that the performance of both kallisto and Kraken degrades with removal of strains from the reference database. However, what is interesting is that in the case of Bacillus, read assignment is still possible at the genus level as reflected in improved performance. In the case of Listeria, the absence of neighboring species means that the reads cannot be assigned. Remarkably, even with a strain completely missing from the reference database, the performance of kallisto is still higher than that of Kraken with all genomes present (compare to Table 1).

Methods
3.1. Illumina100 dataset. We tested kallisto and alternate programs on a set of simulated reads published in [Mende et al., 2012]. The Illumina100 dataset consists of 53.33 million 75bp reads, simulated by the iMESSi metagenomic simulator using an Illumina error model. The reads were simulated from a set of 100 unique bacterial genomes. The set is of genomes from 85 different species and 63 different genuses, over a range of abundances from 0.86% to 2.2%.
Reads were trimmed with the program Trimmomatic (version 0.32) [Bolger et al., 2014] to a minimum length of 40bp, using its adaptive trimming algorithm MAXINFO with a target length of 40 and default strictness. 40 reads were dropped due to quality issues.

Taxonomic identification.
We analyzed each program's output at four taxonomic ranks: phylum, genus, species, and "exact genome" level. The latter tests the abundance estimation of the actual Illumina100 genomes, which are a combination of strains and substrains and thus aren't taxonomically well defined. The other three ranks are as assigned by NCBI's Taxonomy Database, as of November, 2015.
3.3. Count estimation accuracy calculation. Using a simulated dataset with known abundances allowed us to benchmark programs by comparing program outputs with true values for each genome. While kallisto is able to output length-corrected individual genome abundances, most of the programs we compared with only output counts, so for consistency we analyzed the accuracy of assigned or estimated counts for each program. We normalized the estimated counts by the percent of assigned reads in order to be able to compare relative count estimates between programs.
We primarily used the error measures AVGRE (Average Relative Error), which computes the mean of the difference between truth and estimate, and RRMSE (Relative Root Mean Square Error), which computes the root mean square average of the difference between truth and estimate, to judge the accuracy of our estimates. Formally, with n true genomes/species/genera/phyla, true counts τ i (1 ≤ i ≤ n) and estimated counts t i at the rank, and A aligned reads out of T total reads we computed In addition, for comparison of programs using the metrics in [Wood and Salzberg, 2014], we calculated the assignment aggregate precision and sensitivity of kallisto, Kraken and CLARK. Instead of examining the results of assignments of individual reads to specific species or genus, we relaxed the benchmark to instead measure sensitivity and precision based on aggregate counts at taxonomic ranks. Mimicking the computations of [Wood and Salzberg, 2014], we computed the aggregate sensitivity at a given rank R by calculating (# of counts correctly assigned at rank R)/(# of counts assigned at rank R).
Aggregate precision was calculated as (# of counts correctly assigned at or below rank R)/(# of counts assigned at or below rank R + # of reads incorrectly assigned above rank R). Counts at a rank were considered to be correctly assigned when they were less than or equal to the true counts. Our level of granularity in assessing assignment reflects the quantity of interest in reference-based metagenomics, namely the accuracy of the aggregate number of counts assigned to individual members of a rank instead of the correctness of each individual read assignment.
The scripts used to compile the results are available at https://github.com/pachterlab/metakallisto 3.4. Reference Genome Database. In addition to aligning the Illumina100 reads against their originating genomes, we tested the more realistic case of aligning against a large bacterial database. In order to have a consistent, reproducible test that won't change as new bacteria are sequenced, we used a frozen published database for all our i100+Martin tests [Martin et al., 2012]. The database contains 1,751 bacterial genomes and plasmids spread over 1,253 species, in addition to 131 Archaea genomes. The published database also contained 3,683 viral genomes and 326 lower eukaryote genomes, but these were discarded before indexing. For compatibility with Kraken and CLARK, we also discarded 6 genomes that were lacking sequence GI numbers, and one header that didn't contain any sequence information (GI 308222630). This database was then combined with the Illumina100 source genomes for a total of 1,958 individual genomes.

Discussion
The idea of translating RNA-Seq methodology to and from metagenomics was, to our knowledge, first proposed in [Paulson et al., 2013] where statistical methods for identifying differential abundances in microbial marker genes were developed. In that paper, there were comparisons between the proposed metagenomics method and RNA-Seq differential analysis methods implemented in DESeq [Anders and Huber, 2010] and edgeR [Robinson et al., 2010]. Notably, the central idea of the paper, the specific consideration of zero inflated distributions to account for undersampling, is also used in single cell expression analysis [McDavid et al., 2013].
Our results show that RNA-Seq methods for quantification are also applicable in the metagenomics setting, and our results with kallisto demonstrate that it is possible to accurately and rapidly quantify the abundance of individual strains. With a few exceptions, e.g. [Bradley et al., 2015], most metagenomic analyses have focused on higher taxonomy, a point highlighted in the recent benchmarking paper [Lindgreen et al., 2015] which compares predictions at the phylum level because "[comparisons at that level are] less prone to differences". The phylum level is four levels removed from genus, let alone species or strain. Our results suggests that the door is now open to metagenome analyses at the highest possible resolution.
While our benchmarks are based on simulated data, our experiments are much more realistic than previous analyses. For example, the Kraken and CLARK papers report results on simulations with ten genomes, whereas we have simulated from 100 genomes and mapped against nearly 2,000. One of the difficulties we faced in our analyses was the technical issue of taxonomic naming and annotation in collating results. This seemingly trivial matter is complicated by the lack of attention paid to low taxonomic level analysis in previous methods. Hopefully our results will spur others to standardize and organize naming conventions and analysis scripts so that low taxonomic level analysis can become routine.
In addition to quantification, we believe there is a lot of potential for differential analysis tools developed for RNA-Seq to be applied more systematically in the metagenomics setting. In that regard, one of the interesting features of kallisto is the ability to bootstrap to assess uncertainty in assignment, and we have utilized this to develop a new method and tool for differential analysis that takes advantage of the feature. The method should immediately be applicable in comparative metagenomics studies.

Acknowledgments
H.P. was supported by an NSF graduate research fellowship. P.M. was partially supported by a Fulbright fellowship. L.S and L.P. were partially supported by NIH R01 HG006129 and NIH R01 DK094699.