微生物差异丰度分析(LEfSe、ANCOM-BC)
微生物差异丰度分析用于找出在两组或多组样本之间丰度显著不同的菌属、物种或功能通路,例如疾病组中富集的致病菌、干预后减少的产丁酸菌。它比转录组差异分析更棘手:测序得到的是组成性数据,每个样本的相对丰度之和固定,一种菌大量增加会让其他菌看起来减少;同时零值多、测序深度差异大。常用方法包括 LEfSe(非参数检验加 LDA 效应量)、ANCOM 与 ANCOM-BC(基于对数比并校正样本间偏倚)、ALDEx2(CLR 转换)和 MaAsLin2(多变量线性模型)。输入是 16S 或宏基因组的计数表或丰度表加样本元数据。不同方法的结果往往差别很大,稳健的结论通常需要多种方法一致支持。
什么时候用
- 寻找疾病组与对照组之间富集或缺失的菌属或物种
- 调整年龄、BMI、用药等协变量后检验分组与菌群丰度的关联
- 比较干预前后同一个体内菌群的变化
- 筛选与血糖、炎症因子等连续指标相关的类群
- 比较宏基因组功能通路在组间的丰度差异
不适合的情况
- 想检验群落整体组成是否不同,应先做 β 多样性与 PERMANOVA
- 每组只有两三个样本,却想得到稳健的差异菌列表
- 把相对丰度的变化直接解释为绝对菌量的增减
需要准备的数据
- 计数表
- ASV、属或物种水平的原始读段计数;ANCOM-BC、ALDEx2 需要计数,MaAsLin2 也可用相对丰度。
- 样本元数据
- 分组和协变量,连续变量单位统一;纵向或配对数据需有个体 ID 和时间点。
- 分类信息
- 每个特征从门到属或种的注释,用于按分类水平聚合和在图中标注。
- 预过滤标准
- 事先确定流行率和丰度过滤阈值并写入方法,避免反复调整阈值筛出想要的结果。
分析步骤
- 1
确定分析水平
决定在 ASV、属还是种水平分析。ASV 数量多、零值多,检验功效低;属水平聚合后更稳定,但会掩盖同属不同种的反向变化。主要水平应事先确定。
- 2
特征预过滤
去除在大多数样本中为零的稀有特征,例如只保留在至少 10% 样本中出现的类群(经验值)。过滤能减少多重检验负担,但阈值要在看结果之前设定。
- 3
选择主方法
两组比较且无协变量时可用 LEfSe 做探索;需要调整协变量时用 ANCOM-BC2 或 MaAsLin2;重视组成性偏差时用 ANCOM-BC 或 ALDEx2。主方法定好后再选一到两种方法验证。
- 4
运行与校正
对每个特征得到效应量和 P 值,并做 FDR 校正。逐个核对显著类群的效应方向是否与原始相对丰度箱线图一致,方向相反时检查参考组设置。
- 5
多方法一致性
汇总各方法的显著结果,用韦恩图或交集表展示。多种方法都显著且方向一致的类群作为主要结论,只有一种方法支持的类群作为探索性结果单独列出。
- 6
结果可视化
绘制 LDA 柱状图、分类分支图、效应量森林图或火山图,并对关键类群绘制分组箱线图展示原始分布,让读者能看到效应背后的数据。
关键参数与决策
LDA 阈值
LEfSe 默认 LDA 对数得分阈值为 2.0,常见设为 2 到 4。阈值越高结果越少、越保守,但 LDA 值是效应量而不是检验统计量,不能替代多重校正后的 P 值。
显著性与多重校正
LEfSe 默认不做 FDR 校正,Kruskal-Wallis 的 α 默认为 0.05;MaAsLin2 默认以 q 值小于 0.25 判定显著,较为宽松,正式结论常改用 0.05 或 0.1。
结构零处理
ANCOM-BC2 可以识别在某组中完全不存在的结构零类群,单独列出而不进入对数比检验。这类菌更适合描述为「仅在某组中出现」,而不是给出倍数变化。
伪计数
对数转换前需给零值加伪计数,伪计数大小会明显改变低丰度菌的效应量。ANCOM-BC2 提供针对伪计数的敏感性分析,可据此标记不稳健的结果。
协变量与随机效应
抗生素使用、年龄、测序批次是常见混杂因素。纵向或配对数据在 MaAsLin2 或 ANCOM-BC2 中把个体 ID 设为随机效应,忽略配对会高估显著性。
常见错误
- 只用 LEfSe 且不说明未做多重校正:类群数量多时假阳性比例高,结论的可靠性会受到质疑。
- 直接用 DESeq2 分析零值极多的菌群计数:部分评估显示其在菌群数据中假阳性偏高,需与组成性方法交叉验证。
- 在门、科、属多个水平分别检验并分别校正:上下级结果重复计数,检验总数被低估,应明确主要分析水平。
- 把相对丰度差异解释为绝对增减:其他菌大量增加也会让某菌相对丰度下降,需要 qPCR 或流式计数等绝对定量支持。
- 多种方法结果不一致时只挑显著的报告:造成选择性报告偏倚,应同时给出各方法的结果数量和交集。
结果怎么看
LDA 得分柱状图
每根柱子是一个差异类群,长度为 LDA 对数得分,颜色表示富集的组。柱子长只说明效应大,不代表统计证据更强;同一分类链上的父子节点同时出现时,实际是同一个信号。
分类分支图
按分类层级展示差异类群在分类树中的位置。差异集中在某个门或科的分支下,提示该类群整体发生变化;差异零散分布时,结果更可能受随机波动影响。
效应量结果表
列出每个类群的对数倍数变化、标准误、P 值和校正后 q 值。效应量很大但 q 值不显著,常因组内变异大或流行率低,需结合原始数据箱线图判断。
方法一致性图
用韦恩图或热图展示各方法识别出的差异类群交集。交集中的类群是相对稳健的结论;交集为空时,说明效应较弱或样本量不足,应谨慎下结论。
和相近方法怎么选
LEfSe 与 ANCOM-BC
LEfSe 基于非参数检验和 LDA 效应量,速度快、图形直观,但不处理组成性、默认不做多重校正,也难以纳入协变量;ANCOM-BC 估计并校正样本间抽样偏倚,给出带标准误的对数倍数变化,更适合正式推断。探索用 LEfSe,结论用 ANCOM-BC。
MaAsLin2 与 ALDEx2
MaAsLin2 用线性模型或混合模型,可同时纳入多个协变量和随机效应,适合复杂队列;ALDEx2 通过 CLR 转换和蒙特卡洛抽样估计不确定性,对组成性处理较好,但对复杂设计支持有限。
菌群差异丰度与转录组差异表达
两者都是逐个特征检验,但菌群数据零值更多、组成性更强。转录组标准化通常假设多数基因不变,这在菌群中不一定成立,直接套用 DESeq2、edgeR 可能产生偏差,需要结合组成性方法。
常见问题
LEfSe 的 LDA 值多少算显著?
LDA 值本身不是显著性指标,而是效应量。LEfSe 先用 Kruskal-Wallis 检验筛出 P 小于阈值的特征,再计算 LDA 值,默认阈值 2.0,论文中常用 2 或 3。方法部分应同时写明检验阈值和 LDA 阈值。
ANCOM 和 ANCOM-BC 有什么区别?
ANCOM 通过大量对数比检验得到 W 统计量,只判断是否差异,不给出效应量和 P 值。ANCOM-BC 在线性模型框架中估计并校正样本特异的抽样偏倚,输出对数倍数变化、标准误和 P 值;ANCOM-BC2 进一步支持多组比较、重复测量和敏感性分析。
不同方法结果不一致怎么办?
这很常见。建议预先指定一个主方法,再用一到两种原理不同的方法验证,把多种方法一致的类群作为主要发现,其余作为探索性结果。同时报告各方法的显著类群数量,便于读者判断结论的稳健程度。
差异丰度分析前要不要抽平?
通常不需要。抽平会随机丢弃读段、降低检验功效,大多数差异丰度方法自带标准化或偏倚校正,例如 ANCOM-BC 估计样本特异的抽样偏倚,ALDEx2 使用 CLR 转换。抽平主要用于让 α 多样性指数在相同深度下可比。
宏基因组功能通路也能用这些方法吗?
可以。MaAsLin2 常用于 HUMAnN 输出的通路丰度,ANCOM-BC 也可用于计数形式的功能表。功能表同样是组成性数据,且不同通路共享大量基因、彼此高度相关,解释时不要把相关通路当成相互独立的发现。
在 Codar 里怎么做
在 Codar 中提供计数表、分类注释和样本元数据后,会先确认分析的分类水平、预过滤阈值、主方法以及需要调整的协变量,再进行微生物差异丰度分析:运行主方法得到效应量和校正后 P 值,并用原理不同的方法交叉验证,汇总一致的结果。交付 LDA 柱状图或效应量图、结果表、方法一致性图、关键类群箱线图、统计报告和可复现脚本。
我有 80 例结直肠癌和健康对照的属水平计数表,想在调整年龄和 BMI 后找差异菌,以 ANCOM-BC 为主、LEfSe 作补充。