Denoising the Denoisers: an independent evaluation of microbiome sequence error-correction approaches

Sequence acquisition

The HMP mock community (which has even expected abundances) and the ZymoBIOMICS Microbial Community Standard (referred to as the Zymomock community) were sequenced by the Integrated Microbiome Resource at Dalhousie University using an Illumina MiSeq on separate sequencing runs, as previously described using the V4–V5 16S rRNA gene region (Comeau, Douglas & Langille, 2017). Reads were then uploaded to the European Nucleotide Archive (ENA) under accession number PRJEB24409. The Extreme dataset (mock-12) originally presented in the DADA2 paper and the fungal ITS1 dataset (mock-9) were retrieved from the Mockrobiota project (Bokulich et al., 2016). The Extreme dataset was sequenced using an Illumina MiSeq (Callahan et al., 2016) and the fungal mock community was sequenced using an Illumina HiSeq (Bokulich et al., 2016). The soil data set collected from blueberry fields (available under NCBI SRA PRJNA389786) (Yurgel et al., 2017) and the exercise dataset collected from stool of mice that exercised plus controls (ENA accession PRJEB18615) (Lamoureux, Grandy & Langille, 2017) as well as the human-associated dataset of intestinal biopsies of pediatric Crohn’s disease patients plus controls (ENA accession PRJEB21933) (Douglas et al., 2018) were sequenced at the Integrated Microbiome Resource at Dalhousie University.

Filtering

All sample data were filtered using the Microbiome Helper filtering scripts (Comeau, Douglas & Langille, 2017). In summary, primers were trimmed off all reads using Cutadapt (v 1.14) (Martin, 2011) and GNU Parallel (Tange, 2011). Primer-free sequences were then input into the dada2_filter.R script available in Microbiome Helper. This script takes in the maximum expected number of errors allowed as well as a truncation length. The HMP mock community and the Zymomock community were truncated to read lengths of 270 and 210 base pairs for the forward and reverse reads respectively to remove low quality bases at the end of the reads. The shorter length for the reverse reads is a result of their lower overall quality compared to the forward reads. Note that for DADA2, and Deblur to work properly reads need to be of the same length. The single-end reads from the Extreme mock community and the fungal mock community were truncated to 80 base-pairs. The three real datasets: soil, mouse, and human-associated, were truncated to 270 and 210 base-pairs for the forward and reverse reads, respectively. The number of expected errors allowed were defined as three different filtering stringencies: 5 (low), 3 (medium), and 1 (high).

DADA2 pipeline

The DADA2 pipeline was run using Microbiome Helper scripts, which wraps the core algorithms of the DADA2 pipeline (Callahan et al., 2016). Filtered reads were input into the wrapper script dada2_inference.R which runs the DADA2 inference algorithm. Once ASVs are determined, they are passed into DADA2’s chimera-checking algorithm which was run using the wrapper script dada2_chimera_taxa.R to screen out chimeric sequences. The output objects containing ASV sequences and abundances counts were then converted into BIOM table format using convert_dada2_out.R. All DADA2 wrapper scripts were run with default settings.

UNOISE3 pipeline

Filtered reads were input into USEARCH’s (v 10) (Edgar, 2010) fastq_mergepairs command if they were paired-end reads or concatenated together into one FASTQ if they were single-end reads. Next, the merged FASTQ file was converted into a FASTA file using the Microbiome Helper script run_fastq_to_fasta.pl and then used as input for USEARCH’s fastx_uniques command which generated a FASTA containing all the unique sequences and their abundance. Finally, the FASTA containing unique sequences was used as input into USEARCH’s unoise3 (Edgar, 2016) command generating a BIOM table and representative ASVs that were used in subsequent analyses. All USEARCH scripts were run with default settings.

Deblur pipeline

Paired-end filtered reads were stitched together using the Microbiome Helper wrapper script run_pear.pl which wraps the program PEAR (v 0.9.10) (Zhang et al., 2014). This step was skipped for filtered single-end reads. Next, reads were renamed to match a format that was compatible with QIIME2 (Caporaso et al., 2010) and converted into a QIIME2 artifact. Samples were then run through QIIME2’s built-in deblur command using the 16S setting, which uses Greengenes 13_8 (DeSantis et al., 2006) for positive filtering. A separate analysis to test the effect of positive filtering was conducted by using the non-positive filtered output files (“all.biom” and “all.seqs.fa”) from the stand alone version of Deblur. Fungal reads were run using the “other” setting and the UNITE 10.10.2017 database (Kõljalg et al., 2013). Finally, the representative ASV sequences and a BIOM table were exported from the QIIME2 artifact. All other Deblur scripts were run with default settings.

Open reference OTU pipeline

Filtered reads were stitched and imported into a QIIME2 artifact as previously described in the Deblur pipeline. The imported sequences were then dereplicated using the QIIME2 VSEARCH plugin. Dereplicated sequences were filtered for chimeras using the QIIME2 VSEARCH UCHIME reference-based chimera filtering (Edgar et al., 2011) using the Greengenes13_8 97% OTU database for reference. Chimera-checked sequences were then clustered at 97% in an open-reference fashion using the QIIME2 VSEARCH cluster-features-open-reference plugin. The Greengenes 13_8 97% OTU database was used as the reference database during open-reference OTU picking. Finally, singleton OTUs were removed from the OTU table and the table was exported in BIOM format along with its representative sequences in FASTA format.

Run time and memory analysis

The soil dataset was filtered using the low-stringency filter and then individual samples were rarefied to either 5,000, 10,000, 20,000 or 30,000 reads per sample. The different read-depth sets were then run through the three denoising pipelines and user time and maximum memory usage was determined using the GNU time (v 1.7) command.

ASV and OTU analysis of mock communities

ASVs and OTUs were compared against the expected sequences provided with each of the mock communities. This comparison was done using the command-line BLASTN (v 2.7.10) (Altschul et al., 1990) tool against the expected sequences from each community. This allowed us to determine the number of full length 100% matches and 97% matches. All ASVs and OTUs that did not match these criteria were then compared against the SILVA 16S rRNA gene database (v 128) (Pruesse et al., 2007) or the UNITE 10.10.2017 database to find all full length 100% and 97% matches. Any ASVs or OTUs that did not match these databases were then labeled as “Unmatched”. The number of unique expected sequences for each mock community was determined by slicing out the amplified regions using a custom Python (v 3.6.1) script (slice_amplified_region.py, available on GitHub) from the expected sequences of each mock community and then finding the number of unique sequences from this output using USEARCH10’s fastx_uniques command. In addition to analysis done with no abundance filtering, we also compared how filtering of low abundance ASVs/OTUs affected the type and amount of ASVs/OTUs called by each pipeline. This was done by applying a 0.1% minimum abundance filter to each approach and dataset except for the Extreme mock community where a 0.0004% minimum abundance filter was applied. The 0.1% minimum abundance filter was chosen based on the known 0.1% bleed through between Illumina MiSeq runs. A lower minimum abundance filter of 0.0004% was chosen for the Extreme mock community since some of the expected sequences had a lower expected abundance than 0.1%. Histograms with sequence identity to the expected sequences were generated for each pipeline and each mock community except for the fungal mock community due to the community missing some expected sequences from its expected sequence list. All sequences that showed less than 75% identity were binned together.

Abundance data analysis of mock communities

For the HMP, Zymomock and Extreme datasets all ASVs/OTUs that matched at 97% identity or greater with the provided expected sequences based on a BLASTN search were added to the abundance of the corresponding expected taxa. Stacked bar charts of expected taxa relative abundances were created using the ggplot2 (v 2.2.1) (Wickham, 2009) R (v 3.4.3) (R Development Core Team, 2008) package and the cowplot (v 0.9.2) R package (see Data Deposition for exact scripts).

Due to the incomplete nature of the expected sequences for the fungal mock community, UNITE database hits at 97% or greater to an expected species within the community were considered as expected ASVs/OTUs. All other ASVs/OTUs were classified as “Non-Reference” hits.

Analysis of real datasets

Sequences from the three real datasets; soil, human-associated, and exercise were filtered using medium stringencies (allowing up to three expected errors per sequence) for each package and rarified to 5,215, 3,000 and 4,259 reads, respectively. These rarefaction levels were chosen as they were the lowest read count above 2,000, after the DADA2 pipeline was complete. ASV abundance tables output by all approaches were combined into a single table where each biological sample was represented multiple times (once for each pipeline). ASVs/OTUs not called by a specific pipeline were given an abundance of zero in their column (e.g., ASVs only called by Deblur for sampleA were given zero abundances in the columns for DADA2’s and UNOISE3’s outputs of sampleA). Representative sequences were aligned against the Greengenes 13_5 99% OTU aligned sequences and placed on the Greengenes 13_5 99% OTU tree using SEPP (Mirarab, Nguyen & Warnow, 2011). Weighted UniFrac and unweighted UniFrac distances at the ASV/OTU level were then generated using the QIIME1 beta_diversity.py command and principal coordinates were generated using the QIIME1 prinicpal_cordinates.py command.

The Bray–Curtis distance matrix at the genus level was generated by assigning taxonomy to the resulting ASV/OTUs from each package using the RDP classifier (Cole et al., 2014) with the assignTaxonomy function available in the DADA2 package and the rdp_train_set_16 database. Distances were then generated using the summarize_taxa.py and beta_diversity.py commands in QIIME1. Ordination was generated using the metaMDS function in the vegan R package (Oksanen et al., 2018).

Mantel correlations between each distance matrix for each pipeline and each beta-diversity metric (weighted/unweighted UniFrac and Bray–Curtis) were generated using the vegan R package. The number of observed OTUs/ASVs per sample for each real dataset were generated using the rarified combined OTU/ASV tables and the alpha_diversity.py command in QIIME1 using observed OTUs as the metric option.

Differences in methodology and availability between different denoising pipelines

There are important aspects other than accuracy that need to be considered when determining which denoising package a researcher should use for their project. Both DADA2 and UNOISE3 are suggested to be run in a pooled-sample workflow, where all sequences are pooled together during the denoising process (Table 1). This allows them to better account for batch errors across multi-run experiments. Deblur, on the other hand, runs its denoising process sample-by-sample. This approach helps lower Deblur’s computational requirements, but at the cost of reducing its ability to correct multi-run batch effects. Both DADA2 and Deblur are open source projects, whereas UNOISE3 is a closed-source project which has a free 32-bit academic version with a 4 Gb memory cap and a full 64-bit version that costs between $885–1,485 USD (Table 1). Another major difference is that the built-in Deblur function in QIIME2 has a positive filtering process. This default setting causes Deblur to discard reads that do not reach a length-scaled bit score and an e-value threshold to any sequence in the 88% representative sequences Greengenes database. Note the default database can be changed using the “other” version of the Deblur plugin in QIIME2, an important feature when working with fungal or eukaryotic data. It is also important to note that the stand-alone version of Deblur outputs both positively filtered and non-filtered results by default, unlike the QIIME2 plugin which, only outputs the positively filtered results. Currently, the functionality of both DADA2 and Deblur can be accessed through a graphical user interface as plugins in QIIME2 Studio, whereas UNOISE3 does not support a graphical user interface (Table 1).

Total number of ASVs/OTUs vary across approaches in mock communities

We processed four different mock communities with the DADA2, UNOISE3, and Deblur denoising pipelines as well as an open-reference 97% OTU clustering pipeline to compare the resulting ASVs/OTUs from each approach. Focusing on the number of called ASVs from each denoising pipeline, we found no approach consistently called more ASVs. DADA2 called the most ASVs in two mock communities (HMP: 42, Extreme: 78) and UNOISE3 called the most ASVs in the other two mock communities (Fungal: 38, Zymomock: 43) under medium stringency filtering (Fig. 1, Table S1). Overall, open-reference OTU clustering output the most ASVs/OTUs for all four mock communities (HMP: 453, Extreme: 8891, Fungal: 96, Zymomock: 294). None of the approaches were capable at outputting all expected sequences at 100% identity in any of the mock communities that were processed and, in all datasets, at least one denoising pipeline output more ASVs than expected sequences. All four approaches output at least one ASV/OTU at 97% or greater identity from all organisms in the HMP mock community and the Zymomock community (Tables S2, S3). DADA2 output nine more ASVs with 97% or greater identity matches to expected sequences in the Extreme dataset (read depth of 2 million reads) than the other two denoising pipelines (Table S4). Six of the nine taxa that DADA2 called and the other pipelines did not call, had expected relative abundances of only 0.000427% (Table S4). Of the other three taxa not found by Deblur or UNOISE3, one was at an expected abundance of 0.00427% and two had an expected abundance of 0.0427% (Table S4). Interestingly, open reference OTU clustering called OTUs corresponding to these nine expected taxa along with two more expected taxa that DADA2 missed, both of which were in the lowest possible expected abundance range (0.000427%). One organism was missed by all pipelines in the Extreme community (C. methylpentusum), which was also in the lowest possible expected abundance range. Sequence identity histograms were constructed for the HMP, Zymomock and Extreme communities (Figs. S1–S3). In the HMP community, no ASV was found to have below 80% identity by the denoising pipelines, but open reference OTU clustering found 21 OTUs below 80% identity. Furthermore, the OTU pipeline called 324 OTUs below 99% identity whereas all denoising approaches called less than 15 ASVs (Fig. S1). In the Extreme community, most ASVs were found in the 99–100% identity range whereas the majority of OTUs found by open-reference OTU clustering were in the 95–97% identity range (Fig. S2).

Due to the expected sequences for the fungal community being incomplete, we included any UNITE database hits as expected sequences if the matching sequence was from a species that was included in the fungal community. This resulted in almost all the fungi present in the community to be found by all four pipelines except for Penicillium allii and Penicillium commune. Open-reference OTU clustering did find P. allii, although at an extremely low abundance of 0.004% compared to its expected abundance of 6.25%.

Given that some of the above potential spurious ASVs would be removed by sequence bleed-through (Illumina, 2017) or low-abundance filters in typical workflows, we applied an abundance cutoff filter of 0.1% to each mock community except for the Extreme community, where an abundance filter of 0.0004% was applied (Fig. S4). This lower abundance filter was chosen due to some expected sequences being in abundance of 0.000427%. Application of the abundance filter to the HMP community resulted in all 10 unmatched ASVs (those that did not match either the expected or SILVA by 97% or greater) called by DADA2 to be discarded, but none of the four unmatched reads in UNOISE3 to be discarded (Fig. S4). A similar phenomenon was seen in the Zymomock community with all 12 of Deblur’s unmatched reads being discarded (along with one database hit) and UNOISE3 only discarding one of 19 unmatched reads it called. The largest effect that the application of the filters had was the removal of a significant amount of OTUs found in each mock community by the open-reference OTU pipeline. The most drastic change was seen in the Extreme community that started with 8,891 OTUs and was reduced to 1,248 OTUs (Fig. S4). In the fungal community the number of OTUs found by open-reference OTU clustering (27) became less than the number of ASVs found by UNOISE3 (36).

To determine how read quality filtering affects the number of ASVs called by each pipeline, we ran each denoising pipelines using two additional quality filtering stringencies, low and high (see ‘Methods’). The different quality filter stringencies used made only small impacts on the numbers of ASVs called by each pipeline for the HMP, Extreme and fungal datasets (Table S1). A difference of six ASVs was the largest between the high and medium stringencies on the HMP community and was output by UNOISE3 (Table S1). In the Zymomock community, the number of ASVs called by DADA2 only varied by one for all three stringencies, but Deblur varied by as much as 12 ASVs and UNOISE3 varied by as much as 16 ASVs being called between the high and medium filter stringencies (Table S1).

Denoising pipelines are consistent in determining mock community composition

Despite the different ASV counts between each pipeline in the mock community, the relative abundances of the expected taxa are strikingly similar (Fig. 2). In both the HMP and Zymomock communities, only a small portion of ASVs called by DADA2 and Deblur did not match the expected sequences by 97% identity or greater. In contrast, UNOISE3 identified multiple (eight in HMP, 20 in Zymomock) sequences that summed together to make up 2.9% and 4.8% of the relative abundance in the HMP and Zymomock communities, respectively (Tables S1, S2, S3). Open-reference OTU clustering found 6783 non-reference OTUs (Table S1) in the Extreme community that summed together to make up 2.6% of the community, whereas the denoising approaches all found non-reference abundances less than or equal to 0.3% (Table S4). None of the approaches were good at distinguishing the proper abundances of the two Parabacteroides distasonis strains with denoising pipelines finding similar proportions of both strains and the open-reference OTU pipeline finding dominance of the 13,400 strain and not the 13,401 strain (Fig. 2B, Table S4). None of the approaches performed well at matching the expected abundance of the Zymomock community or the fungal community (Fig. 2). All three denoising pipelines called over-abundances of Lactobacillus fermentum in the Zymomock community, with Deblur calling the most (44.5%) and UNOISE3 calling the least (37.6%) (Table S3). Similarly, all approaches called non-reference hits in greater than 9% abundance in the fungal community (Table S5), with UNOISE3 calling the most non-reference hits (12.4%) and Deblur calling the least (9.8%). Due to all four pipelines producing similar proportions of non-reference hits in the fungal community, this could suggest that either the mock compositions are not in the expected proportions, the four pipelines are similarly biased, or that an upstream process during sequencing caused the introduction of unexpected sequences.

Weighted beta-diversity results from different approaches are indistinguishable in real soil and host-associated communities

After comparing each pipeline using mock communities, we next wanted to investigate how comparable the results between pipelines were for real 16S datasets. We compared the pipelines on a soil dataset (soil) due to its high diversity (Fierer & Jackson, 2006) a mouse exercise stool dataset (exercise) and a gut biopsy sample dataset from Crohn’s Disease patients plus controls (human-associated). The intra-sample distances based on weighted UniFrac measurements were comparable between each approach among all three datasets (Fig. 3). In general, we found the intra-samples distances to be small, ranging between median values of ∼0.09 and ∼0.15. There were no consistent differences seen between all three datasets. Looking at the soil dataset, we found that all three denoising pipelines had similarly small intra-sample distances (∼0.11) (Fig. 3A), while the open-reference OTU clustering pipeline was slightly further away from UNOISE3 and Deblur (∼0.13). Despite this observation, it is important to note that this difference is still relatively small. Furthermore, when each sample was plotted onto a PCoA we found that samples tended to group by biological origin rather than the approach used to process the sequences, suggesting that a relatively similar PCoA plot would be generated by each pipeline (Figs. 3D–3F). The Mantel correlations between each weighted UniFrac distance matrix were all highly correlated (correlation values ranging from 0.764 to 0.975), which suggests similar weighted UniFrac profiles between each approach (Tables S6–S8).

Looking at Bray–Curtis dissimilarity, another metric that takes in abundance information, we found that intra-sample distances based on genus-level assignment were similar among different approaches (Fig. S5). In general, the intra-sample distances were relatively low, with median values ranging from ∼0.04 to ∼0.20, indicating a high amount of agreement between approaches. Interestingly, we found that the intra-sample distances seemed to increase across datasets based on how diverse the samples within them were (Fig. S5). In two of the three datasets (soil and BISCUIT), we found that DADA2 and UNOISE3 tended to be closer together in distance than Deblur was to DADA2 or UNOISE3, indicating a slightly higher amount of agreement between DADA2 and UNOISE3 at the genus level. Again, these differences were relatively small and would have minimal impacts on biological results obtained from them. Plotting the Bray-Curtis dissimilarity matrices onto an NMDS plot resulted in similar findings as weighted UniFrac, with samples grouping by biological origin rather than the pipeline used to process them. Furthermore, we found that the Bray-Curtis dissimilarity matrices were extremely well correlated with each other, ranging in values between 0.956-0.995 (Figs. S6–S8).

We next compared unweighted UniFrac distances, which is a metric that considers the presence or absence of ASVs/OTUs and their phylogenetic distance. We found that the median intra-sample distance between each pipeline was much greater, ranging between ∼0.40 to ∼0.79 (Fig. S6). Similar to our results focused on Bray-Curtis dissimilarity, we also found that the median intra-sample distance between samples tended to increase with sample diversity (0.72–0.79 for soil; 0.40–0.55 for human-associated; 0.60–0.70 for exercise). The PCoA plot based on these distances resulted in samples grouping together based on pipeline, rather than biological origin, indicating large differences between the different approaches (Fig. S6).

To look at how much of a difference the positive filtering process had on the profiles generated by Deblur and to see if the slight difference in Bray-Curtis dissimilarity among the denoising approaches was caused by this filter, we ran the stand-alone version of Deblur and examined the non-positive filtered results. We found that the intra-sample distances based on any of the three metrics tested (weighted/unweighted UniFrac and Bray-Curtis dissimilarity) between any of the pipelines did not vary (Fig. S7).

Alpha-diversity metrics vary between denoising pipelines

We next investigated how the number of ASVs/OTUs called in a real dataset differed between processing pipelines. In the soil dataset, DADA2 called 16609 ASVs, UNOISE3 called 11613 ASVs, Deblur called 8270 ASVs and open-reference OTU clustering called 21297 OTUs (Fig. 4A). Across all datasets, we found that DADA2 called the most ASVs among the denoising pipelines and that open-reference OTU clustering called the most ASVs/OTUs overall (Figs. 4A, 4C, 4E). On average, DADA2 called 727 more ASVs than Deblur and 532 more ASVs than UNOISE3 while open-reference OTU clustering called 3135 more OTUs/ASVs than DADA2. In all datasets, Deblur called the least amount of ASVs/OTUs. Looking at the number of OTUs called per sample, we found that DADA2 did not correlate well with UNOISE3 or open-reference OTU clustering and in the soil dataset the correlation between DADA2 and UNOISE3 was found to not be significant (p = 0.054) (Figs. 4B, 4D, 4F). This is a concerning result as it indicated the possibility of different biological results based on the pipeline that was chosen to process the data.

One interesting result was that, despite DADA2 calling the most total ASVs among the denoising pipelines, it called the least amount per sample (Figs. S8–S10) and, among denoising pipelines, UNOISE3 tended to call the most (Figs. S8–S10). However, open-reference OTU clustering called the most OTUs/ASVs overall per sample among all approaches that were tested with the exception of the soil dataset.

Computational requirements are vastly different across denoising pipelines

Knowing that all three of these pipelines resulted in similar relative abundance profiles on mock communities and small intra-sample distances on real 16S communities, we next investigated how the run time and memory usage differed between the denoising approaches. We found that UNOISE3 (4.6 min) was 1,273 times faster than DADA2 (5,834 min) and 15 times faster than Deblur (69.3 min) at a total read count of 1,926,000 reads evenly distributed across 103 samples (Fig. 5A). Run times for all pipelines increased as the number of reads per sample increased. Deblur used a static amount of memory (611 Mb) as reads per sample increased, whereas in general the other two pipelines increased in memory usage as the number of reads per sample increased, with the exception of DADA2 run at 1,926,000 reads (Fig. 5B). Deblur used the smallest amount of memory at the maximum read count of 1,926,000 reads. We found that DADA2 had the highest amount of memory usage (4,071 Mb at 1,287,000 reads) among the three pipelines. Interestingly, this usage was more than the amount used at the maximum read count (3,600 Mb). In addition, none of the runs exceeded the 4 Gb memory cap on the 32-bit free academic version of USEARCH10.