Chapter 2: Amplicon Read Processing, Error Modeling & ASVs with DADA2
Chapter 2: Amplicon Read Processing, Error Modeling & ASVs with DADA2
Historically, amplicon reads were clustered into Operational Taxonomic Units (OTUs) at an arbitrary 97% sequence identity threshold. Today, OTUs have been superseded by Amplicon Sequence Variants (ASVs), which resolve biological variation down to exact single-nucleotide differences.
The DADA2 Paradigm
DADA2 (Divisive Amplicon Denoising Algorithm 2) models Illumina sequencing error rates directly from the data. It calculates the statistical probability that an observed low-abundance sequence is an erroneous artifact derived from an abundant true template versus a genuine novel biological strain.
1. Complete DADA2 Workflow in R
library(dada2)
# File path setup
fnFs <- sort(list.files("raw_fastqs/", pattern="_R1.fastq.gz", full.names=TRUE))
fnRs <- sort(list.files("raw_fastqs/", pattern="_R2.fastq.gz", full.names=TRUE))
sample_names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)
# 1. Filter and Trim (remove primers and low quality tails)
filtFs <- file.path("filtered", paste0(sample_names, "_F_filt.fastq.gz"))
filtRs <- file.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=8)
# 2. Learn Error Rates from the sequencing data
errF <- learnErrors(filtFs, multithread=8)
errR <- learnErrors(filtRs, multithread=8)
# 3. Sample Dereplication and Denoising
dadaFs <- dada(filtFs, err=errF, multithread=8)
dadaRs <- dada(filtRs, err=errR, multithread=8)
# 4. Merge Paired Reads (require minimum 12 bp overlap)
mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose=TRUE)
# 5. Construct Sequence Table (ASV count matrix)
seqtab <- makeSequenceTable(mergers)
# 6. Remove PCR Chimeras using de novo consensus
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=8, verbose=TRUE)
print(paste("Fraction non-chimeric reads:", sum(seqtab.nochim)/sum(seqtab)))