Step 4 — DADA2 Denoising

Script: scripts/mbx_dada2_run.sh

Companion files in this folder: - 4_dada2_run.html — same content with copy buttons on every code block. - 4_dada2_run.pptx — slide deck for the talk.

Why this step matters

By the time we reach Step 4, the data is just reads — short DNA strings from the sequencer, each carrying its own per-base quality scores. Reads are noisy. Illumina makes about 1 wrong base per 1,000 in the best part of a read, and worse near the back. If we treated every distinct read as a separate organism, our results would be hopelessly contaminated by sequencing errors.

DADA2 (Divisive Amplicon Denoising Algorithm) solves this. Instead of clustering reads at some arbitrary similarity threshold (the old OTU approach), DADA2 models the actual error process and asks: given how Illumina makes mistakes, is this distinct read a real biological variant — or is it just a slightly-corrupted copy of one?

The output is a list of Amplicon Sequence Variants (ASVs) — exactly-resolved DNA sequences that survive the error model. ASVs are single-nucleotide resolved, comparable across studies, and far more defensible than OTUs.

This step's job is to actually run DADA2 on our paired-end data, using the truncation parameters Step 3 already computed, and to produce the three artifacts every downstream step needs.

What the script does in one sentence

It reads the four --p-* numbers from Step 3's dada2_parameters.txt, runs qiime dada2 denoise-paired with those numbers and a sensible thread count, then exports the three QZA outputs as human-readable QZV visualisations.

The algorithm, step by step

1. Validate the metadata file

First, the script reads the user's metadata file and confirms it has the structure DADA2 (and every later step) expects:

• A header row.
• A first column named sample-id, sampleid, id, or one of the QIIME2- recognised synonyms.
• Sample IDs that match the artifact (every sample DADA2 will process appears in the metadata).

A mismatch here would let DADA2 succeed but produce results no downstream step could merge against the metadata. So we catch it before DADA2 even starts.

2. Parse the DADA2 parameter file

Then the script reads 3_dada2_parameters/dada2_parameters.txt and pulls out the four numbers Step 3 picked:

• --p-trim-left-f — bases to chop off the front of the forward read.
• --p-trim-left-r — bases to chop off the front of the reverse read.
• --p-trunc-len-f — length to keep on the forward read.
• --p-trunc-len-r — length to keep on the reverse read.

It also reads INSERT_LENGTH_USED and any FALLBACK_REASONS from the same file — those don't change DADA2's behaviour, but they're logged so the final report can show "Step 3 had to relax its filters because …".

3. Auto-detect the thread count

Next the script decides how many CPU threads to give DADA2:

• If the orchestrator exported MBX_THREADS (because the user passed mbXPro --threads N), use that.
• Otherwise call nproc (Linux) or sysctl -n hw.logicalcpu (macOS).
• Last-resort fallback: 1 thread.

DADA2 is embarrassingly parallel across samples — denoising sample A doesn't need anything from sample B — so the more threads, the faster.

4. Run the denoise-paired plugin

Now the actual denoising. The command:

qiime dada2 denoise-paired \
  --i-demultiplexed-seqs   Paired_End_artifact.qza \
  --p-trim-left-f          <from step 3> \
  --p-trim-left-r          <from step 3> \
  --p-trunc-len-f          <from step 3> \
  --p-trunc-len-r          <from step 3> \
  --p-n-threads            <auto-detected> \
  --o-table                feature_table.qza \
  --o-representative-sequences  representative_sequences.qza \
  --o-denoising-stats     dada2_stats.qza

Inside QIIME2 this triggers the DADA2 R algorithm, which:

1. Filters and trims each read to the requested truncation length and removes pairs whose expected error exceeds max_ee (2.0 by default).
2. Learns the error model from the filtered reads themselves — DADA2 estimates the substitution probabilities P(true base → observed base | quality) empirically. No assumptions baked in.
3. Denoises each sample by asking "is this distinct read more likely a real ASV, or a corrupted copy of a more abundant ASV?", iterated until stable.
4. Merges the forward and reverse reads (using the overlap Step 3 guaranteed would exist).
5. Removes chimeras via the consensus method — sequences that look like a PCR-stitched combination of two real ASVs are dropped.

The output is three artifacts:

• feature_table.qza — counts of every ASV in every sample.
• representative_sequences.qza — the actual DNA sequence of every ASV.
• dada2_stats.qza — how many reads each sample lost at each stage (input → filtered → denoised → merged → non-chimeric).

5. Produce three QZV visualisations

Finally the script runs three QIIME2 visualisation commands so the user can inspect what happened without having to extract the QZAs:

• qiime metadata tabulate → dada2_stats.qzv — interactive table of the per-sample retention numbers.
• qiime feature-table summarize → feature_table.qzv — per-sample read-count distribution, alpha-rarefaction preview, sample frequency histogram.
• qiime feature-table tabulate-seqs → representative_sequences.qzv — searchable table of ASV sequences with one-click BLAST links.

All three drop into 4_dada2_outputs/.

Default parameters and why they are what they are

When and why we fall back to defaults

What the output file looks like

dada2_stats.qzv (after qiime metadata tabulate) shows one row per sample:

A healthy run typically keeps ≥ 70 % of the input reads through to non-chimeric. Big drops (≥ 50 %) at any one stage are diagnostic:

• Big drop at filtered → truncation lengths too aggressive or max_ee too strict.
• Big drop at merged → forward + reverse cuts don't leave enough overlap to merge (Step 3 should have caught this).
• Big drop at non-chimeric → PCR was over-cycled, the library has lots of chimeras.

The final report (Step 18) flags any sample that drops below 50 % retention so the user knows where to look.

Takeaway

Step 4 is the step where reads become biology. DADA2 models the actual Illumina error process and resolves single-nucleotide variants instead of clustering at an arbitrary threshold. The mbX Pro wrapper just makes sure the parameters Step 3 computed actually get used, the metadata matches the artifact, and the three follow-up visualisations are produced so the user can sanity-check the per-sample retention before the next step consumes the result.

Sources

• The script: mbXPro/scripts/mbx_dada2_run.sh
• DADA2 paper: Callahan et al. (2016), DADA2: high-resolution sample inference from Illumina amplicon data, Nature Methods 13:581–583.
• QIIME2 q2-dada2 plugin: https://docs.qiime2.org/2025.4/plugins/available/dada2/
• Why ASVs over OTUs: Callahan et al. (2017), Exact sequence variants should replace operational taxonomic units in marker-gene data analysis, ISME Journal 11:2639–2643.