Interpreting Community Composition

By the end of a metagenomics workflow you are holding several separate results: a taxonomic abundance table, alpha and beta diversity metrics, and a predicted or assembled functional profile. Interpretation means weaving these into one biological story rather than reporting three disconnected result sets. This lesson shows how to move from tables to a defensible narrative and where that narrative most often goes wrong.

Three Lenses One Story

Taxonomy tells you who is present, diversity tells you how varied and even the community is, and function tells you what metabolic capacity it might carry. A claim is strongest when the lenses agree. Lower alpha diversity in a disease group is far more convincing when it coincides with depletion of specific short-chain-fatty-acid producers and a measured drop in butyrate-synthesis gene abundance.

Bar Plots And Heatmaps

library(phyloseq)
library(ggplot2)

ps <- readRDS("metagenome_phyloseq.rds")
ps_rel <- transform_sample_counts(ps, function(x) x / sum(x))

# Stacked bar plot of phylum-level relative abundance, grouped by condition
ps_phylum <- tax_glom(ps_rel, taxrank = "Phylum")
plot_bar(ps_phylum, x = "Sample", fill = "Phylum") +
  facet_grid(~ group, scales = "free_x", space = "free_x") +
  labs(y = "Relative abundance")

# Clustered heatmap of the 50 most abundant genera
ps_genus <- tax_glom(ps_rel, taxrank = "Genus")
top50 <- names(sort(taxa_sums(ps_genus), decreasing = TRUE))[1:50]
plot_heatmap(prune_taxa(top50, ps_genus), taxa.label = "Genus", sample.label = "group")

Bar plots are for orientation and heatmaps for pattern discovery, not for statistics. Because each sample sums to one, the bars are compositional: an apparent rise in one taxon can simply reflect a real fall in another. Use the bar plot to spot dominant phyla and gross shifts, the clustered heatmap to reveal blocks of co-varying genera and sample groupings, then confirm any difference with a method built for compositional data.

Differential Abundance

library(ANCOMBC)

# Compositionally-aware differential abundance at genus level
res <- ancombc2(
  data        = ps,
  tax_level   = "Genus",
  fix_formula = "group",
  group       = "group",
  p_adj_method = "BH",
  prv_cut     = 0.10,
  struc_zero  = TRUE,
  alpha       = 0.05
)

# Effect size (log fold change), adjusted p-value, and significance flag
head(res$res[, c("taxon", "lfc_groupIBD", "q_groupIBD", "diff_groupIBD")])

             taxon lfc_groupIBD q_groupIBD diff_groupIBD
1      Bacteroides        -1.82    0.00061          TRUE
2 Faecalibacterium        -2.41    0.00013          TRUE
3      Escherichia         1.97    0.00420          TRUE
4        Roseburia        -1.55    0.02900          TRUE
5          Blautia        -0.34    0.31000         FALSE
6      Akkermansia         0.88    0.12000         FALSE

The lfc columns are natural-log fold changes relative to the reference group and diff flags the taxa whose Benjamini-Hochberg adjusted q-value passes the alpha threshold. Report the effect size and its direction, not just the flag: here Faecalibacterium is strongly depleted in IBD while Escherichia is enriched, consistent with the diversity and functional picture. Phrase this as an association. ANCOM-BC tells you a genus differs between groups; it does not tell you the genus caused the phenotype.