Variation in 16S ribosomal gene sequences has been used since the mid-1980's to characterize microbial diversity (Stahl et al., 1984
). Interest in sequence-based surveys of environmental microbes has exploded in recent years with the availability of sequencing technologies that produce ever-larger data sets at ever-decreasing cost; in particular, the Illumina platform is attractive because of throughput, despite its short reads (Sogin et al., 2006
; Lazarevic et al., 2009
; Claesson et al., 2010
; Caporaso et al., 2011
; Degnan and Ochman, 2012
). Here, we examine the reliability of assignment of novel sequences to known taxa under thousands of simulated scenarios, varying primer choice, read length and environment.
Previous work on taxonomic classification of environmental 16S rRNA gene sequences has focused on whether reference sequences matching a given query share taxonomic annotations (Jonasson et al., 2002
; Desantis et al., 2006
; Sogin et al., 2006
; Wang et al., 2007
). Validations of taxonomic classifiers have typically compared a limited range of primers, read lengths and environments (Sundquist et al., 2007
; Huse et al., 2008
; Liu et al., 2008
; Wu et al., 2008
; Hamp et al., 2009
). Reference databases contain many sequences derived from some environments and few associated with others, however (Supplementary Figure S1
), leading to substantial variation in classification quality. In addition, the use of ‘leave-one-out' cross-validation at the sequence level (Sundquist et al., 2007
; Wang et al., 2007
; Liu et al., 2008
; Wu et al., 2008
)—where a single sequence with a known annotation is held out from a reference database and classified using the remainder—is problematic: reference sequences matching held-out query sequences are likely to originate from the same sample, because natural environments contain ‘microdiverse' clusters of closely related strains (Acinas et al., 2004
We addressed these issues by simulating truncated reads from eight large environmental data sets of near-full-length 16S rRNA gene sequences extracted from GreenGenes (Supplementary Table S1
), using pairs of 44 ‘universal' primers commonly found in the literature (Supplementary Tables S2 and S3
). These were selected from an initial set of 94 primers by the criterion that each primer had to match at least 40% of the sequences in at least one of the chosen environmental samples. Single-end reads were tested from each primer with all viable amplification partners (794 combinations), and paired-end reads were tested using all 374 viable pairings of the 22 forward and 22 reverse primers. Simulations using 11 read lengths (32
nt and full-length), with the constraint that read length could not exceed amplicon length for each primer pair, produced 6617 single-end and 3061 paired-end datasets per environment.
Reference databases were constructed by holding out each entire study in turn from GreenGenes, clustering the remainder at 99% using UCLUST (Edgar, 2010
), and selecting one representative sequence per cluster (see Supplementary Methods
for details). Each query fragment was then matched against remaining representative sequences using USEARCH (Edgar, 2010
), configured to penalize indels and mismatches equally. Clusters were then selected that matched within 0.5% identity of the best hit (hits <80% identity were disregarded). For paired-end query sequences, our Rtax procedure (Supplementary Methods
) selected those reference clusters that matched both reads simultaneously with an average percent identity within 0.5% identity of the maximum. Taxonomic classifications were made at each level by retaining annotations agreeing among >50% of the clusters (including those with no annotation in the denominator); these generally extended at best to the genus level, because the reference database provides few species-level annotations.
Sequences from novel taxa (or sequences that appear novel due to sequencing error or chimerism) clearly cannot be correctly classified; however, such sequences may constitute a substantial proportion of a given sample (Supplementary Figure S1
and Supplementary Table S1
). The version of GreenGenes that we used excluded taxa (defined by 97% identity) that were unique to a single sample, as one of the several strategies to remove chimeras. These unique sequences were therefore excluded from our query sets. Thus the classification rates we report represent the proportion of non-chimeric, non-unique sequences that can be classified to each rank. If an environmental sample is not similarly filtered before classification, then the classifiable proportion (that is, taken with respect to the total sample) will be correspondingly lower.
Classification rate and accuracy vary widely among environments and sequence regions, for several reasons: (1) the reference database provides different levels of coverage of each environment, (2) no primer is truly ‘universal' and different primers (and pairs) hit different proportions of sequences in each environment and (3) the targeted regions are variably informative. shows proportions of sequences from each environment classified to each rank, for all 9678 single-ended and paired-end primer and read length combinations. Horizontal panels compare unfiltered results to classifications passing 80% and 95% estimated accuracy filters (see Supplementary Methods
), showing that most classifications can be made with high accuracy when optimal primers are chosen. Remarkably, only 96
nt of sequence (taken as a single read or as a pair of 48
nt reads) can provide 82–100% of the 80% accurate genus classifications available from any read length (Supplementary Table S4
). Paired-end sequencing can provide substantial gains in classification rate for some—but not all—environments and read lengths. Paired-end classifications are typically more accurate than those made from single reads, and so are more likely to pass the 95% estimated accuracy filter (Supplementary Table S5
). Another surprise is that hypervariable regions need not be specifically targeted, as there is no obvious relationship between taxonomic informativeness of a region and the extent to which it overlaps any of the classical ‘V-regions'.
Figure 1 Classification performance, at three levels of estimated accuracy (Supplementary Methods), of 6617 possible choices of amplification primer, sequencing primer and read length for single-ended reads from different environments (left portion of each panel) (more ...)
No one combination of primers and read length works best in all environments, but near-optimal performance in six out of the eight environments is available using paired-end 80
nt reads from primers such as E517F, U515F or E341F paired with E1406R or closely related primers (). However, practical considerations such as ability to amplify low-biomass samples will sometimes influence which primers are used. For instance, short amplicons may be preferred because these are less subject to length heterogeneity biases and chimera formation. Similarly, short single-ended sequences are less subject to errors due to chimeras, simply because they are less likely to span a breakpoint. Classification performance for experimental choices matching such constraints can be found in the supplementary data
Genus classification rates for optimal choices of primers, grouped by total read length
The choice of reference database and taxonomy can have a dramatic impact on the resulting classification accuracy. In this study, we used the current GreenGenes taxonomy, which has been filtered to remove chimeras and where the taxonomic annotations are comprehensive and consistent with the phylogenetic tree (McDonald et al., 2011
). Experiments using a previous version of the GreenGenes taxonomy lacking these features yielded far poorer accuracy (data not shown). In addition, bolstering areas of low coverage in reference databases will substantially improve classifier performance. For instance, taxa in the hypersaline mat, coral and grassland soil samples were underrepresented in the reference database (Supplementary Figure S1
), and—presumably as a consequence—classifications of sequences from those samples were less likely to prove correct (). Additional data sets from poorly sampled environments will also help to distinguish chimeric from legitimate but novel sequences.
In combination, these results indicate that taxonomic classifications of short reads—especially genus-level classifications—should be treated with skepticism, unless the specific combination of primer, read length, environmental source, reference database and assignment method has been thoroughly validated. At the same time, optimal choices of these parameters allow high classification rates and high accuracy. Thus, large-scale projects such as the Earth Microbiome Project (Gilbert et al., 2010
), which aims to collect and analyze samples from tens of thousands of microbial habitats around the globe, may reasonably proceed with standardized primer choices and short reads.