跳转至

微生物组差异丰度分析

一句话概述

微生物组差异丰度分析用于鉴定不同条件/分组间显著差异的微生物类群，核心方法包括ALDEx2、ANCOM-BC、MaAsLin2和LEfSe，各方法在组成性数据处理、假阳性控制和适用场景上各有优势。

核心知识点表格

知识点说明
组成性数据(Compositional data)微生物相对丰度之和为1，各特征间非独立，不能直接用标准统计方法
ALDEx2基于CLR变换+Monte Carlo模拟的差异分析，假阳性控制好
ANCOM-BC基于对数变换+偏差校正，估计采样分数(sampling fraction)差异
MaAsLin2多变量关联分析，支持混合效应模型和多种数据变换
LEfSeLinear Discriminant Analysis Effect Size，结合统计显著性和生物效应大小
CLR变换Centered Log-Ratio，处理组成性数据的标准方法
零膨胀(Zero-inflation)微生物数据中大量零值的问题
假阳性率(FDR)差异分析中需要严格控制的多重比较问题
效应量(Effect size)衡量差异的生物学意义（不仅看p-value）
混杂因素年龄、性别、BMI等协变量需要在模型中校正

各步骤详解

第一步：理解组成性数据问题

白话解释： 16S或宏基因组测序数据是"相对丰度"——每个样本的所有物种丰度加起来是100%。这意味着一个物种增加了，其他物种的相对丰度就会"被迫"下降，即使它们的绝对量没有变。这个"总和限制"会让传统统计方法（如t检验）产生虚假的差异结果。

技术细节： - 测序数据是从一个固定库容（如100万reads）中抽样 - Aitchison几何：组成性数据应在对数比空间中分析 - 常见变换：CLR（centered log-ratio）、ALR（additive log-ratio）、ILR（isometric log-ratio） - 零值处理：pseudocount、贝叶斯估计、multiplicative replacement

为什么不能直接用Wilcoxon/t-test：

样本A：物种X = 1000 reads (10%), 总reads = 10,000
样本B：物种X = 1000 reads (5%), 总reads = 20,000
结论：物种X的绝对量相同，但相对丰度不同！
直接比较相对丰度会产生spurious associations

第二步：ALDEx2分析

白话解释： ALDEx2先把count数据做CLR变换（对数中心化比），然后用蒙特卡洛模拟（Dirichlet分布）来估计数据的不确定性，最后在这些模拟实例上做统计检验。核心优势是假阳性控制最好。

技术细节： - 输入：count矩阵（行=特征，列=样本） - 使用Dirichlet-multinomial模型生成128/1000个实例 - 对每个实例做CLR变换 - 在CLR空间中做Welch's t-test或Wilcoxon检验 - 效应量用within/between组方差比衡量 - 两组和多组设计均支持

代码示例：

library(ALDEx2)

# 输入：count矩阵（行=OTU/species，列=样本）
# 注意：必须是整数count，不是相对丰度
count_matrix <- as.matrix(read.csv("otu_table.csv", row.names = 1))
# 分组信息
conditions <- c(rep("Control", 10), rep("Disease", 10))

# === 两组比较 ===
# 生成CLR实例（mc.samples=128是最低推荐，发表用1000）
aldex_clr <- aldex.clr(count_matrix, conditions, mc.samples = 128, denom = "all")

# 统计检验
aldex_ttest <- aldex.ttest(aldex_clr, paired.test = FALSE)

# 效应量
aldex_effect <- aldex.effect(aldex_clr)

# 合并结果
aldex_all <- data.frame(aldex_ttest, aldex_effect)

# 筛选显著差异特征
# 推荐使用 we.eBH（Welch's t BH校正后p值）< 0.05
# 且 |effect| > 1 表示有生物学意义的差异
sig_features <- aldex_all[aldex_all$we.eBH < 0.05 & abs(aldex_all$effect) > 1, ]
print(sig_features)

# 可视化
# MW plot（效应量 vs. 差异）
aldex.plot(aldex_all, type = "MW", test = "welch")

# MA plot
aldex.plot(aldex_all, type = "MA", test = "welch")

# === 多组比较（ANOVA-like） ===
conditions_multi <- c(rep("Healthy", 8), rep("Mild", 8), rep("Severe", 8))
aldex_clr_multi <- aldex.clr(count_matrix, conditions_multi, mc.samples = 128)
aldex_kw <- aldex.kw(aldex_clr_multi)

第三步：ANCOM-BC分析

白话解释： ANCOM-BC认为不同样本之间测序深度不同，导致一个"采样偏差"。它通过估计这个偏差并在统计模型中校正，使得结果近似于比较绝对丰度差异。支持多组比较和协变量校正。

技术细节： - 基于线性回归模型，加入采样分数(sampling fraction)偏差项 - ANCOM-BC2是最新版本，增强了对结构零值的处理 - 输出log fold change和校正后p值 - 内置多重比较校正（Holm, BH等） - 支持连续变量和分类变量混合模型 - 灵敏度和特异度平衡较好

代码示例：

library(ANCOMBC)
library(phyloseq)
library(TreeSummarizedExperiment)

# 输入：phyloseq对象
# 假设已有physeq对象包含otu_table + sample_data + tax_table
# physeq <- phyloseq(otu_table(count_matrix, taxa_are_rows=TRUE),
#                    sample_data(metadata),
#                    tax_table(taxonomy))

# === ANCOM-BC2（推荐新版） ===
ancombc2_result <- ancombc2(
  data = physeq,
  fix_formula = "Group + Age + Sex",  # 固定效应（校正协变量）
  rand_formula = NULL,                 # 随机效应（如有纵向数据）
  p_adj_method = "holm",
  pseudo_sens = TRUE,                  # 伪计数灵敏度分析
  prv_cut = 0.10,                      # 最低流行率阈值（去除太稀有的）
  lib_cut = 1000,                      # 最低文库大小
  s0_perc = 0.05,
  group = "Group",
  struc_zero = TRUE,                   # 结构零值检测
  neg_lb = TRUE,
  alpha = 0.05,
  n_cl = 4,                            # 并行核数
  verbose = TRUE,
  global = TRUE,                       # 全局检验（多组ANOVA-like）
  pairwise = TRUE,                     # 配对比较
  dunnet = TRUE,                       # Dunnett检验（与对照组比）
  trend = FALSE
)

# 提取结果
res_primary <- ancombc2_result$res

# 查看差异特征
# log fold change
head(res_primary[, c("taxon", "lfc_GroupDisease", "q_GroupDisease", "diff_GroupDisease")])

# 筛选显著
sig_taxa <- res_primary[res_primary$diff_GroupDisease == TRUE, ]
print(sig_taxa$taxon)

# 可视化
library(ggplot2)
plot_df <- res_primary %>%
  filter(diff_GroupDisease == TRUE) %>%
  arrange(lfc_GroupDisease)

ggplot(plot_df, aes(x = reorder(taxon, lfc_GroupDisease), y = lfc_GroupDisease)) +
  geom_bar(stat = "identity", aes(fill = lfc_GroupDisease > 0)) +
  coord_flip() +
  labs(x = "Taxon", y = "Log Fold Change (Disease vs Control)") +
  theme_minimal()

第四步：MaAsLin2分析

白话解释： MaAsLin2是一个灵活的多变量关联框架，可以同时考虑多个自变量（分组、年龄、BMI等），支持固定效应和随机效应（如纵向数据），可选择不同的数据变换和统计模型。

技术细节： - 支持变换：LOG、LOGIT、AST（arcsine square root）、NONE - 支持模型：LM（线性模型）、CPLM（复合泊松线性）、NEGBIN（负二项）、ZINB（零膨胀负二项） - 支持归一化：TSS（Total Sum Scaling）、CSS、TMM、CLR、NONE - 随机效应适合纵向/重复测量数据 - 输出：coefficient（效应大小）、p-value、q-value

代码示例：

library(Maaslin2)

# 输入准备
# features: 特征矩阵（行=样本，列=特征），TSV或data.frame
# metadata: 样本元数据（行=样本，列=变量）

# 注意MaAsLin2要求：行=样本，列=特征（与很多工具相反！）
features_df <- as.data.frame(t(count_matrix))  # 转置
metadata_df <- data.frame(
  SampleID = colnames(count_matrix),
  Group = c(rep("Control", 10), rep("Disease", 10)),
  Age = rnorm(20, 50, 10),
  Sex = sample(c("M", "F"), 20, replace = TRUE),
  row.names = colnames(count_matrix)
)

# 运行MaAsLin2
maaslin2_result <- Maaslin2(
  input_data = features_df,
  input_metadata = metadata_df,
  output = "maaslin2_output",
  fixed_effects = c("Group", "Age", "Sex"),  # 固定效应
  random_effects = NULL,                      # 随机效应（如Subject ID）
  reference = c("Group,Control"),             # 参考水平
  normalization = "TSS",                      # 归一化方法
  transform = "LOG",                          # 数据变换
  analysis_method = "LM",                     # 统计模型
  correction = "BH",                          # 多重比较校正
  min_abundance = 0.001,                      # 最低相对丰度
  min_prevalence = 0.1,                       # 最低流行率
  max_significance = 0.25,                    # q-value阈值
  cores = 4,
  plot_heatmap = TRUE,
  plot_scatter = TRUE
)

# 查看结果
sig_results <- maaslin2_result$results[maaslin2_result$results$qval < 0.25, ]
print(sig_results[, c("feature", "metadata", "coef", "pval", "qval")])

# 纵向数据示例
maaslin2_longitudinal <- Maaslin2(
  input_data = features_df,
  input_metadata = metadata_df,
  output = "maaslin2_longitudinal",
  fixed_effects = c("Timepoint", "Treatment"),
  random_effects = c("SubjectID"),  # 受试者随机效应
  normalization = "TSS",
  transform = "LOG",
  analysis_method = "LM"
)

第五步：LEfSe分析

白话解释： LEfSe先用Kruskal-Wallis检验找组间差异特征，再用Wilcoxon检验验证亚组一致性，最后用LDA（线性判别分析）估计效应量大小。LDA score > 2（或3、4）的特征被认为是biomarker。

技术细节： - 三步流程：KW检验 → Wilcoxon pairwise → LDA效应量 - 输入：相对丰度表 + 分组信息 - 非参数检验，不需要正态性假设 - LDA score反映该特征区分各组的能力 - 可生成cladogram（进化树形式展示差异物种） - 注意：LEfSe因不处理组成性问题而受批评，建议与其他方法交叉验证

代码示例：

# === Galaxy版本（在线）===
# 访问 https://huttenhower.sph.harvard.edu/galaxy/
# 上传TSV格式数据，运行LEfSe模块

# === 命令行版本 ===
# 安装
conda install -c bioconda lefse

# 输入格式准备（LEfSe特殊格式）
# 第1行：类别（Class）
# 第2行：亚类（Subclass，可选）
# 第3行：样本ID
# 第4行起：特征丰度
# 示例见下

# 1. 格式化输入
lefse_format_input.py input_table.tsv formatted_input.in \
  -c 1 \    # class行号
  -s 2 \    # subclass行号
  -u 3 \    # subject行号
  -o 1000000  # 归一化因子

# 2. 运行LEfSe
lefse_run.py formatted_input.in lefse_result.res \
  -a 0.05 \          # KW检验alpha
  -w 0.05 \          # Wilcoxon检验alpha
  -l 2.0 \           # LDA score阈值
  --min_c 10 \       # 最少样本数
  -e 1               # 多组LDA策略(1=one-vs-all)

# 3. 可视化
# 条形图
lefse_plot_res.py lefse_result.res lefse_barplot.pdf --format pdf

# 进化分支图（cladogram）
lefse_plot_cladogram.py lefse_result.res lefse_cladogram.pdf \
  --format pdf --dpi 300

# R中实现LEfSe（microbiomeMarker包）
library(microbiomeMarker)

# 从phyloseq对象运行
lefse_result <- run_lefse(
  physeq,
  group = "Group",
  subgroup = NULL,
  taxa_rank = "Genus",
  kw_cutoff = 0.05,
  wilcoxon_cutoff = 0.05,
  lda_cutoff = 2,
  norm = "CPM"           # Counts Per Million归一化
)

# 查看结果
marker_table(lefse_result)

# 可视化
plot_ef_bar(lefse_result)           # 效应量条形图
plot_cladogram(lefse_result)        # 进化分支图
plot_abundance(lefse_result, group = "Group")  # 丰度箱线图

第六步：方法比较与选择策略

白话解释： 没有一种方法是完美的。最佳实践是使用多种方法，取交集作为高置信度结果。

方法对比表：

特性ALDEx2ANCOM-BCMaAsLin2LEfSe
组成性校正CLR采样分数估计可选CLR无
假阳性控制最好好中等较差
灵敏度中等高高高
多变量支持有限支持最灵活不支持
纵向数据不支持支持支持不支持
效应量有（effect）有（LFC）有（coef）有（LDA）
计算速度慢(MC模拟)快快快
适用场景严格控制假阳性标准两/多组比较复杂实验设计快速探索性分析

代码示例（多方法取交集）：

# 取多方法交集
aldex_sig <- rownames(aldex_all[aldex_all$we.eBH < 0.05, ])
ancombc_sig <- res_primary$taxon[res_primary$diff_GroupDisease == TRUE]
maaslin_sig <- sig_results$feature[sig_results$metadata == "Group"]

# 两方法交集
consensus_2 <- Reduce(intersect, list(aldex_sig, ancombc_sig))

# 三方法交集（高置信度）
consensus_3 <- Reduce(intersect, list(aldex_sig, ancombc_sig, maaslin_sig))

# Venn图
library(VennDiagram)
venn.diagram(
  x = list(ALDEx2 = aldex_sig, ANCOMBC = ancombc_sig, MaAsLin2 = maaslin_sig),
  filename = "DA_methods_venn.png",
  fill = c("red", "blue", "green"),
  alpha = 0.5
)

实战命令（可复制）

# ===== 环境安装 =====
# R包安装
# install.packages("BiocManager")
# BiocManager::install(c("ALDEx2", "ANCOMBC", "phyloseq"))
# install.packages("Maaslin2", repos = "https://cran.r-project.org")
# BiocManager::install("microbiomeMarker")

# LEfSe安装
conda install -c bioconda -c conda-forge lefse

# ===== 数据准备 =====
# 从QIIME2导出count表
qiime tools export --input-path table.qza --output-path exported/
biom convert -i exported/feature-table.biom -o otu_table.tsv --to-tsv

# ===== R中完整流程 =====
# 见上文各方法的R代码

# ===== LEfSe命令行 =====
lefse_format_input.py input.tsv formatted.in -c 1 -u 2 -o 1000000
lefse_run.py formatted.in result.res -l 2.0 -a 0.05
lefse_plot_res.py result.res barplot.pdf --format pdf
lefse_plot_cladogram.py result.res cladogram.pdf --format pdf

面试常问点

Q1: 为什么微生物组数据不能直接用t检验/Wilcoxon检验？

A: 因为微生物组数据是组成性数据（compositional），各特征的相对丰度和恒为1。一个特征的增加必然导致其他特征相对减少（假性负相关），这违反了传统统计方法中特征独立性的假设。此外还有零膨胀、过度离散、测序深度不均等问题。需要使用对数比变换（如CLR）或专用方法。

Q2: ALDEx2的Monte Carlo模拟目的是什么？

A: 目的是估计count数据的技术不确定性。测序count值本质上是从一个多项式分布中抽样的结果，同一样本重测可能得到不同count。ALDEx2从Dirichlet分布生成多次模拟实例（默认128次），对每个实例分别做统计检验，最后取统计量的中位值或均值，这样结论对单次抽样波动更稳健。

Q3: ANCOM-BC的"BC"是什么意思？

A: BC = Bias Correction。ANCOM-BC估计不同样本间由于测序深度/效率差异导致的观测偏差（sampling fraction），并在线性模型中作为偏移量(offset)校正。这使得比较结果近似于比较绝对丰度的log fold change，而非仅仅比较相对丰度。

Q4: MaAsLin2在什么场景下是最佳选择？

A: (1) 需要校正多个协变量（年龄、性别、BMI、用药等）；(2) 纵向研究（重复测量）需要随机效应；(3) 连续变量与微生物丰度的关联分析；(4) 需要尝试不同变换和模型组合。MaAsLin2是最灵活的框架。

Q5: LEfSe目前还推荐使用吗？

A: LEfSe仍被广泛使用（引用量高），但有已知缺陷：(1) 不处理组成性问题；(2) 假阳性率较高；(3) 不支持协变量校正；(4) KW+Wilcoxon+LDA流程缺乏统一的统计框架。建议仅用于探索性分析和可视化（cladogram很直观），正式结论应用ALDEx2/ANCOM-BC验证。

Q6: 差异丰度分析需要多少样本量？

A: 取决于效应大小和期望检验力。经验值：(1) 每组至少10-15个样本才有合理统计力；(2) 如果预期差异大（如健康vs IBD），每组8-10个可能足够；(3) 效应小的差异（如饮食干预）需要20-50+/组。可用micropower包做功效分析。

Q7: 如何处理微生物组数据中的大量零值？

A: (1) ALDEx2：Dirichlet先验自然处理零（加入均匀先验后再生成实例）；(2) ANCOM-BC2：有structural zero检测模块区分"真零"和"采样零"；(3) MaAsLin2：ZINB模型专门处理零膨胀；(4) 通用：设置最低流行率阈值（如>10%样本中出现）过滤太稀有的特征。

Q8: 效应量和统计显著性哪个更重要？

A: 两者都需要关注。p-value/q-value告诉你差异是否"真实"（不是随机波动），效应量告诉你差异有多大（是否有生物学意义）。推荐同时报告：ALDEx2的effect>1、ANCOM-BC的|LFC|>1、MaAsLin2的|coef|配合q-value，避免报告统计显著但效应微小的特征。

易错点

1. 将相对丰度表当作count表输入ALDEx2

问题： ALDEx2需要整数count矩阵，输入相对丰度（0-1之间的小数）会导致Dirichlet采样出错。 正确做法： 使用原始read count或rarefied count，不做任何归一化。

2. 未过滤低丰度/低流行率特征

问题： 太稀有的特征（只在1-2个样本中出现）检验力极低且增加多重比较负担。 正确做法： 过滤：流行率<10-20%的特征移除；最低count/相对丰度阈值。

3. ANCOM-BC中未正确指定参考水平

问题： 多组比较时不指定参考组，结果解读困难。 正确做法： 在formula或reference参数中明确指定对照组（如Control）。

4. MaAsLin2的输入矩阵方向搞反

问题： MaAsLin2要求行=样本，列=特征。很多OTU表是行=特征、列=样本。 正确做法： 用t()转置OTU表后再输入。

5. LEfSe的输入格式不正确

问题： LEfSe有特殊的TSV格式要求（前几行是元数据），格式错误会静默产生错误结果。 正确做法： 严格按照格式规范准备输入或使用R包microbiomeMarker从phyloseq对象直接运行。

6. 多重比较校正后无显著结果就放弃

问题： 小样本量+大量特征导致FDR校正后无显著结果。 正确做法： (1) 减少检验数量：先过滤稀有特征或提高分类学层级（属→科）；(2) 使用更灵敏的方法；(3) 报告未校正p-value但明确标注；(4) 增加样本量。

7. 忽略生物学重复和技术重复的区分

问题： 把技术重复（同一样本重复测序）当作独立样本，人为增加样本量。 正确做法： 技术重复应合并（sum counts），只有生物学重复才算独立样本。

补充知识

新兴方法

方法特点
LinDA与ANCOM-BC类似但更简洁的线性模型
ZicoSeq基于排列检验的零膨胀组成性方法
corncobBeta-binomial模型，直接建模count和离散度
DESeq2/edgeRRNA-seq差异分析工具，可用于微生物组但需谨慎（不处理组成性）
songbird/Qurro多项式回归建模差异排名

可视化最佳实践

  1. 火山图：x轴=效应量/LFC，y轴=-log10(q-value)
  2. 热图：显著特征在各样本中的丰度
  3. 箱线图/小提琴图：top差异特征在各组的分布
  4. Cladogram：LEfSe特色，展示差异物种的系统发育关系
  5. Venn图：多方法结果取交集

发表论文建议

  • 至少使用两种方法交叉验证
  • 明确报告所用参数（阈值、归一化方法等）
  • 报告效应量而非仅报p-value
  • 对协变量进行校正
  • 描述零值处理策略
  • 代码和数据公开（可复现性）

本教程系统对比了四种主流微生物组差异丰度分析方法，建议根据实验设计选择合适方法并进行多方法交叉验证。

📚 相关文章推荐

191_微生物组常见统计错误
36% 相关
181_微生物组功能通路差异
36% 相关
微生物组关联分析MWAS（菌群-表型关联分析）
35% 相关
LEfSe差异分析详解
34% 相关
MaAsLin2多变量关联
33% 相关