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.

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.

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.

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.