Selection of primers for optimal taxonomic classification of environmental 16S rRNA gene sequences.

Microbial community profiling using 16S rRNA gene sequences requires accurate taxonomy assignments. 'Universal' primers target conserved sequences and amplify sequences from many taxa, but they provide variable coverage of different environments, and regions of the rRNA gene differ in taxonomic informativeness--especially when high-throughput short-read sequencing technologies (for example, 454 and Illumina) are used. We introduce a new evaluation procedure that provides an improved measure of expected taxonomic precision when classifying environmental sequence reads from a given primer. Applying this measure to thousands of combinations of primers and read lengths, simulating single-ended and paired-end sequencing, reveals that these choices greatly affect taxonomic informativeness. The most informative sequence region may differ by environment, partly due to variable coverage of different environments in reference databases. Using our Rtax method of classifying paired-end reads, we found that paired-end sequencing provides substantial benefit in some environments including human gut, but not in others. Optimal primer choice for short reads totaling 96 nt provides 82-100% of the confident genus classifications available from longer reads.

Variation in 16S ribosomal gene sequences has been used since the mid-1980's to characterize microbial diversity (Stahl ). 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 ; Lazarevic ; Claesson ; Caporaso ; 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 ; Desantis ; Sogin ; Wang ). Validations of taxonomic classifiers have typically compared a limited range of primers, read lengths and environments (Sundquist ; Huse ; Liu ; Wu ; Hamp ). 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 ; Wang ; Liu ; Wu )—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 ).

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, 48 nt, 64 nt, 80 nt, 96 nt, 112 nt, 128 nt, 260 nt, 400 nt, 800 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. Figure 1 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) and 3061 possible choices of primer pair and read length for paired-end reads (right portion). Combinations of primers and read lengths are sorted on the x axis according to a measure of overall classification performance (Supplementary Methods). Stacked bars show the proportion of non-chimeric, non-unique sequences from each sample—not the proportion of the total sample—that can be classified to each taxonomic level for each combination. See Supplementary Figure S1 and Supplementary Table S1 for the excluded proportion of novel (and thus a priori unclassifiable) sequences in each sample. The top of each colored section indicates how much of the sample can be classified to the given level or better. ‘Primer miss' (black) indicates sequences that did not match a given primer and so would not be amplified. Classifications more specific than the genus level are exceedingly rare and so are not visible here. Horizontal lines indicate the maximum proportion of each sample classifiable to the genus level using 96 nt or less of sequence (i.e., with an optimal choice of primer or primer pair; see also Supplementary Tables S4 and S5), showing that short reads from the best primers frequently—but not always—provide taxonomic information nearly matching that obtained from longer read lengths. Full-size versions of these panels are available in the supplementary data.

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 (Table 1). 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.

Table 1

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 ). 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 (Figure 1). 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 ), 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.