16S amplicon analysis · illustrated

From reads to a result

A 16S study passes through six or seven programs before it produces a sentence anyone would publish. Each one has defaults, each default encodes an assumption, and almost none of them announce themselves. This site walks the whole chain — the SILVA reference database, QIIME 2 and DADA2, alpha and beta diversity, and MaAsLin — with the decision points marked.

Every figure is animated and steppable. The read-attrition numbers are real, from a 50-sample run. There is also a game about the decisions.

The whole chain

Eight stages. The animation builds them in order; the caption for each names the decision it hides.

1 · raw reads run 1 · run 2 — different error profiles 2 · import → .qza .qza action parameters versions 3 · DADA2 denoise 100 → 92.3 → 91.2 → 83.9 → 62.7 % 4 · ASV table ASVs × samples 5 · taxonomy SILVA d__Bacteria p__Bacillota g__Blautia train on your amplicon region 6 · phylogeny UniFrac Faith's PD only needed for phylogenetic metrics 7 · diversity + stats PERMANOVA + PERMDISP 8 · MaAsLin coefficient · q-value default is q < 0.25

Six places the answer is decided

Not bugs — defaults. Each of these silently determines part of the result, and each has a defensible setting that is not the one you get for free.

  1. Which reference, and trained on which region A naive-Bayes classifier trained on full-length 16S and applied to a 250 bp V4 read is measurably worse than one trained on the V4 region itself. SILVA now ships region-specific classifiers, and RESCRIPt will build one for any primer pair. → SILVA
  2. Truncation lengths The single most destructive parameter in the pipeline. Truncate to clean up quality and you can leave the read pair without enough overlap to merge — losing most of your data at a step that reports success. → DADA2
  3. One error model or one per run DADA2 learns error rates from the data. Error profiles differ between runs, so a multi-run study needs one model per run and a merged table afterwards. → DADA2
  4. How sequencing depth is handled In the real run below, final depth ranges from 66k to 174k reads — a 2.6× spread that maps directly onto observed richness if you do nothing about it. → Diversity
  5. Whether you tested dispersion PERMANOVA returns a significant p-value both when two groups sit in different places and when one group is simply more variable. Only PERMDISP distinguishes them. → Diversity
  6. Normalisation, transform, and what counts as significant MaAsLin's defaults are TSS + log + a linear model, and its default significance threshold is q < 0.25. All three are reasonable choices that need to be stated out loud. → MaAsLin

And one that is decided before the pipeline starts at all: whether contamination was controlled. A constant reagent DNA load becomes a larger fraction of a smaller sample, so in any study where biomass differs between groups the contaminant correlates with the treatment. → The kitome

Where the reads actually go

From a real 50-sample V4 run: 9,973,604 input reads, 6,256,824 surviving to the feature table. Both series are percentages of each sample's own input, so they are directly comparable.

Read survival through DADA2

50 samples, QIIME 2 2024.10.1. Percent of each sample's own input reads.

median sample worst sample

Two things worth taking from this. First, chimera removal is the largest single loss — 21% of all input reads, more than quality filtering and merging combined. That is normal for V4 and not a sign of anything wrong. Second, the gap between the median and the worst sample widens at every stage: 0.8 points after filtering, 6.5 after merging, 17.7 after chimera removal. Attrition is not a constant tax, and the samples that lose most are not random — which is precisely why the depth you end up with needs handling rather than ignoring.