Metagenomics Amplicon

16S/ITS Amplicon Metagenomics

When to Use

• Analyzing a 16S rRNA (or ITS fungal) amplicon feature table (OTU/ASV counts x samples) from QIIME2, DADA2, or mothur.
• Comparing microbial community diversity within samples (alpha) or between samples/groups (beta).
• Testing whether community composition differs significantly across experimental groups (PERMANOVA).
• Deciding between OTU clustering and ASV denoising, or interpreting a DADA2/QIIME2 pipeline's output.
• Visualizing community structure via PCoA/NMDS ordination or taxonomic bar plots.

Version Compatibility

• QIIME2 ≥ 2024.2, DADA2 ≥ 1.30 (R/Bioconductor), scikit-bio ≥ 0.6, Python ≥ 3.10
• SILVA 138 / Greengenes2 as reference taxonomy databases

Prerequisites

• pip install scikit-bio pandas numpy scipy scikit-learn (scikit-bio has native Shannon/Simpson/UniFrac/PCoA/PERMANOVA — prefer it over hand-rolled code in production; the functions below are for when scikit-bio isn't available or you need to see the math)
• A feature table: rows = OTUs/ASVs, columns = samples, values = read counts
• Sample metadata (grouping variable) and, for phylogenetic metrics, a rooted tree of the ASVs

Background: OTU vs ASV

• OTU (97% similarity clustering, e.g. VSEARCH/UCLUST): loses within-cluster variation, not reproducible across studies.
• ASV (DADA2 exact denoising): single-nucleotide resolution, error-corrected, reproducible — the current standard (Callahan et al. 2017, "Exact sequence variants should replace OTUs").

Alpha Diversity (within-sample)

Goal: quantify richness and evenness of a single sample's community. Approach: compute from relative abundances; always check whether samples need rarefying to equal depth first (raw alpha diversity is confounded by sequencing depth).

import numpy as np
import pandas as pd

def observed_species(counts):
    """Richness: count of features with count > 0."""
    return int(np.sum(np.asarray(counts) > 0))

def shannon_diversity(counts):
    """Shannon index H' = -sum(p_i * ln(p_i)); richness + evenness."""
    counts = np.asarray(counts, dtype=float)
    p = counts[counts > 0] / counts[counts > 0].sum()
    return -np.sum(p * np.log(p))

def simpson_diversity(counts):
    """Simpson's 1-D = 1 - sum(p_i^2); robust to rare taxa, higher = more diverse."""
    counts = np.asarray(counts, dtype=float)
    p = counts[counts > 0] / counts[counts > 0].sum()
    return 1 - np.sum(p ** 2)

def pielou_evenness(counts):
    """Pielou's J' = H' / ln(S); 1.0 means perfectly even abundances."""
    s = observed_species(counts)
    if s <= 1:
        return np.nan
    return shannon_diversity(counts) / np.log(s)

def rarefy(counts, depth, rng=None):
    """Subsample a sample's counts down to `depth` total reads without replacement."""
    rng = rng or np.random.default_rng(0)
    counts = np.asarray(counts, dtype=int)
    reads = np.repeat(np.arange(len(counts)), counts)
    if len(reads) < depth:
        raise ValueError(f"sample has only {len(reads)} reads, cannot rarefy to {depth}")
    sub = rng.choice(reads, size=depth, replace=False)
    return np.bincount(sub, minlength=len(counts))

# Apply across a feature_table DataFrame (features x samples)
def alpha_diversity_table(feature_table):
    """Compute Observed, Shannon, Simpson, Pielou for every sample column."""
    return pd.DataFrame({
        "observed": feature_table.apply(observed_species, axis=0),
        "shannon": feature_table.apply(shannon_diversity, axis=0),
        "simpson": feature_table.apply(simpson_diversity, axis=0),
        "pielou_evenness": feature_table.apply(pielou_evenness, axis=0),
    })

Beta Diversity, Ordination, and PERMANOVA

Goal: measure between-sample dissimilarity and test whether groups differ. Approach: build a distance matrix (Bray-Curtis for abundance, Jaccard/UniFrac for presence-absence or phylogeny-aware), reduce with PCoA for visualization, then test group separation with PERMANOVA.

import numpy as np
import pandas as pd
from itertools import combinations

def bray_curtis(s1, s2):
    """Bray-Curtis dissimilarity: 0 = identical, 1 = no shared taxa."""
    s1, s2 = np.asarray(s1, float), np.asarray(s2, float)
    denom = np.sum(s1 + s2)
    return np.sum(np.abs(s1 - s2)) / denom if denom > 0 else 0.0

def distance_matrix(feature_table, metric=bray_curtis):
    """Build a symmetric sample x sample distance matrix from a features x samples table."""
    samples = feature_table.columns
    dm = pd.DataFrame(0.0, index=samples, columns=samples)
    for a, b in combinations(samples, 2):
        d = metric(feature_table[a], feature_table[b])
        dm.loc[a, b] = dm.loc[b, a] = d
    return dm

def pcoa(dm_df):
    """Classical MDS (Principal Coordinates Analysis) on a distance matrix DataFrame.
    Returns (coords DataFrame with PC1/PC2, proportion of variance explained per axis)."""
    dm = dm_df.values
    n = len(dm)
    H = np.eye(n) - np.ones((n, n)) / n
    B = -0.5 * H @ (dm ** 2) @ H
    eigvals, eigvecs = np.linalg.eigh(B)
    idx = np.argsort(eigvals)[::-1]
    eigvals, eigvecs = eigvals[idx], eigvecs[:, idx]
    pos = eigvals > 1e-10
    coords = eigvecs[:, pos] * np.sqrt(eigvals[pos])
    prop = eigvals[pos] / eigvals[pos].sum()
    return pd.DataFrame(coords[:, :2], index=dm_df.index, columns=["PC1", "PC2"]), prop

def permanova(dm, grouping, n_perm=999, seed=0):
    """Permutational MANOVA: tests whether group centroids differ in the distance matrix.
    `grouping` is a Series indexed by sample id -> group label. Returns (pseudo-F, p-value)."""
    rng = np.random.default_rng(seed)
    dm_v = dm.values
    groups = grouping.loc[dm.index].values
    unique, n = np.unique(groups), len(groups)

    def f_stat(g):
        ss_tot = np.sum(dm_v ** 2) / n
        ss_w = sum(
            np.sum(dm_v[np.ix_(g == grp, g == grp)] ** 2) / (2 * (g == grp).sum())
            for grp in unique if (g == grp).sum() > 1
        )
        df_b, df_w = len(unique) - 1, n - len(unique)
        return ((ss_tot - ss_w) / df_b) / (ss_w / df_w) if df_w > 0 and ss_w > 0 else 0.0

    obs = f_stat(groups)
    perm_stats = [f_stat(rng.permutation(groups)) for _ in range(n_perm)]
    p = (np.sum(np.array(perm_stats) >= obs) + 1) / (n_perm + 1)
    return obs, p

QIIME2/DADA2 Pipeline (command reference)

Goal: go from raw paired-end FASTQ to a taxonomy-annotated, phylogeny-aware feature table. Approach: import → denoise → classify → build tree → compute core diversity metrics.

# 1. Import paired-end reads via a manifest CSV (sample-id,absolute-filepath,direction)
qiime tools import --type 'SampleData[PairedEndSequencesWithQuality]' \
    --input-path manifest.csv --input-format PairedEndFastqManifestPhred33V2 \
    --output-path demux.qza
qiime demux summarize --i-data demux.qza --o-visualization demux.qzv

# 2. Denoise with DADA2 -> ASV feature table + representative sequences
qiime dada2 denoise-paired \
    --i-demultiplexed-seqs demux.qza \
    --p-trim-left-f 0 --p-trim-left-r 0 \
    --p-trunc-len-f 240 --p-trunc-len-r 200 \
    --o-table table.qza --o-representative-sequences rep-seqs.qza \
    --o-denoising-stats denoising-stats.qza

# 3. Classify ASVs against SILVA 138 (pre-trained Naive Bayes classifier for your primer pair)
qiime feature-classifier classify-sklearn \
    --i-classifier silva-138-99-nb-515-806-classifier.qza \
    --i-reads rep-seqs.qza --o-classification taxonomy.qza

# 4. Build phylogenetic tree (needed for Faith's PD and UniFrac)
qiime phylogeny align-to-tree-mafft-fasttree \
    --i-sequences rep-seqs.qza \
    --o-alignment aligned.qza --o-masked-alignment masked.qza \
    --o-tree unrooted-tree.qza --o-rooted-tree rooted-tree.qza

# 5. Core diversity metrics (rarefies to --p-sampling-depth internally)
qiime diversity core-metrics-phylogenetic \
    --i-phylogeny rooted-tree.qza --i-table table.qza \
    --p-sampling-depth 10000 --m-metadata-file metadata.tsv \
    --output-dir core-metrics-results

# 6. Differential abundance (compositionally-aware, corrects for sampling fraction bias)
qiime composition ancombc --i-table table.qza --m-metadata-file metadata.tsv \
    --p-formula group --o-differentials ancombc-diff.qza

Differential Abundance in R (compositional data)

Goal: find taxa that differ between groups without ignoring the compositional nature of relative-abundance data. Approach: ANCOM-BC (Bioconductor) models sampling-fraction bias explicitly; a Kruskal-Wallis + BH fallback works for a quick screen.

library(ANCOMBC)
library(phyloseq)

# ps: a phyloseq object built from the QIIME2 feature table, taxonomy, tree, and metadata
out <- ancombc2(
  data = ps, fix_formula = "group", p_adj_method = "BH",
  group = "group", struc_zero = TRUE, neg_lb = TRUE
)
sig_taxa <- out$res[out$res$diff_group == TRUE, ]

Pitfalls

• Alpha diversity without rarefying — deeper-sequenced samples appear more diverse; rarefy to a common depth (or use rarefaction curves to confirm plateau) before comparing across samples.
• PCoA axes without variance-explained — PC1 may only explain 15-20%; always report the prop values, don't over-interpret a 2D plot.
• PERMANOVA detects dispersion, not just centroids — pair with a betadisper (homogeneity of multivariate dispersion) test in R's vegan to rule out a dispersion artifact.
• Standard t-tests/ANOVA on relative abundances — compositional data (sums to 1) violates independence assumptions; use ANCOM-BC, ALDEx2, or at minimum a non-parametric test with BH correction.
• OTU (97%) tables from old pipelines are not comparable to ASV tables — don't merge them across studies; re-process raw reads through DADA2 if you need cross-study comparability.
• Missing chimera/mitochondria/chloroplast filtering — DADA2 removes chimeras, but classify and filter out host mitochondrial/chloroplast 16S hits before diversity analysis.

See Also

• bio-microbiome-qiime2-workflow — full QIIME2 CLI pipeline details
• bio-microbiome-diversity-analysis — deeper alpha/beta diversity statistics
• bio-microbiome-differential-abundance — ANCOM-BC/ALDeX2 in depth
• metagenomics-shotgun — shotgun (non-amplicon) taxonomic and functional profiling