Sequencing 16S rRNA gene fragments using the PacBio SMRT DNA sequencing system

Community DNA

We utilized genomic DNA isolated from four communities. These same DNA extracts were previously used to develop an Illumina MiSeq-based sequencing strategy (Kozich et al., 2013). Briefly, we used a “Mock Community” composed of genomic DNA from 21 bacterial strains: Acinetobacter baumannii ATCC 17978, Actinomyces odontolyticus ATCC 17982, Bacillus cereus ATCC 10987, Bacteroides vulgatus ATCC 8482, Clostridium beijerinckii ATCC 51743, Deinococcus radiodurans ATCC 13939, Enterococcus faecalis ATCC 47077, Escherichia coli ATCC 70096, Helicobacter pylori ATCC 700392, Lactobacillus gasseri ATCC 33323, Listeria monocytogenes ATCC BAA-679, Neisseria meningitidis ATCC BAA-335, Porphyromonas gingivalis ATCC 33277, Propionibacterium acnes DSM 16379, Pseudomonas aeruginosa ATCC 47085, Rhodobacter sphaeroides ATCC 17023, Staphylococcus aureus ATCC BAA-1718, Staphylococcus epidermidis ATCC 12228, Streptococcus agalactiae ATCC BAA-611, Streptococcus mutans ATCC 700610, Streptococcus pneumoniae ATCC BAA-334. The mock community DNA is available through BEI resources (v3.1, HM-278D). Genomic DNAs from the three other communities were obtained using the MO BIO PowerSoil DNA extraction kit. The human and mouse fecal samples were obtained using protocols that were reviewed and approved by the University Committee on Use and Care of Animals (Protocol #PRO00004877) and the Institutional Review Board at the University of Michigan (Protocol #HUM00057066). The human stool donor provided informed consent.

Library generation and sequencing

The DNAs were each amplified in triplicate using barcoded primers targeting the V4, V1–V3, V3–V5, V1–V5, V1–V6, and V1–V9 variable regions (Table 1). The primers were synthesized so that the 5′ end of the forward and reverse primers were each tagged with paired 16-nt symmetric barcodes (https://github.com/PacificBiosciences/Bioinformatics-Training/wiki/Barcoding-with-SMRT-Analysis-2.3) to allow multiplexing of samples within a single sequencing run. Methods describing PCR, amplicon cleanup, and pooling were described previously (Kozich et al., 2013). The SMRTbell adapters were ligated onto the PCR products and the libraries were sequenced by Pacific Biosciences using the P6-C4 chemistry on a PacBio RS II SMRT DNA Sequencing System. Diffusion Loading was used for regions V4, V1–V3, and V3–V5 and MagBead loading was used for regions V1–V5, V1–V6, and V1–V9. Each region was sequenced separately using movies ranging in length between 180 and 360 min. The sequences were processed using pbccs (v.3.0.1; https://github.com/PacificBiosciences/pbccs), which generates predicted error rates using a proprietary algorithm.

Data analysis

All sequencing data were curated using mothur (v1.36) (Schloss et al., 2009) and analyzed using the R programming language (R Core Team, 2016). The raw data can be obtained from the Sequence Read Archive at NCBI under accession SRP051686, which are associated with BioProject PRJNA271568. This accession and bioproject also contain data from the same samples sequenced using the P4-C2 chemistry. Several specific features were incorporated into mothur to facilitate the analysis of PacBio sequence data. First, because non-ambiguous base calls are assigned to Phred quality scores of zero, the consensus fastq files were parsed so that scores of zero were interpreted as corresponding to an ambiguous base call (i.e., N) in the fastq.info command using the pacbio=T option. Second, because the consensus sequence can be generated in the forward and reverse complement orientations, a checkorient option was added to the trim.seqs command in order to identify the proper orientation. These features were incorporated into mothur v.1.30. Because chimeric molecules can be generated during PCR and would artificially inflate the sequencing error, it was necessary to remove these data prior to assessing the error rate. Because we knew the true sequences for the strains in the mock community we could calculate all possible chimeras between strains in the mock community (in silico chimeras). If a sequence read was 3 or more nucleotides more similar to an in silico chimera than it was to a non-chimeric reference sequence, it was classified as a chimera and removed from further consideration. Identification of in silico chimeras and calculation of sequencing error rates was performed using the seq.error command in mothur (Schloss, Gevers & Westcott, 2011). De novo chimera detection was also performed on the mock and other sequence data using the abundance-based algorithm implemented in UCHIME (Edgar et al., 2011). Sequences sequences were aligned against a SILVA-based reference alignment (Pruesse et al., 2007) using a profile-based aligner (Schloss, 2009) and were classified against the SILVA (v123) (Pruesse et al., 2007), RDP (v10) (Cole et al., 2013), and greengenes (v13_8_99) (Werner et al., 2011) reference taxonomies using a negative Bayesian classifier implemented within mothur (Wang et al., 2007). Sequences were assigned to operational taxonomic units using the average neighbor clustering algorithm with a 3% distance threshold ( Schloss & Westcott, 2011 ). Detailed methods including this paper as an R markdown file are available as a public online repository (http://github.com/SchlossLab/Schloss_PacBio16S_PeerJ_2016).

The PacBio error profile and a basic sequence curation procedure

To build a sequence curation pipeline, we first needed to characterize the error rate associated with sequencing the 16S rRNA gene. Using the consensus sequence obtained from at least 3 sequencing passes of each fragment, we observed an average sequencing error rate of 0.65%. Insertions, deletions, and substitutions accounted for 31.2, 17.9, and 50.9% of the errors, respectively. The substitution errors were equally likely and all four bases were equally likely to cause insertion errors. Interestingly, guanines (39.4%) and adenines (24.3%) were more likely to be deleted than cytosines (18.3%) or thymidines (18.0%). The PacBio quality values varied between 2 and 93. Surprisingly, the percentage of base calls that had the maximum quality value did not vary among correct base calls (80.5%), substitutions (80.0%), and insertions (80.4%). It was not possible to use the individual base quality scores to screen sequence quality as has been possible in past studies using the Phred-based quality scores that accompany data generated using the 454 and Illumina technologies. We did observe a strong correlation between our observed error rate and the predicted error rate as calculated by the PacBio software (Pearson’s R: −0.67; Fig. 1A).

We established a simple curation procedure by culling any sequence that had a string of the same base repeated more than eight times or did not start and end at the expected alignment coordinates for that region of the 16S rRNA gene. This reduced the experiment-wide error rate from 0.68 to 0.65%. This basic procedure resulted in the removal of between 0.714 (V1–V3) and 9.47 (V1–V9)% of the reads (Table 2). Although the percentage of reads removed increased with the length of the fragment, there was no obvious relationship between fragment length and error rate (Fig. 2).

Identifying correlates of increased sequencing error

In contrast to the 454 and Illumina-based platforms where the sequencing quality decays with length, the consensus sequencing approach employed by the PacBio sequencer is thought to generate a uniform distribution of errors. This makes it impossible to simply trim sequences to high-quality regions. Therefore, we sought to identify characteristics within sequences that would allow us to identify and remove those sequences with errors using three different approaches. First, we hypothesized that errors in the barcode and primer would be correlated with the error rate for the entire sequence. We observed a strong relationship between the number of mismatches to the barcodes and primers and the error rate of the rest of the sequence fragment (Fig. 1B). Although allowing no mismatches to the barcodes and primers yielded the lowest error rate, that stringent criterion removed a large fraction of the reads from the dataset. Allowing at most one mismatch only marginally increased the error rate while retaining more sequences in the dataset (Fig. 2). Second, we hypothesized that increased sequencing coverage should yield lower error rates. We found that once we had obtained 10-fold coverage of the fragments, the error rate did not change appreciably (Fig. 1C). When we compared the error rates of reads with at least 10-fold coverage to those with less coverage, we reduced the error rate by 8.48–37.08% (Fig. 2). Third, based on the observed correlation between the predicted and observed error rates, we sought to identify a minimum predicted error rate that would allow us to reduce the observed error rate. The average observed error rate for sequences with predicted error rates between 0.01 and 0.10% was linear. We decided to use a threshold of 0.01% because a large number of sequence reads were lost when we used a smaller threshold. When we used this threshold, we were able to reduce the error rate by 51.4–70.0% (Fig. 2). Finally, we quantified the effect of combining filters. We found that any combination of filters that included the predicted error rate threshold had the most significant impact on reducing the observed error rate. Furthermore, the inclusion of the mismatch and coverage filters had a negligible impact on error rates, but had a significant impact on the number of sequences included in the analysis. For instance, among the V1–V9 data, requiring sequences to have a predicted error rate less than 0.01% resulted in a 69.2% reduction in error and resulted in the removal of 53.5% of the sequences. Adding the mismatch or coverage filter had no effect on the reduction of error, but resulted in the removal of 56.2 and 56.6% of the sequences, respectively. Use of all filters had no impact on the reduction in the observed error rate, but resulted in the removal of 59.1% of the sequences. The remainder of this paper only uses sequences with a predicted error rate less than 0.01%.

Pre-clustering sequences to further reduce sequencing noise

Previously, we implemented a pre-clustering algorithm where sequences were sorted by their abundance in decreasing order and rare sequences are clustered with a more abundant sequence if the rare sequences have fewer mismatches than a defined threshold when compared to the more abundant sequence (Huse et al., 2010; Schloss, Gevers & Westcott, 2011). The recommended threshold was a 1-nt difference per 100-nt of sequence data. For example, the threshold for 250 bp fragment from the V4 region would be 2 nt or 14 for the 1,458 bp V1–V9 fragments. This approach removes residual PCR and sequencing errors while not overwhelming the resolution needed to identify OTUs that are based on a 3% distance threshold. The tradeoff of this approach is that one would be unable to differentiate V1–V9 sequences that truly differed by less than 14 nt. When we applied this approach to our PacBio data, we observed a reduction in the error rate between 33.0 (V1–V3) and 48.7% (V1–V9). The final error rates varied between 0.02 (V1–V5) and 0.2% (V4). The full-length (i.e., V1–V9) fragments had an error rate of 0.03% (Fig. 2; Table 2), this is similar to what we have previously observed using the 454 and Illumina MiSeq platforms (0.02%) (Schloss, Gevers & Westcott, 2011; Kozich et al., 2013).

Effects of error rates on OTU assignments

The sequencing error rate is known to affect the number of OTUs that are observed (Schloss, Gevers & Westcott, 2011). For each region, we determined that if there were no chimeras or PCR or sequencing errors, then we would expect to find 19 OTUs. When we achieved perfect chimera removal, but allowed for PCR and sequencing errors, we observed between 0.5 (V1–V9) and 14.4 (V4) extra OTUs (Table 2). The range in the number of extra OTUs was largely explained by the sequencing error rate (Pearson’s R = 1.0). Next, we determined the number of OTUs that were observed when we used UCHIME to identify chimeric sequences. Under these more realistic conditions, we observed between 8.2 (V1–V5) and 29.9 (V4) extra OTUs. Finally, we calculated the number of OTUs in the soil, mouse, and human samples using the same pipeline with chimera detection and removal based on the UCHIME algorithm. Surprisingly, there was not a clear relationship across sample type and region. Again, we found that there was a strong correlation between the number of observed OTUs and the error rate for the mouse (R = 0.95) and human samples (R = 0.60). These results underscore the effect of sequencing error on the inflation of the number of observed OTUs.

Classification varies by region, environment, and database

We classified all of the sequence data we generated using the naïve Bayesian classifier using the RDP, SILVA, and greengenes reference taxonomies (Fig. 3). In general, increasing the length of the region improved the ability to assign the sequence to a genus or species. Interestingly, each of the samples we analyzed varied in the ability to assign its sequences to the depth of genus or species. Furthermore, the reference database that did the best job of classifying the sequences varied by sample type. For example, the SILVA reference did the best for the human feces and soil samples and the RDP did the best for the mouse feces samples. An advantage of the greengenes database is that it contains information for 2,514 species-level lineages for 11% of the reference sequences; the other databases only provided taxonomic data to the genus level. There was a modest association between the length of the fragment and the ability to classify sequences to the species-level for the human samples; there was no such association for the mouse and soil samples. In fact, at most 6.2% of the soil sequences and 4.3% of the mouse sequences could be classified to a species. These results indicate that the ability to classify sequences to the genus or species level is a function of read length, sample type, and the reference database.

Sequencing errors are not random

Above, we described that although there was no obvious bias in the substitution or insertion rate, we did observe that guanines and adenines were more likely to be deleted than cytosines or thymidines. This lack of randomness in the error profile suggested that there might be a systematic non-random distribution of the errors across the sequences. This would manifest itself by the creation of duplicate sequences with the same error. We identified all of the mock community sequences that had a 1-nt difference to the true sequence (Fig. 4). For these three regions, between 70.7 and 88.9 of the sequences with 1-nt errors were only observed once. We found that the frequency of the most abundant 1-nt error paralleled the number of sequences. Surprisingly, the same 1-nt error appeared 1,954 times (0.02%) in the V1–V6 mock data and another 1-nt error appeared 1,070 times (0.03%) in the V1–V9 mock data. Contrary to previous reports (Carneiro et al., 2012; Koren et al., 2012), these results indicate that reproducible errors occur with the PacBio sequencing platform and that they can be quite abundant. Through the use of the pre-clustering step described above these 1-nt errors would be ameliorated; however, this result indicates that caution should be used when attempting to use fine-scale OTU definitions.