OTUs vs ASVs: Denoising with DADA2

Turning raw amplicon reads into a feature table begins with one decision: what counts as a feature. Historically that unit was the operational taxonomic unit, or OTU; modern pipelines instead resolve amplicon sequence variants, or ASVs. The choice propagates into every diversity metric and taxonomy assignment you compute downstream.

OTUs versus ASVs

Classic clustering groups reads that are at least 97 percent identical into OTUs, which lumps genuine biological variation together with sequencing error inside fuzzy bins whose membership shifts with the reference database and the dataset. DADA2 instead learns a per-run error model and uses it to infer exact, error-corrected sequences that can differ by a single nucleotide, so an ASV is a precise, reproducible unit that is directly comparable across studies without re-clustering. That resolution and portability are why ASVs have largely replaced OTUs, and the workflow that produces them chains filterAndTrim, learnErrors, dada, mergePairs, makeSequenceTable, and removeBimeraDenovo.

The DADA2 workflow

library(dada2)

path <- "fastq"
fnFs <- sort(list.files(path, pattern = "_R1_001.fastq.gz", full.names = TRUE))
fnRs <- sort(list.files(path, pattern = "_R2_001.fastq.gz", full.names = TRUE))
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))

out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs,
                     truncLen = c(240, 160),
                     maxN = 0, maxEE = c(2, 2), truncQ = 2,
                     rm.phix = TRUE, compress = TRUE, multithread = TRUE)
head(out)

                          reads.in reads.out
sample1_R1_001.fastq.gz       7793      7113
sample2_R1_001.fastq.gz       5869      5299
sample3_R1_001.fastq.gz       5958      5463
sample4_R1_001.fastq.gz       3183      2914

errF <- learnErrors(filtFs, multithread = TRUE)
errR <- learnErrors(filtRs, multithread = TRUE)

dadaFs <- dada(filtFs, err = errF, multithread = TRUE)
dadaRs <- dada(filtRs, err = errR, multithread = TRUE)

mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose = TRUE)

seqtab <- makeSequenceTable(mergers)
dim(seqtab)

seqtab.nochim <- removeBimeraDenovo(seqtab, method = "consensus",
                                    multithread = TRUE, verbose = TRUE)
sum(seqtab.nochim) / sum(seqtab)

[1]  20 293
Identified 61 bimeras out of 293 input sequences.
[1] 0.9643574

The same run in QIIME 2

qiime dada2 denoise-paired \
  --i-demultiplexed-seqs demux.qza \
  --p-trunc-len-f 240 \
  --p-trunc-len-r 160 \
  --p-trim-left-f 0 \
  --p-trim-left-r 0 \
  --p-n-threads 0 \
  --o-representative-sequences rep-seqs.qza \
  --o-table table.qza \
  --o-denoising-stats denoising-stats.qza

Saved FeatureTable[Frequency] to: table.qza
Saved FeatureData[Sequence] to: rep-seqs.qza
Saved SampleData[DADA2Stats] to: denoising-stats.qza