Bio Microbiome Differential Abundance

Version Compatibility

Reference examples tested with: ALDEx2 1.34+, ANCOMBC 2.4+, Maaslin2 1.16+, MicrobiomeStat 1.2+ (LinDA), GUniFrac 1.8+ (ZicoSeq), phyloseq 1.46+.

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

• R: packageVersion('<pkg>') then ?function_name to verify parameters

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

ANCOM-BC2 changed argument names between ancombc() and ancombc2(), and its default p_adj_method is holm, not BH - confirm both against the installed version. The MaAsLin3 maaslin3() API differs from MaAsLin2's Maaslin2().

Differential Abundance Testing

"Find which taxa differ between my groups" -> Run two or more compositionally-aware DA tools and report their consensus - because the significant-taxa list is a property of the tool as much as of the sample, and a relative-abundance change is not an absolute change.

• R: ALDEx2::aldex(counts, conds, test='t', effect=TRUE, denom='all') then a second tool (ANCOMBC::ancombc2() or MicrobiomeStat::linda())

Scope: per-taxon DA on an amplicon feature table. Whole-community alpha/beta/PERMANOVA -> diversity-analysis. Shotgun profiler-table DA -> metagenomics/metagenome-visualization. Shared compositional/closure/CLR/zero theory -> metagenomics/abundance-estimation. Collapse ASVs to genus/species first -> taxonomy-assignment. QIIME2 CLI route -> qiime2-workflow.

The Single Most Important Modern Insight -- Which Taxa Are "Significant" Depends More on the Tool Than on the Biology

Run ALDEx2, ANCOM-BC2, MaAsLin2, and LinDA on the same ASV table and the four significant-taxa lists overlap but disagree (Nearing 2022 Nat Commun 13:342, across 38 datasets). So the deliverable is NOT "the differential taxa" - it is the CONSENSUS of >=2 compositionally-aware tools, every tool NAMED: the intersection is high-confidence, the union is exploratory, and a single-tool hit is tentative. Picking the tool with the prettiest volcano is p-hacking by software (uncorrected multiplicity hidden in the method menu). Three corollaries:

1. A relative-abundance increase is not an absolute increase. Microbiome counts are compositional - the sequencer fixes the total, so one taxon blooming forces every other taxon's proportion down (the blooming-taxon illusion). "Taxon X increased" is a statement about its SHARE unless an external load anchor (spike-in / flow cytometry / qPCR, see metagenomics/abundance-estimation) or MaAsLin3's absolute-abundance mode licenses an absolute claim.
2. Uncorrected Wilcoxon/t-test on raw relative abundances is wrong twice in one line - closure (reference-frame) AND multiple testing. But a BH-corrected simple test, honestly labelled as relative, can replicate BETTER than a sophisticated model (Pelto 2025): the forbidden thing is the uncorrected, closure-blind form, not simple tests per se.
3. There is no settled best tool. The benchmarks optimize different criteria, so they rank tools differently. Consensus-of-tools is the only stance that survives all of them.

The Benchmark Landscape (no settled winner)

Report the disagreement AS the result; verify current best practice against the latest tool docs rather than hard-coding one method.

Tool Taxonomy

Decision Tree by Scenario

Filter Before Testing (a modeling knob, not housekeeping)

Goal: Drop rare features before testing so the BH denominator is not crushed and log/CLR transforms are well-behaved.

Approach: Keep features present in at least 10-25% of samples (and optionally a mean-abundance floor); declare the threshold and confirm the headline result is not knife-edge-sensitive to it. Every tool exposes this (prv_cut, min_prevalence, prev.filter).

library(phyloseq)
ps <- readRDS('phyloseq_object.rds')
# prv_cut 0.10: a feature must appear in >= 10% of samples; raising to 0.25 removes more tests
# (smaller BH correction, more power on survivors) but discards rare-but-real taxa - a declared choice
keep <- filter_taxa(ps, function(x) sum(x > 0) >= 0.10 * nsamples(ps), TRUE)

ALDEx2: The Conservative Floor of the Consensus

Goal: Identify taxa that differ between two groups while propagating the sampling uncertainty of low-count features.

Approach: Draw mc.samples Monte-Carlo instances from a Dirichlet posterior of the counts (this IS the zero handling - no explicit pseudocount), CLR-transform each instance against the geometric mean of all features (denom='all'), run the test on every draw, and report the EXPECTED effect size and BH-adjusted p over the draws.

library(ALDEx2)
counts <- as.matrix(otu_table(ps))            # integer counts, taxa in ROWS
if (!taxa_are_rows(ps)) counts <- t(counts)
groups <- as.character(sample_data(ps)$Group)

# mc.samples 128: standard Monte-Carlo draws; 256+ for publication (more stable expected p)
res <- aldex(counts, groups, mc.samples = 128, test = 't', effect = TRUE, denom = 'all')
# we.eBH = Welch expected BH-adjusted p (report this, NOT we.ep); wi.eBH = Wilcoxon equivalent
# effect = median standardized effect = median(diff.btw / max(diff.win)); the primary decision variable
hits <- res[res$we.eBH < 0.05 & abs(res$effect) > 1, ]   # q AND effect floor (Gloor: gate on effect, not p alone)

Gate on effect size AND q, not p alone: with large n trivially small CLR differences become "significant," and Gloor's own guidance is that |effect| > 1 is a strong ~2-SD signal. For >2 groups use aldex.kw(); for covariates the aldex.glm() + model.matrix route works but ALDEx2 is weakest here - prefer ANCOM-BC2/LinDA/MaAsLin2 for serious covariate or random-effect modeling.

ANCOM-BC2: Bias-Corrected LFC With a Sensitivity Safeguard

Goal: Estimate an interpretable bias-corrected log-fold-change per taxon, with covariate adjustment, structural-zero handling, and a flag for hits that are hostage to the pseudo-count.

Approach: Model log(observed count) as a function of covariates, estimate each sample's log sampling fraction as an offset and subtract it, then refit across a range of pseudo-counts and record how often each q-value flips (passed_ss).

library(ANCOMBC)
out <- ancombc2(data = ps, fix_formula = 'Group + Age + Sex',
                rand_formula = NULL,        # '(1 | SubjectID)' for repeated measures - see Failure Modes
                p_adj_method = 'BH',        # DEFAULT is 'holm'; set 'BH' deliberately for FDR
                prv_cut = 0.10, lib_cut = 1000,
                group = 'Group', struc_zero = TRUE, pseudo_sens = TRUE,
                global = FALSE, pairwise = FALSE, n_cl = 2)
res <- out$res
# a confident hit is BOTH significant AND robust to the pseudo-count. ANCOM-BC2 suffixes the
# diff_/passed_ss_ columns with the literal model-matrix coefficient (variable + factor level,
# verbatim case, e.g. 'Grouptreated') - match it by pattern rather than hard-coding the case.
dcol <- grep('^diff_Group', names(res), value = TRUE)[1]
robust <- res[res[[dcol]] & res[[sub('^diff_', 'passed_ss_', dcol)]], ]

passed_ss is the most valuable ANCOM-BC2-specific feature: a CLR/log model on sparse data is hostage to the zero-replacement constant, and passed_ss quantifies that per taxon. A hit with passed_ss == FALSE depends on the arbitrary pseudo-count - do not report it as confident. For >2 groups set global=TRUE (omnibus), pairwise=TRUE (mdFDR-controlled pairs), dunnet=TRUE, or trend=TRUE; results land in out$res_global/res_pair/res_dunn/res_trend.

LinDA: Fast CLR Regression With Native Mixed Models

Goal: Get FDR-controlled log2-fold-changes on a large cohort, including repeated-measures designs, without Monte-Carlo or EM cost.

Approach: Fit ordinary linear regression on the CLR-transformed table covariate by covariate, estimate the compositional bias as the mode of the per-feature coefficients and subtract it; a random effect in the formula makes it a linear mixed model.

library(MicrobiomeStat)
otu <- as.data.frame(otu_table(ps)); if (!taxa_are_rows(ps)) otu <- t(otu)
meta <- as.data.frame(sample_data(ps))
fit <- linda(feature.dat = otu, meta.dat = meta,
             formula = '~ Group + Age + (1 | SubjectID)',   # random effect -> mixed model
             feature.dat.type = 'count', prev.filter = 0.10, alpha = 0.05)
fit$output[[1]]   # names(fit$output) are the model-matrix coefficient columns (e.g. 'Grouptreated' - the factor level keeps its case); per-feature: log2FoldChange, lfcSE, stat, pvalue, padj, reject

LinDA is the natural fast modern entry in a consensus panel and the cleanest route to mixed models. Yang & Chen rate it among the best FDR-power trade-offs.

MaAsLin2 / MaAsLin3 and ZicoSeq (the rest of the panel)

Goal: Fit covariate-rich or longitudinal differential-abundance models, or add a permutation-based panel member, when ALDEx2/ANCOM-BC2/LinDA do not cover the design.

Approach: Use MaAsLin2/3 for multivariable GLMs with random effects, or ZicoSeq for a non-parametric permutation-FDR test against empirically selected reference taxa.

MaAsLin2 fits a flexible per-feature GLM; its package DEFAULT is TSS + LOG + LM (not CLR), and random_effects is the canonical route for longitudinal designs. NOTE the orientation gotcha: it expects features in COLUMNS, samples in rows.

library(Maaslin2)
fit <- Maaslin2(input_data = as.data.frame(t(otu)), input_metadata = meta,
                output = 'maaslin2_out', fixed_effects = c('Group', 'Age'),
                random_effects = c('SubjectID'),
                normalization = 'TSS', transform = 'LOG', analysis_method = 'LM',
                min_prevalence = 0.10, max_significance = 0.05)
# writes all_results.tsv / significant_results.tsv with columns feature, metadata, coef, pval, qval

MaAsLin3 (maaslin3()) splits each feature into an abundance model (level when present) and a logistic prevalence model (present/absent) tested jointly, and can ingest total-load measurements for absolute-abundance inference. ZicoSeq (GUniFrac::ZicoSeq()) winsorizes, posterior-samples, normalizes against empirically selected reference taxa, and returns permutation FDR (zc$p.adj.fdr) - a non-parametric panel member that accepts covariates via adj.name.

Consensus: Intersect the Tools

Goal: Convert two or more per-tool hit sets into a confidence-graded result instead of one tool's answer.

Approach: Collect the significant feature SETS (BH within each tool), then report the intersection as high-confidence, the union as exploratory, and tabulate, per taxon, how many of N tools agree and which ones. Never pool p-values across tools.

sig_aldex <- rownames(res)[res$we.eBH < 0.05 & abs(res$effect) > 1]
sig_linda <- rownames(fit$output[[1]])[fit$output[[1]]$reject]   # [[1]] = the group coefficient (named 'Grouptreated')
confident  <- intersect(sig_aldex, sig_linda)   # high-confidence
exploratory <- union(sig_aldex, sig_linda)       # report with the tool that found each

Per-Method Failure Modes

Cherry-picking the tool with the prettiest result

Trigger: running several tools and reporting only the one(s) that flag the favored taxon. Mechanism: that is uncorrected multiplicity hidden in the method menu (p-hacking by software). Symptom: "the recommended method found X" with no mention of the tools that disagreed. Fix: decide the panel a priori, report ALL tools, intersect for confident hits, disclose disagreement.

Uncorrected Wilcoxon/t-test on relative abundances

Trigger: a per-taxon Wilcoxon/t-test on TSS proportions with no FDR correction. Mechanism: closure makes a naive test call every taxon "decreased" when one blooms, and hundreds of uncorrected tests inflate false positives. Symptom: dozens of "significant" taxa, all in the same direction, no q-values. Fix: use a CoDA/reference-frame tool; if a simple test is used, BH-correct it and label the comparison as relative (Pelto 2025).

Pseudo-replication of repeated measures

Trigger: longitudinal/paired samples treated as independent rows. Mechanism: fewer independent units than rows inflates significance. Symptom: implausibly small p-values on a small subject count. Fix: a random effect - ANCOM-BC2 rand_formula='(1|SubjectID)', MaAsLin2 random_effects='SubjectID', LinDA (1|SubjectID) in the formula. Cross-check ANCOM-BC2 mixed-model output against LinDA/MaAsLin2 (GitHub issue #111 reported rand_formula correctness problems in some versions).

Prevalence filter set blindly

Trigger: an undeclared prevalence cut, or none at all. Mechanism: the cut decides which taxa are even tested and thus the BH landscape - it is a modeling choice. Symptom: the hit list changes materially between prv_cut=0.1 and 0.25. Fix: declare and justify the threshold; confirm the headline result survives moving it.

Relative change reported as absolute

Trigger: "taxon X doubled" from a closed table with no load data. Mechanism: one taxon blooming compresses every other proportion. Symptom: whole-community "depletion" that is really one taxon rising. Fix: anchor to load (spike-in/flow/qPCR) or MaAsLin3 absolute mode; otherwise state the claim is relative.

DESeq2/edgeR on a sparse 16S table

Trigger: RNA-seq median-of-ratios / TMM on a zero-heavy ASV table. Mechanism: the geometric-mean size-factor reference collapses on zeros and the "most features unchanged" assumption is violated. Symptom: degenerate size factors, errors, or inflated hit counts that disagree with CoDA tools (Nearing). Fix: use a compositional tool; if DESeq2 is unavoidable, the poscounts estimator is the minimum mitigation - present as a caveat, not a recipe.

ANCOM-BC2 hit held hostage by the pseudo-count

Trigger: reporting diff_* == TRUE without checking passed_ss_*. Mechanism: significance depends on the arbitrary zero-replacement constant. Symptom: a hit that vanishes when the pseudo-count changes. Fix: require diff_* & passed_ss_* for a confident call.

Quantitative Thresholds

Common Errors

References

• Fernandes AD, Reid JNS, Macklaim JM, McMurrough TA, Edgell DR, Gloor GB. 2014. Unifying the analysis of high-throughput sequencing datasets: characterizing RNA-seq, 16S rRNA gene sequencing and selective growth experiments by compositional data analysis. Microbiome 2:15.
• Gloor GB, Macklaim JM, Fernandes AD. 2016. Displaying variation in large datasets: plotting a visual summary of effect sizes. J Comput Graph Stat 25:971-979.
• Lin H, Peddada SD. 2020. Analysis of compositions of microbiomes with bias correction (ANCOM-BC). Nat Commun 11:3514.
• Lin H, Peddada SD. 2024. Multigroup analysis of compositions of microbiomes with covariate adjustments and repeated measures (ANCOM-BC2). Nat Methods 21:83-91.
• Mallick H, Rahnavard A, McIver LJ, et al. 2021. Multivariable association discovery in population-scale meta-omics studies. PLoS Comput Biol 17:e1009442.
• Nickols WA, Kuntz T, Shen J, et al. 2026. MaAsLin 3: refining and extending generalized multivariable linear models for meta-omic association discovery. Nat Methods 23:554-564.
• Zhou H, He K, Chen J, Zhang X. 2022. LinDA: linear models for differential abundance analysis of microbiome compositional data. Genome Biol 23:95.
• Yang L, Chen J. 2022. A comprehensive evaluation of microbial differential abundance analysis methods: current status and potential solutions. Microbiome 10:130.
• Yang L, Chen J. 2023. Benchmarking differential abundance analysis methods for correlated microbiome sequencing data. Brief Bioinform 24:bbac607.
• Pelto J, Auranen K, Kujala JV, Lahti L. 2025. Elementary methods provide more replicable results in microbial differential abundance analysis. Brief Bioinform 26(2):bbaf130.
• Nearing JT, Douglas GM, Hayes MG, et al. 2022. Microbiome differential abundance methods produce different results across 38 datasets. Nat Commun 13:342.
• Segata N, Izard J, Waldron L, et al. 2011. Metagenomic biomarker discovery and explanation. Genome Biol 12:R60.
• Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15:550.
• Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B 57:289-300.

Related Skills

• diversity-analysis - Whole-community alpha/beta/PERMANOVA; answer "do the communities differ" before "which taxa differ"
• taxonomy-assignment - Collapse ASVs to genus/species before per-taxon testing
• amplicon-processing - Produces the ASV feature table tested here
• qiime2-workflow - The qiime composition ancombc CLI route
• metagenomics/abundance-estimation - Shared compositional/closure/CLR/zero/load-anchor theory
• metagenomics/metagenome-visualization - The same DA mechanics on shotgun profiler tables
• experimental-design/multiple-testing - FDR control and multiplicity across taxa