DADA2 ASV Pipeline — Skrip & Penjelasan Langsung

Dokumen ini menyajikan skrip R lengkap untuk pipeline DADA2 sekaligus penjelasan langsung di dalam MDX.
Dengan format ini, Anda bisa mempelajari setiap blok kode beserta konteks dan alasannya.

0) Setup

library(dada2)       # core pipeline
library(ShortRead)   # optional: peek fastq
library(Biostrings)  # sequences

Paket yang dipanggil:

dada2: inti pipeline ASV inference.
ShortRead: untuk membaca/mengecek isi FASTQ (opsional).
Biostrings: memanipulasi sekuens DNA.

projdir <- "./../.."
fastqdir <- file.path(projdir, "data")
outdir   <- file.path(projdir, "output/dada2_outputs")
dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
refdir <- file.path(projdir, "reference")

projdir: akar proyek (diset naik dua folder).
fastqdir: tempat FASTQ hasil Cutadapt (sudah bebas adaptor/primer).
outdir: folder output DADA2, dibuat otomatis bila belum ada.
refdir: tempat menyimpan file referensi (SILVA).

1) Filtering & Trimming

fnFs <- sort(list.files(fastqdir, pattern="_R1_001\.trim\.fastq\.gz$", full.names=TRUE))
fnRs <- sort(list.files(fastqdir, pattern="_R2_001\.trim\.fastq\.gz$", full.names=TRUE))
stopifnot(length(fnFs) == length(fnRs), length(fnFs) > 0)
samples <- gsub("_S.*_L001_R1_001\.trim\.fastq\.gz$", "", basename(fnFs))
filtFs <- file.path(outdir, paste0(samples, "_F_filt.fastq.gz"))
filtRs <- file.path(outdir, paste0(samples, "_R_filt.fastq.gz"))

fnFs & fnRs: daftar file FASTQ hasil trimming (paired-end).
samples: nama sampel bersih dengan regex (ubah sesuai pola file Anda).
filtFs, filtRs: nama file output hasil filter.

truncLen <- c(270, 220)
maxEE    <- c(2, 2)
truncQ   <- 2
minLen   <- 100
 
set.seed(1)
out <- filterAndTrim(fnFs, filtFs,
                     fnRs, filtRs,
                     truncLen=truncLen,
                     maxEE=maxEE,
                     truncQ=truncQ,
                     minLen=minLen,
                     rm.phix=TRUE,
                     compress=TRUE,
                     multithread=TRUE)

Parameter penting:

truncLen: panjang yang dipertahankan (disesuaikan dengan kualitas dari FastQC). Contoh: R1=270, R2=220 → overlap ±30 bp untuk V3–V4 (amplicon ~460 bp).
maxEE: batas error per read. c(2,2) artinya ≤2 error pada masing-masing R1/R2.
truncQ: potong bila kualitas dibawah 2.
minLen: buang reads dibawah 100 nt.

Output out: ringkasan jumlah reads masuk & lolos filter.

2) Error Model

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

DADA2 membuat profil error dari data untuk mengoreksi kesalahan sekuensing.
Cek dengan plotErrors(errF) untuk validasi visual.

3) Dereplikasi & Inferensi DADA

derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)
names(derepFs) <- samples
names(derepRs) <- samples
 
dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)

Dereplikasi: menggabungkan reads identik → hemat memori/CPU.
DADA inference: menghasilkan ASV dengan resolusi nukleotida tunggal.

4) Merge Pairs

mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)

Menggabungkan pasangan F/R menjadi sekuens penuh.
Butuh overlap ≥20–30 nt. Jika gagal merge, periksa truncLen atau kualitas ujung reads.

5) Tabel ASV & Chimera Removal

seqtab <- makeSequenceTable(mergers)
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
 
getN <- function(x) sum(getUniques(x))
track <- cbind(data.frame(input=rowSums(out)),
               filtered=out[,2],
               denoisedF=sapply(dadaFs, getN),
               denoisedR=sapply(dadaRs, getN),
               merged=sapply(mergers, getN),
               nonchim=rowSums(seqtab.nochim))
rownames(track) <- samples
saveRDS(seqtab.nochim, file.path(outdir, "seqtab_nochim.rds"))

makeSequenceTable: membangun matriks sampel × ASV.
removeBimeraDenovo: menghapus chimera (artefak PCR).
track: pelacakan jumlah reads dari input hingga nonchim.

6) Assign Taxonomy

silva_train <- file.path(refdir, "silva_nr99_v138.1_train_set.fa.gz") 
silva_species <- file.path(refdir, "silva_species_assignment_v138.1.fa.gz") 
 
taxa <- assignTaxonomy(seqtab.nochim, silva_train, multithread=TRUE, tryRC=TRUE)
taxa <- addSpecies(taxa, silva_species, verbose=TRUE)
saveRDS(taxa, file.path(outdir, "taxa_silva_v1381.rds"))

Database SILVA v138.1 digunakan untuk memberi label taksonomi.
assignTaxonomy: sampai level genus.
addSpecies: coba menambahkan level spesies (sering kosong di 16S).

7) END-OF-RUN: Konsolidasi Output

Semua keluaran disatukan dalam blok 99 agar mudah dicari.

# 99.1 Per-sample summary
per_sample <- data.frame(
  Sample        = rownames(track),
  Input         = track[, "input"],
  Filtered      = track[, "filtered"],
  DenoisedF     = track[, "denoisedF"],
  DenoisedR     = track[, "denoisedR"],
  Merged        = track[, "merged"],
  Nonchim       = track[, "nonchim"],
  RetentionPct  = round(100 * track[, "nonchim"] / pmax(track[, "input"], 1), 2),
  ChimeraFracPct= round(100 * pmax(track[, "merged"] - track[, "nonchim"], 0) / pmax(track[, "merged"], 1), 2),
  ASVsPerSample = rowSums(seqtab.nochim > 0),
  stringsAsFactors = FALSE
)

RetentionPct: persentase reads yang bertahan sampai akhir.
ChimeraFracPct: fraksi chimera per sampel.

# 99.2 Global stats
asv_ids     <- paste0("ASV", seq_len(ncol(seqtab.nochim)))
asv_seq_chr <- colnames(seqtab.nochim)
asv_len_vec <- width(DNAStringSet(asv_seq_chr))
chim_frac_global_pct <- round(100 * (1 - (sum(seqtab.nochim) / sum(seqtab))), 2)

Menghitung jumlah ASV, distribusi panjang, fraksi chimera global.

# 99.3 ASV length distribution
asv_len_dist <- as.data.frame(table(ASV_Length_bp = as.integer(asv_len_vec)))

Distribusi panjang sekuens ASV.

# 99.4 Taxonomy tables
seq_to_asv <- setNames(asv_ids, asv_seq_chr)
tax_df <- as.data.frame(unclass(taxa))
tax_df$ASV <- unname(seq_to_asv[rownames(tax_df)])
tax_df <- tax_df[, c("ASV", setdiff(colnames(tax_df), "ASV"))]

Membuat tabel taksonomi dengan ID ASV ringkas.
Versi dengan bootstrap confidence juga dibuat.

# 99.5 Abundance & fasta
asv_tab <- as.data.frame(seqtab.nochim)
colnames(asv_tab) <- asv_ids
asv_seqs <- DNAStringSet(asv_seq_chr)
names(asv_seqs) <- asv_ids

asv_tab: matriks kelimpahan (samples × ASV).
asv_seqs: FASTA berisi sekuens tiap ASV dengan nama ASV.

# 99.6 Write outputs
write.csv(per_sample,      file.path(outdir, "99_per_sample_summary.csv"), row.names = FALSE)
write.csv(asv_stats,       file.path(outdir, "99_asv_length_and_chimera_stats.csv"), row.names = FALSE)
write.csv(asv_len_dist,    file.path(outdir, "99_asv_length_distribution.csv"),    row.names = FALSE)
write.csv(tax_df,          file.path(outdir, "99_taxonomy_table.csv"),             row.names = FALSE)
write.csv(tax_with_boot,   file.path(outdir, "99_taxonomy_table_with_boot.csv"),   row.names = FALSE)
write.csv(asv_tab,         file.path(outdir, "ASV_table.csv"))
Biostrings::writeXStringSet(asv_seqs, file.path(outdir, "ASV_sequences.fasta"))

Semua output ditulis ke folder dada2_outputs/.

Ringkasan Output Utama

99_per_sample_summary.csv: ringkasan pipeline per sampel.
99_asv_length_and_chimera_stats.csv: statistik global panjang ASV & chimera.
99_asv_length_distribution.csv: distribusi panjang ASV.
99_taxonomy_table.csv & 99_taxonomy_table_with_boot.csv: tabel taksonomi.
ASV_table.csv: matriks kelimpahan ASV (samples × ASVs).
ASV_sequences.fasta: sekuens FASTA tiap ASV.

📊

Dengan output ini, Anda bisa langsung lanjut ke phyloseq atau mia/miaViz untuk analisis lanjutan.

🎉 Pipeline DADA2 ini menghasilkan ASV dengan resolusi nukleotida tunggal, taksonomi, dan ringkasan statistik siap untuk analisis ekologi mikroba.

2 Trimming Adaptor dan Filtering (Cutadapt)4 Phyloseq