# Bio Applied Microbial Diversity > Compute alpha/beta diversity (Shannon, Simpson, Bray-Curtis, UniFrac) from a 16S/ASV feature table with scikit-bio; PCoA ordination, PERMANOVA/ANOSIM. Use for microbiome diversity or community composition questions. - Skill: `pavel-kravchenko/bio-applied-microbial-diversity` (Agent Skill) - Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-microbial-diversity` - Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-microbial-diversity/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/bio-applied-microbial-diversity --- # Microbial Diversity Analysis ## When to Use - Computing within-sample diversity (Shannon, Simpson, observed richness, Faith's PD) from an OTU/ASV feature table. - Comparing community composition between groups (body site, treatment, timepoint) with Bray-Curtis or UniFrac distances. - Visualizing sample relationships via PCoA/NMDS ordination or hierarchical clustering. - Testing whether groups differ significantly in composition (PERMANOVA, ANOSIM) or whether sequencing depth is sufficient (rarefaction curves). - NOT for the QIIME2 CLI pipeline itself (raw reads -> DADA2 -> feature table) — see `bio-applied-qiime2-16s` for that; this skill covers the downstream diversity math once you already have a feature table. ## Version Compatibility scikit-bio >= 0.6 (API used here verified against 0.7.0), pandas >= 1.5, numpy >= 1.24, scikit-learn >= 1.3 (for NMDS via `sklearn.manifold.MDS`), Python >= 3.9. ## Prerequisites - `pip install scikit-bio pandas numpy scipy scikit-learn matplotlib` - A feature table (rows = samples, columns = OTUs/ASVs, values = read counts) and sample metadata with a grouping column. scikit-bio's diversity functions expect **samples as rows** — transpose a taxa-by-sample table (`table.T`) first if needed. - Concepts: relative abundance, rarefaction, distance matrices. Related: `bio-applied-qiime2-16s` (upstream pipeline), `bio-applied-taxonomic-profiling` (shotgun metagenomics). ## Alpha Diversity: Within-Sample Richness and Evenness **Goal:** quantify how many taxa are present and how evenly they are distributed, per sample. **Approach:** call `skbio.diversity.alpha_diversity` with a named metric over the counts matrix; compare groups with Kruskal-Wallis since diversity indices are rarely normally distributed. ```python import numpy as np import pandas as pd from scipy import stats from skbio.diversity import alpha_diversity def alpha_diversity_by_group(feature_table: pd.DataFrame, metadata: pd.DataFrame, group_col: str) -> pd.DataFrame: """Compute Shannon, Simpson, and observed-richness alpha diversity per sample and test for group differences with Kruskal-Wallis. feature_table: samples (rows) x OTUs/ASVs (columns), integer counts. metadata: sample-indexed DataFrame containing `group_col`. Returns a DataFrame of per-sample metrics with the group label attached. """ counts = feature_table.values.astype(int) ids = feature_table.index.tolist() metrics = { "shannon": alpha_diversity("shannon", counts, ids=ids), "simpson": alpha_diversity("simpson", counts, ids=ids), "sobs": alpha_diversity("sobs", counts, ids=ids), # observed richness (formerly 'observed_otus') } result = pd.DataFrame(metrics) result[group_col] = metadata.loc[result.index, group_col] for name in metrics: groups = [result.loc[result[group_col] == g, name] for g in result[group_col].unique()] stat, p = stats.kruskal(*groups) print(f"{name}: Kruskal-Wallis H={stat:.2f}, p={p:.4f}") return result ``` For phylogenetic richness, use `alpha_diversity("faith_pd", counts, ids=ids, otu_ids=feature_table.columns.tolist(), tree=tree)` where `tree` is a `skbio.TreeNode` with tip names matching the OTU/ASV ids. ## Beta Diversity and Ordination: Between-Sample Comparison **Goal:** measure and visualize how different samples' communities are from one another. **Approach:** build a distance matrix with `skbio.diversity.beta_diversity`, reduce it to 2D with PCoA (classical MDS on the distance matrix) for a quick, reproducible plot. ```python import pandas as pd from skbio.diversity import beta_diversity from skbio.stats.ordination import pcoa def ordinate_samples(feature_table: pd.DataFrame, metric: str = "braycurtis"): """Build a beta-diversity distance matrix and run PCoA on it. metric: any skbio beta-diversity metric, e.g. 'braycurtis', 'jaccard', or 'unweighted_unifrac'/'weighted_unifrac' (the latter two require otu_ids= and tree= kwargs pointing to a phylogenetic tree). Returns (distance_matrix, ordination_results) where ordination_results.samples holds PC1/PC2/... coordinates and .proportion_explained the variance per axis. """ counts = feature_table.values.astype(int) ids = feature_table.index.tolist() dm = beta_diversity(metric, counts, ids=ids) ord_results = pcoa(dm) explained = ord_results.proportion_explained print(f"PC1 explains {explained.iloc[0] * 100:.1f}%, PC2 {explained.iloc[1] * 100:.1f}%") return dm, ord_results ``` ## Statistical Testing: Do Groups Differ? **Goal:** decide whether community composition significantly differs between groups (e.g. body sites), not just "looks different" on an ordination plot. **Approach:** run PERMANOVA (and/or ANOSIM) directly on the distance matrix with permutation-based p-values. ```python import pandas as pd from skbio.stats.distance import permanova, anosim def test_group_differences(distance_matrix, metadata: pd.DataFrame, group_col: str, n_permutations: int = 999) -> dict: """Test whether groups differ in community composition using PERMANOVA and ANOSIM. distance_matrix: a skbio.DistanceMatrix (e.g. from beta_diversity()). metadata: sample-indexed DataFrame; must be reduced/aligned to the distance matrix's sample ids and the grouping Series must be named `group_col` (skbio looks up the column by that name). """ grouping = metadata.loc[list(distance_matrix.ids), group_col] grouping.name = group_col # skbio requires a named Series/column perm = permanova(distance_matrix, grouping, permutations=n_permutations) ano = anosim(distance_matrix, grouping, permutations=n_permutations) return { "permanova_F": perm["test statistic"], "permanova_p": perm["p-value"], "anosim_R": ano["test statistic"], "anosim_p": ano["p-value"], } ``` ## Pitfalls - **Rarefaction depth**: unequal library sizes bias alpha diversity comparisons (deeper samples look artificially "richer"); rarefy to the minimum library size or use depth-invariant metrics before comparing across samples. - **Skewed distributions**: diversity indices are rarely normal — use Kruskal-Wallis/Mann-Whitney, not t-tests/ANOVA, unless you've checked normality. - **Multiple testing**: apply FDR correction (Benjamini-Hochberg) when testing thousands of taxa for differential abundance. - **UniFrac needs a real tree**: `unweighted_unifrac`/`weighted_unifrac` require a `skbio.TreeNode` whose tip names exactly match your OTU/ASV ids — a mismatched or missing tree silently breaks the metric. - **PERMANOVA is sensitive to dispersion, not just location**: a significant PERMANOVA can reflect different within-group variance (heterogeneous dispersion) rather than different centroids — pair it with `skbio.stats.distance.permdisp` (Betadisper equivalent) to check. - **Batch effects**: always check for batch confounding (extraction kit, sequencing run) before interpreting biological signal. ## See Also - `bio-applied-qiime2-16s` — upstream pipeline: raw reads to feature table (DADA2, taxonomy). - `bio-applied-taxonomic-profiling` — shotgun metagenomics taxonomic profiling (Kraken2/Bracken). - `bio-applied-statistics-for-bioinformatics` — general permutation testing and multiple-testing correction. - `bio-applied-population-genetics` — related distance/ordination methods (Fst, PCA) for population-level data.