# Metagenomics Amplicon

> Compute 16S/ITS amplicon diversity (Shannon, Simpson, Bray-Curtis, UniFrac), PCoA/NMDS ordination, PERMANOVA on OTU/ASV tables from QIIME2/DADA2. Use for 16S microbiome analysis.

- Skill: `pavel-kravchenko/metagenomics-amplicon` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/metagenomics-amplicon`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/metagenomics-amplicon/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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).

```python
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.

```python
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.

```bash
# 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.

```r
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