Bio Metagenomics Visualization

Version Compatibility

Reference examples tested with: phyloseq 1.46+, vegan 2.6+, microViz 0.12+, ALDEx2 1.34+, pandas 2.2+, scikit-bio 0.6+, matplotlib 3.8+.

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

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• Python: pip show <package> then help(module.function) to check signatures
• CLI: ktImportTaxonomy (no args) to confirm Krona column flags on the installed build

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

The input profiler decides the toolchain: MetaPhlAn gives relative abundance (cannot rarefy; ships an SGB tree so UniFrac is available); Bracken gives counts (rarefaction and count models valid; no tree); HUMAnN gives gene-family/pathway features. Record the profiler, its filtering (Kraken --confidence, Bracken threshold, MetaPhlAn --stat_q), and every modeling choice below - they are not recoverable from the figure.

Metagenome Visualization

"Show me how my communities differ." -> Choose a transform, a distance, and a test - each a modeling choice that can create or erase the difference - then declare them and show the conclusion survives them.

• R: phyloseq + vegan + microViz on a parsed profiler table
• Python: pandas + scikit-bio + matplotlib (plotting and wrangling; R is primary for community stats)

Scope: visualizing and testing a shotgun profiler table. Input generation -> kraken-classification, metaphlan-profiling, abundance-estimation, functional-profiling. Compositional theory and absolute load -> abundance-estimation. Amplicon/QIIME2 stats -> the microbiome category. Generic plotting primitives -> data-visualization.

The Single Most Important Modern Insight -- A Figure Is a Modeling Choice, Not an Observation

The distance metric, the transform, the rarefaction depth, the confidence filter, and the top-N cutoff are knobs turned before the answer is seen; turning them differently gives a different paper. A shotgun table is a compositional, depth-confounded, false-positive-laden estimate, and the figure inherits all of it. The job is to declare the modeling choices and show they did not fabricate the conclusion. Memory hooks:

• A stacked bar of relative abundance hides absolute load and visually inflates whatever is already dominant - it shows the relative race, not whether the community bloomed or collapsed.
• Bray-Curtis on relative abundance is the field default and is compositionally incoherent; Aitchison (Euclidean on CLR) is correct and almost nobody runs it.
• Shotgun richness is a readout of the database and the confidence threshold, not biology.
• A significant PERMANOVA may be a difference in variability, not composition.
• The differentially abundant taxa reported depend more on which DA tool was run than on biology.

The fork that decides everything downstream: MetaPhlAn = percent, Bracken = counts. Rarefaction is only meaningful for counts; CLR needs a pseudocount on percentages; richness estimators (Chao1/Observed) require integer counts and are meaningless on MetaPhlAn relative abundances.

Honest Composition Plots

Krona gives the full drillable hierarchy - often the honest answer to "what is in it," versus a bar that pre-collapses to the top 10:

kreport2krona.py -r kraken.kreport -o krona.txt && ktImportText krona.txt -o krona.html
# Confirm ktImportTaxonomy column flags with no-arg usage; they drift across versions.

For a stacked bar: collapse rare taxa to "Other" but label how many taxa and what percent that hides; state n per group (never stack one representative sample); use a colorblind-safe palette of <=12 colors; show absolute load alongside if total-biomass data exist. microViz comp_barplot() handles Other and palettes.

Alpha Diversity: Hill Numbers and Richness Honesty

Goal: Report within-sample diversity in interpretable units without letting shotgun richness masquerade as biology.

Approach: Use Hill numbers (effective species) at q=0/1/2; prefer evenness-weighted q=1/q=2 for shotgun because richness (q=0) is dominated by database size and false positives; only compute richness estimators on integer counts.

library(phyloseq); library(vegan)
# estimate_richness Observed/Chao1 assume INTEGER COUNTS - valid for Bracken, meaningless on MetaPhlAn %.
alpha <- estimate_richness(ps_counts, measures = c('Shannon', 'InvSimpson'))
hill_q1 <- exp(alpha$Shannon)        # effective species (q=1), interpretable units
hill_q2 <- alpha$InvSimpson          # q=2, dominance-weighted; robust to rare-taxon noise

Richness moves drastically with the classifier's confidence/threshold and 25-70% of shotgun species can be false positives - report the filtering, show a rarefaction or coverage curve, and gate richness behind unique-k-mer evidence (KrakenUniq). The rarefaction debate is unresolved: McMurdie & Holmes 2014 (PLoS Comput Biol 10:e1003531) call rarefying inadmissible for differential abundance; Schloss 2024 (mSphere 9:e00355-23) defends it for diversity. Decide per analysis - rarefy/coverage-standardize for diversity, model-based for DA - and show the curve either way.

Beta Diversity and Ordination

PCoA decomposes any distance; PCA on CLR is the compositional-coherent ordination and gives taxon loadings (which taxa drive an axis) that PCoA cannot; NMDS reports stress (not variance); UMAP does not preserve global distances and is for spotting clusters, not measuring dissimilarity. Pick the geometry for a stated reason and show the conclusion survives at least one alternative.

library(vegan)
dist_bc <- vegdist(otu_matrix, method = 'bray')   # samples as rows
pm <- adonis2(dist_bc ~ group, permutations = 999, by = 'terms')
# ALWAYS pair PERMANOVA with a dispersion test - a significant adonis2 can be a spread difference,
# not a location shift (Anderson & Walsh 2013), especially for unbalanced designs.
bd <- betadisper(dist_bc, group); permutest(bd, permutations = 999)

If betadisper is significant, the PERMANOVA is ambiguous - report both.

Differential Abundance: Consensus, Not a Single Tool

DA methods disagree wildly across datasets (Nearing 2022 Nat Commun 13:342); the taxa called significant depend more on the tool than on biology. Run at least two compositionally aware methods, report the intersect as high-confidence and the union as exploratory, and name every tool. ALDEx2 and ANCOM-II were the most conservative/consistent in Nearing; LinDA and ANCOM-BC were well FDR-controlled in Yang & Chen 2022. Prevalence-filter first and BH-correct across taxa.

The anti-pattern to forbid: an uncorrected Wilcoxon or t-test on raw relative abundances - wrong on compositionality and on multiple testing in one line. Most tools want counts (Bracken); for MetaPhlAn percentages use methods that accept proportions (MaAsLin2/LinDA) or convert to pseudo-counts.

Compositional Ordination in Python

Goal: Produce a compositionally coherent ordination in Python instead of the incoherent StandardScaler-then-PCA on raw relative abundance.

Approach: Replace zeros, CLR-transform, then PCA on the CLR coordinates (this is Aitchison-PCA); the loadings are interpretable as taxa.

import pandas as pd
from skbio.stats.composition import clr, multi_replace   # renamed from multiplicative_replacement in skbio 0.6
from sklearn.decomposition import PCA

ab = pd.read_csv('merged_abundance.txt', sep='\t', index_col=0)
ab = ab[ab.index.str.contains(r'\|s__') & ~ab.index.str.contains(r'\|t__')]  # species rows only
proportions = (ab.T.values / ab.T.values.sum(axis=1, keepdims=True))
clr_mat = clr(multi_replace(proportions))                # zeros replaced, then CLR (NOT StandardScaler on relab)
pca = PCA(n_components=2).fit(clr_mat)
coords = pca.transform(clr_mat)