microclass: an Rpackage for 16S taxonomy classification
 Kristian Hovde Liland^{1, 2},
 Hilde Vinje^{1} and
 Lars Snipen^{1}Email authorView ORCID ID profile
DOI: 10.1186/s1285901715832
© The Author(s) 2017
Received: 21 August 2016
Accepted: 3 March 2017
Published: 16 March 2017
Abstract
Background
Taxonomic classification based on the 16S rRNA gene sequence is important for the profiling of microbial communities. In addition to giving the best possible accuracy, it is also important to quantify uncertainties in the classifications.
Results
We present an R package with tools for making such classifications, where the heavy computations are implemented in C++ but operated through the standard R interface. The user may train classifiers based on specialized data sets, but we also supply a readytouse function trained on a comprehensive training data set designed specifically for this purpose. This tool also includes some novel ways to quantify uncertainties in the classifications.
Conclusions
Based on input sequences of varying length and quality, we demonstrate how the output from the classifications can be used to obtain high quality taxonomic assignments from 16S sequences within the R computing environment. The package is publicly available at the Comprehensive R Archive Network.
Keywords
R Taxonomy 16SBackground
The profiling of microbial communities by the sequencing of the 16S rRNA gene has become a standard approach in metagenomics [1]. This means that collected DNA is subject to a targeted sequencing to extract a selected region of the 16S gene from all organisms in the sample. The actual content of the sample can then be described by performing a large scale taxonomic classification of these sequences, i.e. assign them to the proper taxonomic bin, also referred to as binning [2]. Since 16Sbased microbial profiling has become such a widely adopted approach, it is also important that the bioinformatics tools involved are optimized to the highest standard. A widely used tool for this job is the RDPclassifier [3]. It is beyond doubt a good tool for this job, but at the same time it is not perfect, and in a systematic testing of this and other approaches we found there were always other methods that performed better [4]. Alternative tools are a benefit to the scientific community, and here we present a software to be used within the popular R computing environment [5].
There are some issues that must be considered when it comes to making tools for the binning of 16S sequences. First, the pattern recognition algorithm itself must be capable of recognizing the, sometimes small, differences in DNA that separates taxa. It must also handle the huge amount of bins or categories we are facing here, thousands rather than 2–3 which is often the case in textbook literature. Precision also very much depends on the quality of the training data [6]. Due to the ever expanding taxonomy of prokaryotes, there are no comprehensive gold standard data sets available. Along with the microclass package described here, we also supply the microcontax R data package, containing designed data sets based on a consensus taxonomy assignment among several data repositories. This is probably the closest we get to a gold standard today.
Speed is another issue. With today’s sequencing technology and low prices, a data set may easily contain millions of reads. Some procedures for OTU (Operational Taxonomic Unit) picking will start out by classifying reads to predefined taxa [7]. Thus, the number of sequences to classify may be huge. Other procedures cluster the reads before taxonomy assignments, defining OTU’s as ’spherical’ clusters in a space of evolutionary distances approximated by alignment percentage identity, and then classify only the cluster centroids [8]. In some applications, e.g. in forensic applications [9], we are more interested in recognizing specific taxonomic profiles rather than discovering new taxa. In such cases the classification of all reads into predefined bins is clearly what we seek.
Uncertainty is the third issue. In any collection of reads there will be a number of sequences that cannot be given a highconfidence classification. There are several reasons for this. First, the taxonomy itself is not always well defined, and sometimes even highquality sequences fall on the border between existing taxa, making the classification uncertain. Second, due to sequencing errors and chimeras some reads may be difficult to recognize, and third, some microbial communities will contain new taxa not previously seen.
In the presented Rpackage we have implemented some algorithms that have proved efficient and/or are much used for 16S taxonomic classification. Efforts have been made to make them both fast and memoryefficient. All methods can be trained on new data, but we have also supplied the package with a readytouse tool that is already trained and optimized for the contax.trim data set from [10]. This tool also quantifies uncertainties in a new way. The microclass Rpackage, as well as its symbiotic data package microcontax, are freely available at the Comprehensive R Archive Network (CRAN, [11]).
Implementation
The multinomial method
Based on our previous testing of Kmer based classification methods in [4] we found that the best overall results were obtained by the algorithm denoted the multinomial method [12]. Thus, we have focused the attention on this method in this package. The function multinomTrain is used to train a model of this type on any data set containing FASTAformatted sequences along with the correct taxon assignments for each sequence. The function multinomClassify is then used to classify new sequences based on a trained model.
Both training of a multinomial model and classification of new sequences involves counting a large number of Kmers (overlapping words of length K) in the sequences. The overhead when doing such operations is large, and efficient vectorization is difficult to achieve. A direct implementation would also require the computation of a matrix product of size (N×4^{ K }) · (4^{ K }×M), where N and M are the number of sequences to classify and the number of taxa in the training data, respectively. This is a time consuming task for large N, M and K. Therefore, these computations have been implemented in C++ through the Rcpp [13] interface in R, and some shortcuts are made, which will be explained in the following paragraphs.
The nucleotide sequences are first converted to integer vectors my mapping A, C, G and T (or U) to 0, 1, 2, and 3, while all other letters are mapped to  4^{15}. The latter is done to easily discard Kmers including alien symbols when counting. For training of a multinomial model, all Kmers of each taxon are counted. The counting itself is done by sliding a window along the integer vector of each sequence and computing a position as the inner product between [ 4^{ K−1},4^{ K−2},…,4,1] and the integers in the window. For each of the inner products, this position in the taxon’s counting vector is increased by 1. The result is a matrix, X, of size (M×4^{ K }) that holds the counts for all Kmers in all taxa. Finally, each position in the matrix is rescaled to \(\log _{2}\left (\frac {x_{ij}P/4^{K}}{\sum {x_{i\cdot }}P}\right)\), where P is the number of pseudocounts added. This is stored in an (M×4^{ K }) matrix named Q to represent multinomial logprobabilities with pseudocounts.
When classifying new sequences using the multinomial method, we avoid the mentioned matrix product by combining the Kmer counting with summing of multinomial logprobabilities. For each counted Kmer, the corresponding column in the Q matrix is added to the result, thus never explicitly creating the Kmer count matrix or performing the product with Q. As such we reduce from (4^{ K }·M) operations to ((n−K)·M) for a new sequence of length n. For full 16S sequences (with n≈1500 bases) the number of calculations will be lower for K>5.27 and is easily parallelized.
The taxMachine
Users often want a readytouse tool to classify (many) 16S sequences without having to perform all the training. Based on the work in [4] we have arrived at an optimized tool for classifying 16S sequences, called taxMachine in this package. The taxMachine is based on using the multinomial method with a word length of K=8 and a pseudo count of 100. It has been trained on fulllength 16S sequences to recognize full or partial (reads) sequences at the genus level, using the designed and optimized contax.trim data set for training. See the microcontax data package for details. The taxMachine includes computations of classification uncertainties that requires a detailed explanation.
Classification uncertainty
Uncertainty in a taxonomic classification can be split into two types. The first type is when a sequence happens to be very close to the decision boundary between two or more taxa. We can be fairly certain it belongs to one of these taxa, but it lacks the final discriminative power to safely assign it to one of them. The second type of uncertainty occurs when something completely new is seen. This is not uncommon in metagenome samples, and should be flagged separately since it may indicate sequencing errors, chimeras or some novel type of organism.
The dscore
where \(\hat {p}_{l}\) and \(\hat {s}_{l}\) are the predicted mean and standard deviation at sequence length l, using the fitted regression models. Note that p _{ i,2} (and any other posterior logprobability) can be normalized in the same way, using the same fitted regression model.
If we are near a decision boundary we expect d _{ i }≈0 since the second best genus is almost as good as the best. On the other hand, if d _{ i }>>0 it means the predicted genus is much more likely than any other, and we have a high confidence classification.
The rscore
where \(\bar {s}_{g}\) is a smoothed version of s _{ g } as explained below. Thus, the rscore is a standardized measure of how different a sequence is from its predicted genus centre.
When a new sequence is classified, we do not know its true genus. The predicted genus is then used as a plugin in (3), i.e. we use \(\bar {p}_{g}\) and \(\bar {s}_{g}\) where g is the predicted genus. If the resulting r _{ i } has a large negative value, it means the computed \(\tilde {p}_{i,1}\) is much smaller than the average \(\bar {p}_{g}\) for genus g, and sequence i is unlikely to belong to this genus even if this is where it maximizes the posterior probability.
Exactly how negative is the rscore for an unrecognized sequence? To guide this decision we computed the rscores for all sequences in the contax.full data set [10], and from this we computed the empirical cumulative distribution function. For any given r _{ i } value this gives us the probability of having an rscore this small, or smaller, given that the sequence was from the training data. A very small probability means the sequence is very unusual compared to the training data.
Other methods
The package also contains some alternatives to the multinomial method, mostly for comparisons. The RDPclassifier [3] is a popular tool used in many metagenome applications. The version implemented here is a stripped version without the bootstrapping effort to quantify uncertainties in the classifications. It has been implemented in C++ and accelerated similarly to the multinomial method, see above for details.
A classification using BLAST is also included, since this approach has been common. It is both slower and less precise then the other methods. It requires the BLAST+ software to be installed on the system.
Results and discussion
The microclass package provides optimized tools for taxonomic classification of 16S sequence data in the R computing environment. Some well established and proven methods are available to all users of R, with the possibility to train all methods on new and specialized data sets. However, a readymade classification tool, taxMachine, is also supplied as an Rfunction. This has been optimized in several ways to produce the most accurate classifications at the genus level, without consuming too much memory. Specifically, it employs Kmers of length 8, where an increase to K=9 or K=10 comes at high cost in computation time and memory consumption compared to the small gain in accuracy for genus classifications. Pseudo counts have been set to 100 in the taxMachine as a robust compromise regardless of sequence length (see Additional file 1: Figure S1).
The classification of 16S is the most fundamental approach to profiling a microbial community, and due to the explosion in metagenomic research activities, tools for recognizing taxa from 16S sequences (reads) should be tuned to their optimal performance. The taxMachine Rfunction builds on a parallelized sparsematrix implementation of the multinomial method that makes it efficient both with respect to speed and memory usage. It has been trained on the contax.trim data set, containing 38 871 fulllength highquality sequences covering 1774 genera, where all sequences have a consensus taxonomy, making it the closest we get to a wellbalanced gold standard training set.
The proposed implementation of Kmer counting simply discards a word if it contains an ambiguous character. The main reason for this is the added overhead to the computations by introducing another layer of logic to handle these symbols. For instance the occurrence of the letter D in a sequence means that the base in question could be a G, A or T. One could add 1/3 count to each of the three resulting words, but this would require a substantially slower Kmer counting logic. Since informative ambiguous characters (not N) are rarely seen in reads, we chose to disregard these words and keep the speed advantage of the integer logic.
As described in the Implementation section the taxMachine provides information about classification uncertainty, based on the posterior probabilities of the multinomial model. The very first step needed in these computations is to remove the bias from sequence length in the logprobabilities, as suggested in Eq. (1). The right panel of Fig. 1 shows how the normalized posterior logprobabilities have no apparent trends over sequence length, as opposed to the rawvalues in the left panel. This normalization makes it possible to compute uncertainty/reliability scores to sequences regardless of their exact lengths.
The proposed dscore for a sequence is the difference in score between the most likely and the second most likely taxon. A dscore close to 0 means the sequence is close to a decision border, being almost equally similar to both taxa, and more likely to be misclassified. To visualize this, we classified fragments of all sequences in the contax.trim data set using the taxMachine. We considered fragments of typical readlengths; 120–150 and 270–300 bases, which is typical for Illumina HiSeq and MiSeq raw data, and 450–500 bases, which is typical for Roche 454 and merged (pairedend) Illumina MiSeq data. From each of the original 38 871 sequences we sampled 10 such fragments at random locations along each sequence.
Chimera sequences will also result in sequences that are different in Kmer composition from its source sequences. In Additional file 3: Figure S2 we show an example of such a mixture, including d and rscores.
The rscore, and/or its corresponding probability, may be used to discard sequences that appear unusual. As always, the strictness of this procedure will depend on the application. For most applications we would not discard reads unless they are in the lower 1% or 0.1% quantile, at least (probabilities smaller than 10^{−2}−10^{−3}). Instead of fixing some threshold, and discarding reads, one may also use these probabilities as weights, and give reads with small rscores less weight. When tabulating readcounts into a taxonomic profile, this seems like a natural procedure. Conservative estimates of the expected success rates in classifying new reads and full sequences can be found in Additional file 2: Table S1.
The heavy computations of the microclass package are performed in optimized, parallelized C++. This means that the users can comfortably work in R, knowing that reasonably large data can be processed on a personal computer. Larger problems can be tackled in a further parallelized fashion on computational clusters, simply by splitting data into blocks for separate processing. As the package has an open GPL 2/3 licence, reuse of the code in other, possibly pure C++, implementations is allowed as long as the licencing is correct and proper acknowledgements are used.
Conclusions
The package microclass offers tools for taxonomic classification based on 16S rRNA sequence data to the R community. There are function for training classifiers on your own, specialised data sets, and for using these classifiers to classify new sequences. The taxMachine function has synthesized the designed training data from the microcontax data package with the methods of this package, and is our suggested tool for general classifications. It also implements some novel ways to express uncertainties in the classifications, indicating if the input sequences are difficult to recognize.
Abbreviations
 BLAST:

Basic local alignment search tool
 CRAN:

Comprehensive R archive network
 OTU:

Operational taxonomic unit
 RDP:

Ribosomal database project
Declarations
Acknowledgements
Not applicable.
Funding
This paper is a part of the PhDproject of HV, and this project has been 100% financed by the Norwegian University of Life Sciences.
Availability of data and material
The R package is available for free from The Comprehensive R Archive Network [11]. It is most easily obtained by starting R and running install.packages(~microclass~,repos=~http://cran.rproject.org/~) in the console window. All data used in this paper are also publicly available, in the Rpackage microcontax, at [11].
Authors’ contributions
Authors KHL, HV and LS have all contributed significantly to the programming and documentation of the software, and to the preparation and writing of this manuscript. All authors read and approved the final manuscript.
Competing interests
The nonprofit corporation Nofima  Norwegian Institute of Food, Fisheries and Aquaculture Research has no competing interests related to this publication or the presented software. The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Özlem TaştanBishop, (ed).Bioinformatics and Data Analysis in Microbiology. Rhodes University Bioinformatics, Department of Biochemistry, Microbiology and Biotechnology, Rhodes University, South Africa: Caister Academic Press; 2014.Google Scholar
 Schloss PD, Westcott SL. Assessing and improving methods used in operational taxonomic unitbased approaches for 16s rrna gene sequence analysis. Appl Environ Microbiol. 2011; 77(10).
 Wang Q, Garrity GM, Tiedje JM, Cole JR. Naïve bayesian classifier for rapid assignment of rrna sequences into the new bacterial taxonomy. Appl Environ Microbiol. 2007; 73:5261–7.View ArticlePubMedPubMed CentralGoogle Scholar
 Vinje H, Liland KH, Almøy T, Snipen L. Comparing kmer based methods for improved classification of 16s sequences. BMC Bioinformatics. 2015; 16(1):205.View ArticlePubMedPubMed CentralGoogle Scholar
 R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. http://www.Rproject.org/.
 Werner JJ, Koren O, Hugenholtz P, DeSantis TZ, Walters WA, Caporaso JG, Angenent LT, Knight R, Ley RE. Impact of training sets on classification of highthroughput bacterial 16s rrna gene surveys. ISME J. 2012; 6:94–103.View ArticlePubMedGoogle Scholar
 Caporaso J, Kuczynski J, Stombaugh J, Bittinger K, Bushman F, Costello E, Fiere N, Pena A, Goodrich J, Gordon J, Huttley, S GA and Kelley, Knights D, Koenig J, Lozupone C, McDonald D, Muegge B, Pirrung M, Reeder J, Sevinsky J, Turnbaugh P, Walters W, Widmann J, Yatsunenko T, Zaneveld J, Knigh R. Qiime allows analysis of highthroughput community sequencing data. Nat Methods. 2010. doi:10.1038/nmeth.f.303.
 Edgar R. Uparse: highly accurate otu sequences from microbial amplicon reads. Nat Methods. 2013; 10:996–8.View ArticlePubMedGoogle Scholar
 Leake S, Pagni M, Falquet L, Taroni F, Greub G. The salivary microbiome for differentiating individuals: proof of principle. Microbes Infect. 2016; 18:399–405.View ArticlePubMedGoogle Scholar
 Microcontax R Package. https://cran.rproject.org/package=microcontax.
 Comprehensive R Archive Network. https://cran.rproject.org/.
 Liu K, Wong T. Naïve bayesian classifiers with multinomial models for rrna taxonomic assignment. IEEE/ACM Trans Comput Biol Bioinformatics. 2013; 10(5):1334–9.View ArticleGoogle Scholar
 Eddelbuettel D, Francois R. Rcpp: Seamless r and c++ integration. J Stat Softw. 2011; 40(8):1–18.View ArticleGoogle Scholar
 Cleveland WS, Grosse E, Shyu WM. Statistical Models in S In: Chambers JM, Hastie TJ, editors.. Wadsworth & BrooksCole: 1992. p. 8.
 Schirmer M, Ijaz UZ, D’Amore R, Hall N, Sloan WT, Quince. Insight into biases and sequencing errors for amplicon sequencing with the illumina miseq platform. Nucl acid Res. 2015; 43(6).
 Schloss PD, Jenior ML, Koumpouras CC, Westcott SL, Highlander SK. Sequencing 16s rrna gene fragments using the pacbio smrt dna sequencing system. PeerJ. 2016; 4:e1869.View ArticlePubMedPubMed CentralGoogle Scholar