Bio microbiome diversity analysis skills diversity analysis

Version Compatibility

Reference examples tested with: phyloseq 1.46+, vegan 2.6+, picante 1.8+, GUniFrac 1.8+, scikit-bio 0.6+, QIIME2 2024.2+.

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: qiime <plugin> <action> --help to confirm flags
• Python: pip show scikit-bio then help(skbio.diversity.beta_diversity) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

scikit-bio 0.6.0 renamed OTU to taxon across the API and drifted metric kwargs (otu_ids= vs newer forms) - discover names with skbio.diversity.get_beta_diversity_metrics() before hard-coding. UniFrac/Faith PD results inherit the tree (de novo vs SEPP vs Greengenes2 reference build) AND the chosen sampling depth - record both alongside the QIIME2 release that produced the .qza artifacts.

Diversity Analysis

"Compare microbial diversity across my samples" -> Summarize within-sample richness/evenness (alpha) and between-sample dissimilarity (beta) - but only after declaring the rarefaction depth, the tree, and the metric, because each is a knob that sets the answer before it is seen.

• CLI: qiime diversity core-metrics-phylogenetic --i-phylogeny rooted-tree.qza --i-table table.qza --p-sampling-depth N --m-metadata-file md.tsv --output-dir cm/
• R: estimate_richness(ps_rare) for alpha; UniFrac(ps_rare, weighted=) / vegdist() then adonis2() + betadisper() for beta

Scope: whole-community summary (a number or ordination per sample) of an amplicon ASV/OTU table plus a tree. Per-taxon between-group testing -> differential-abundance. Shotgun profiler tables (MetaPhlAn/Bracken) -> metagenomics/metagenome-visualization. The shared CoDA and rarefaction-debate theory lives in metagenomics/abundance-estimation; the Hill-number and PERMANOVA-dispersion theory in metagenomics/metagenome-visualization - cross-referenced here, not re-derived. Tree handling -> phylogenetics/tree-io.

The Single Most Important Modern Insight -- A Diversity Number Is the Output of Three Knobs Turned Before the Answer Appears

An alpha or beta diversity value is not a measurement of the community; it is the output of three choices made before the number appears - the rarefaction DEPTH, the TREE, and the METRIC. Turn them differently and the conclusion can change. The job is to declare all three and show the result survives a second reasonable choice, not to run core-metrics-phylogenetic and read the p-value. The quietest and most dangerous knob is the depth:

1. --p-sampling-depth is a sample-deletion knob in a normalization costume. core-metrics rarefies every sample to the depth and, per the QIIME2 docs, silently drops every sample whose total count is below it - no warning, just fewer points in the PCoA. The dropped samples are the lowest-yield ones (the lowest-biomass swab, the sickest patient, the failed extraction), so the loss is almost never random. Pick the depth from the feature-table summary plus the alpha-rarefaction plateau, report the depth AND the dropped samples, and confirm the conclusion at a nearby depth. Rarefying to min(sample_sums) is the worst of both worlds - one tiny library drags everyone to noise.
2. UniFrac and Faith PD are only as real as the tree, and the tree is a model not a property of the data. A de novo MAFFT+FastTree tree from ~250 bp reads is poorly resolved and arbitrarily midpoint-rooted (Janssen 2018); SEPP fragment-insertion into a full-length reference, or Greengenes2 placement, gives stable topology and correct associations - and SEPP/GG2 align 16S with shotgun (McDonald 2024). SEPP also drops fragments that fail to insert, a second silent table-shrink.
3. Rarefy for diversity, never for differential abundance. Rarefaction-to-even-depth is defensible for alpha/beta (Schloss 2024); for DA it discards count information a compositional model needs (McMurdie 2014). Keep the raw counts; rarefy only into the diversity branch; route DA to differential-abundance on the unrarefied table.

Tool / Metric Taxonomy

Metric / tool	Citation	What it measures / does	When
Observed features	-	ASV richness (Hill q=0); most depth-sensitive; an ASV count, not species	richness, but report denoising params; prefer Hill q1/q2
Shannon	-	entropy = richness+evenness (Hill q=1 = exp(H')); QIIME2 log2/bits, R ln/nats	balanced diversity; report exp(H') to dodge the base
Pielou evenness	Pielou 1966 J Theor Biol 13:131	H'/ln(S); 0-1; isolates evenness from richness	when evenness is the question
Faith PD	Faith 1992 Biol Conserv 61:1	sum of branch lengths spanning observed taxa; phylogenetic q=0	amplicon-native richness; needs a tree
Jaccard	-	presence/absence dissimilarity; no tree	membership turnover; depth/rare-ASV sensitive
Bray-Curtis	-	abundance dissimilarity; no tree; compositionally incoherent	abundance default; intuitive, label the caveat
Unweighted UniFrac	Lozupone 2005 Appl Environ Microbiol 71:8228	branch length unique to one community (presence/absence)	rare/divergent lineages + topology; needs a tree
Weighted UniFrac	Lozupone 2007 Appl Environ Microbiol 73:1576	branch length weighted by abundance difference	abundant-lineage shifts; needs a tree
Generalized UniFrac	Chen 2012 Bioinformatics 28:2106	alpha in [0,1] interpolating unweighted-weighted	alpha=0.5 compromise; powerful for moderately abundant lineages
Aitchison / RPCA	Martino 2019 mSystems 4:e00016-19	CLR + matrix completion; ordination with feature loadings	compositionally coherent; sparse data; no pseudocount
SEPP insertion	Janssen 2018 mSystems 3:e00021-18	places ASVs into a full-length reference tree	the preferred tree for short reads
Greengenes2	McDonald 2024 Nat Biotechnol 42:715	unified genome+16S reference tree	makes 16S UniFrac comparable to shotgun

Decision Tree by Scenario

Scenario	Recommended	Why
Need a phylogenetic metric (UniFrac, Faith PD)	SEPP-into-reference or Greengenes2 tree	de novo from short reads is unstable (Janssen 2018)
De novo tree is the only option	treat unweighted UniFrac with suspicion	topology noise on ~250 bp reads dominates it
Change is in rare/low-abundance lineages	unweighted UniFrac, observed features	presence/absence + topology see rare taxa
Change is a bloom of dominant taxa	weighted UniFrac, Bray-Curtis	abundance-weighted metrics see dominant shifts
Do not want to metric-shop	generalized UniFrac alpha=0.5 + report both un/weighted	Chen 2012 compromise; single-metric hit is tentative
Richness vs evenness question	observed/Faith (q0) AND Shannon-exp (q1) / InvSimpson (q2)	span the richness-evenness spectrum
Compositional, want axis-driving taxa	RPCA (DEICODE/gemelli)	CLR ordination with interpretable loadings
Picking a rarefaction depth	feature-table summarize + alpha-rarefaction plateau	depth must retain samples AND saturate richness
Per-taxon "which bug changed"	-> differential-abundance	diversity is whole-community; DA is per-feature
Shotgun profiler table, not amplicon	-> metagenomics/metagenome-visualization	no per-feature tree; different idiom

Choosing the Sampling Depth (the biggest lever)

Goal: Pick a rarefaction depth that saturates richness while retaining an acceptable fraction of samples, and know exactly which samples were dropped.

Approach: Read the per-sample frequency distribution from the feature-table summary, find where the alpha-rarefaction curve plateaus, set the depth there, then declare the depth and the dropped-sample list.

qiime feature-table summarize --i-table table.qza --o-visualization table.qzv   # per-sample frequencies; the depth lives here

qiime diversity alpha-rarefaction \
    --i-table table.qza --i-phylogeny rooted-tree.qza \
    --p-max-depth 20000 \   # set near the median sample depth; the curve panel shows survivors per depth
    --m-metadata-file metadata.tsv --o-visualization alpha-rarefaction.qzv

qiime diversity core-metrics-phylogenetic \
    --i-phylogeny rooted-tree.qza --i-table table.qza \
    --p-sampling-depth 10000 \   # on the observed-features plateau; SILENTLY DROPS samples below this
    --m-metadata-file metadata.tsv --output-dir core-metrics-results

core-metrics-phylogenetic rarefies the table, computes the four alpha vectors (faith_pd_vector, observed_features_vector, shannon_vector, evenness_vector) and four beta matrices (unweighted_unifrac_, weighted_unifrac_, jaccard_, bray_curtis_distance_matrix), and produces a PCoA + Emperor plot for each beta metric. The non-phylogenetic twin qiime diversity core-metrics drops Faith PD and both UniFracs and needs no tree.

Building the Tree (a modeling choice, not a fixed step)

Goal: Obtain a phylogeny over the ASVs that does not inject topology noise into UniFrac/Faith PD.

Approach: Prefer SEPP fragment-insertion into a full-length reference (or Greengenes2 placement) over a de novo build from short reads; for de novo, mask the alignment and accept that unweighted UniFrac will be shaky.

qiime fragment-insertion sepp \
    --i-representative-sequences rep-seqs.qza \
    --i-reference-database sepp-refs-gg-13-8.qza \
    --p-threads 4 \
    --o-tree insertion-tree.qza --o-placements insertion-placements.qza

qiime fragment-insertion filter-features \
    --i-table table.qza --i-tree insertion-tree.qza \
    --o-filtered-table table-sepp.qza --o-removed-table removed-table.qza   # fragments that failed to insert are DROPPED

De novo is qiime phylogeny align-to-tree-mafft-fasttree (MAFFT align -> mask -> FastTree2 -> midpoint root) - acceptable only when no reference package fits the marker/region, and unweighted UniFrac on it must be treated as suspect.

Alpha Diversity in R (counts on the rarefied table)

Goal: Compute richness and evenness per sample and test for a group difference without confounding by sequencing depth.

Approach: Rarefy to a chosen depth, estimate Hill-spanning metrics, test with a non-parametric test (escalate to a linear/mixed model for covariates), and report effective species exp(H').

library(phyloseq); library(vegan)

ps_rare <- rarefy_even_depth(ps, sample.size = chosen_depth, rngseed = 42, replace = FALSE)
alpha <- estimate_richness(ps_rare, measures = c('Observed', 'Shannon', 'InvSimpson'))   # q0, exp gives q1, q2
alpha$Group <- sample_data(ps_rare)$Group
alpha$Shannon_eff <- exp(alpha$Shannon)   # effective species; base-invariant in interpretation (Hill q=1)

kruskal.test(Shannon ~ Group, data = alpha)   # non-parametric; escalate to lme4/nlme for covariates or repeated measures

Faith PD in R uses picante::pd(otu_matrix, tree, include.root = TRUE). The Shannon from estimate_richness is in natural log (nats); QIIME2 reports log2 (bits) - report exp(Shannon) to compare across the two.

Beta Diversity in R (report weighted AND unweighted)

Goal: Quantify between-sample dissimilarity with phylogenetic and abundance-weighted views, then test the group effect while ruling out a dispersion artifact.

Approach: Compute both UniFrac variants (and generalized UniFrac alpha=0.5), ordinate by PCoA, run adonis2 for location, and ALWAYS pair it with betadisper for spread.

wu  <- UniFrac(ps_rare, weighted = TRUE)    # abundant-lineage view
uwu <- UniFrac(ps_rare, weighted = FALSE)   # rare-lineage + topology view
# generalized UniFrac alpha=0.5 (Chen 2012 compromise):
gu  <- as.dist(GUniFrac::GUniFrac(t(as(otu_table(ps_rare), 'matrix')), phy_tree(ps_rare), alpha = 0.5)$unifracs[, , 'd_0.5'])

meta <- data.frame(sample_data(ps_rare))
adonis2(wu ~ Group, data = meta, permutations = 999)   # >=999 permutations; significance = LOCATION
permutest(betadisper(wu, meta$Group))                  # MANDATORY: is it dispersion, not location?

If betadisper is significant the adonis2 result is ambiguous (location vs spread) - state it. The PERMANOVA-dispersion theory is shared; see metagenomics/metagenome-visualization. For a compositionally coherent ordination with feature loadings use RPCA (DEICODE qiime deicode rpca / gemelli). The Python engine is scikit-bio (skbio.diversity.beta_diversity, skbio.stats.ordination.pcoa, skbio.stats.distance.permanova).

Per-Method Failure Modes

Sampling-depth sample-massacre

Trigger: a --p-sampling-depth higher than some samples' totals. Mechanism: core-metrics drops every sample below the depth with no warning. Symptom: fewer points in the PCoA than samples in the metadata; the lost ones skew low-biomass. Fix: pick the depth from the rarefaction plateau, report the dropped-sample list, confirm at a nearby depth.

De novo tree noise

Trigger: UniFrac/Faith PD on a MAFFT+FastTree tree from short reads. Mechanism: ~250 bp reads give an unstable topology and arbitrary midpoint root. Symptom: unweighted-UniFrac separation that vanishes under SEPP insertion or weighted UniFrac. Fix: use SEPP-into-reference or Greengenes2; treat de novo unweighted UniFrac as suspect.

Unweighted-vs-weighted flip

Trigger: reporting only the UniFrac variant that gives p<0.05. Mechanism: unweighted listens to rare/short branches, weighted to abundant lineages. Symptom: the two disagree and the chosen one is the significant one. Fix: report both plus generalized alpha=0.5; state which lineage axis each implicates.

Rarefy-then-reuse-for-DA

Trigger: feeding the rarefied table to a differential-abundance tool. Mechanism: rarefaction discards count information the DA model needs. Symptom: underpowered or distorted DA. Fix: keep raw counts; rarefy only into the diversity branch; route DA to differential-abundance.

Observed-features-as-species

Trigger: comparing raw ASV counts across runs/studies as "richness". Mechanism: ASV count tracks DADA2 truncation/maxEE/pooling and intragenomic 16S copy variants, not just biology. Symptom: richness shifts with denoising settings. Fix: prefer Hill q1/q2; report observed features with the denoising parameters stated.

Shannon base mismatch

Trigger: comparing a QIIME2 Shannon to an R Shannon. Mechanism: QIIME2 uses log2 (bits), R diversity/estimate_richness natural log (nats). Symptom: numbers differ by a constant factor and look like a real effect. Fix: state the base, convert, or report exp(H').

PERMANOVA dispersion

Trigger: a significant adonis2 read as a composition shift. Mechanism: pseudo-F responds to within-group spread, not only centroid location (shared theory; metagenomics/metagenome-visualization). Symptom: significant adonis2 with significant betadisper. Fix: always run betadisper/permutest alongside; report both.

Quantitative Thresholds

Threshold	Source	Rationale
Sampling depth on the observed-features plateau	Janssen 2018; QIIME2 docs	depth must saturate richness while retaining samples; report dropped list
Do NOT use min(sample_sums) as the depth	McMurdie 2014	one tiny library drags every sample to under-saturated noise
Generalized UniFrac alpha = 0.5	Chen 2012 Bioinformatics 28:2106	most powerful for moderately abundant lineages; beats running un/weighted jointly
Report Hill q = 0, 1, 2 together	(shared; metagenomics/metagenome-visualization)	spans richness (q0) -> evenness-weighted (q2)
PERMANOVA permutations >= 999	vegan docs	resolution floor for p ~ 0.001; use 9999 for publication
Pair adonis2 with betadisper	Anderson & Walsh 2013 (shared)	distinguishes a location shift from a dispersion difference
Rarefy for diversity, not for DA	McMurdie 2014; Schloss 2024	per-analysis decision, not a global switch

Common Errors

Error / symptom	Cause	Solution
PCoA has fewer points than samples	--p-sampling-depth dropped low-count samples	lower the depth or report the loss; never assume zero drops
UniFrac errors / Faith PD missing	no phy_tree slot in the phyloseq object	attach a SEPP/GG2 (preferred) or de novo tree
Unweighted UniFrac significant, weighted not	change is in rare lineages, or de novo tree noise	report both; verify the tree; treat single-metric hit as tentative
R and QIIME2 Shannon disagree	log base differs (nats vs bits)	report exp(H'); convert by log2(e)
adonis2 p<0.001 but groups visually overlap	dispersion difference, not location	run betadisper; report it
scikit-bio otu_ids= deprecation warning	0.6 renamed OTU to taxon; otu_ids= kept as a deprecated alias	get_beta_diversity_metrics() and help() to find current kwargs
Diversity tracks host/plant content	host mitochondria/chloroplast 16S not removed	filter Mitochondria/Chloroplast features (see taxonomy-assignment) before computing diversity
"Community" in a near-sterile/low-biomass sample	reagent kitome not removed	sequence controls + run decontam upstream (amplicon-processing; metagenomics/contamination-controls)

References

• Faith DP. 1992. Conservation evaluation and phylogenetic diversity. Biol Conserv 61:1-10.
• Pielou EC. 1966. The measurement of diversity in different types of biological collections. J Theor Biol 13:131-144.
• Lozupone C, Knight R. 2005. UniFrac: a new phylogenetic method for comparing microbial communities. Appl Environ Microbiol 71:8228-8235.
• Lozupone CA, Hamady M, Kelley ST, Knight R. 2007. Quantitative and qualitative beta diversity measures lead to different insights into factors that structure microbial communities. Appl Environ Microbiol 73:1576-1585.
• Chen J, Bittinger K, Charlson ES, Hoffmann C, Lewis J, Wu GD, Collman RG, Bushman FD, Li H. 2012. Associating microbiome composition with environmental covariates using generalized UniFrac distances. Bioinformatics 28:2106-2113.
• Janssen S, McDonald D, Gonzalez A, et al. 2018. Phylogenetic placement of exact amplicon sequences improves associations with clinical information. mSystems 3:e00021-18.
• Mirarab S, Nguyen N, Warnow T. 2012. SEPP: SATe-enabled phylogenetic placement. Pac Symp Biocomput 2012:247-258.
• McDonald D, Jiang Y, Balaban M, et al. 2024. Greengenes2 unifies microbial data in a single reference tree. Nat Biotechnol 42:715-718.
• Martino C, Morton JT, Marotz CA, Thompson LR, Tripathi A, Knight R, Zengler K. 2019. A novel sparse compositional technique reveals microbial perturbations. mSystems 4:e00016-19.
• McDonald D, Vazquez-Baeza Y, Koslicki D, et al. 2018. Striped UniFrac: enabling microbiome analysis at unprecedented scale. Nat Methods 15:847-848.
• McMurdie PJ, Holmes S. 2014. Waste not, want not: why rarefying microbiome data is inadmissible. PLoS Comput Biol 10:e1003531.
• Schloss PD. 2024. Rarefaction is currently the best approach to control for uneven sequencing effort in amplicon sequence analyses. mSphere 9:e00354-23.
• McMurdie PJ, Holmes S. 2013. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 8:e61217.

Related Skills

• amplicon-processing - Generate the ASV table and representative sequences upstream
• taxonomy-assignment - Label the ASVs summarized here
• differential-abundance - Per-taxon between-group testing on the unrarefied counts
• qiime2-workflow - The QIIME2 CLI home for core-metrics and tree building
• phylogenetics/tree-io - Read, write, and root the UniFrac/Faith PD tree
• metagenomics/abundance-estimation - Shared CoDA and rarefaction-debate theory
• metagenomics/metagenome-visualization - Shared Hill-number and PERMANOVA-dispersion theory; diversity/ordination on shotgun profiler tables
• data-visualization/ggplot2-fundamentals - Custom ordination and diversity plots