做微生物组研究的人大概率都绕不开 LEfSe 这四个字母。不管是16S rRNA测序还是宏基因组测序在“找差异物种”这一步几乎都会看到 LEfSe 的三件套结果一张五颜六色的LDA条形图、一圈一圈往外扩的Cladogram分支图以及若干张差异物种的丰度分布图。LEfSe全称 Linear discriminant analysis Effect Size直译是线性判别分析效应大小最早由哈佛公共卫生学院的Huttenhower团队在2011年发布。它解决的问题非常具体在样本量不大、特征维度却高达几千上万且大多数物种丰度都是零的微生物数据里怎么找出真正能区分两组或多组样本的物种。这篇文章我就从分析原理、数据准备、命令行实操到结果解读和常见报错把LEfSe整个流程讲透适合正在做16S、ITS、宏基因组差异分析的同学参考。1. 项目概述LEfSe 到底是什么为什么大家都在用1.1 一个不算新但一直能打的差异分析工具LEfSe 的定位很明确它是一套“先筛选、再计算效应量、最后可视化”的生物标志物发现流程主要用于发现不同分组之间丰度存在显著差异的微生物类群。在微生物组论文里LEfSe的结果几乎等于“差异物种”的代名词很多期刊审稿人也默认你要给出LDA score和Cladogram图。从技术实现上看LEfSe并不是一个单一的统计方法而是一套流水线内部串联了Kruskal-Wallis秩和检验、Wilcoxon秩和检验和线性判别分析。它接受一个高维的物种丰度表格输出哪些物种在哪个分组里显著富集以及富集程度有多大。整个过程对数据分布没有严格要求也不需要样品量达到几十个因此在微生物组这类“特征多、样本少、零值多”的数据上特别适用。为什么它十多年了还在被大量使用我个人的看法是它的结果非常直观LDA条形图一眼就能看出哪些物种在哪个组更占优势Cladogram又能把不同分类层级的富集情况展示在一张图里很适合放在论文结果部分。1.2 LEfSe 和普通差异分析有什么不一样很多人一开始会跑DESeq2、edgeR或者直接做Wilcoxon检验然后把显著性p值筛出来就算完事。LEfSe跟这些方法最大的区别在于它不只看p值还要算一个“效应量”。p值只能说明“这个物种在组间有没有统计显著差异”但没法回答“这个差异到底有多大、方向是什么”。LEfSe通过LDA计算出每个物种对区分不同组别的贡献大小用LDA score来表示这样就不容易陷入“p值很小但实际差异微弱”的尴尬。另外LEfSe自带了一套“阶梯式过滤”逻辑先用非参数检验排除大部分无关特征再通过两两比较和LDA筛选进一步压缩候选物种。相比之下单独做几十上百次Wilcoxon检验往往需要自己做多重校正而且结果很零散还要自己画图。LEfSe把这些环节整合起来输出格式相对规范。不过也要实话实说LEfSe不是唯一的选择ANCOM、ALDEx2、ANCOM-BC这些工具在某些场景下可能更严谨。但从实用主义的角度看LEfSe依然是目前“性价比”最高的一套流程尤其适合做探索性分析和论文插图。1.3 什么场景适合用 LEfSe16S/ITS扩增子测序得到的OTU/ASV丰度表或者宏基因组物种注释丰度表。样本有明确分组比如疾病组和对照组、处理组和空白组、不同时间点或不同环境梯度。需要比较两组或多组并且想找出具体富集的物种。想要输出LDA条形图、Cladogram图这类论文常用图。不适合的场景也很明显没有分组的样本、样本量每组只有1到2个、数据质量特别差且分组混杂因素巨大这些情况下LEfSe结果不一定可靠。结论要结合生物学背景不能只靠一个工具拍板。2. 原理不搞懂结果容易误读2.1 第一关Kruskal-Wallis 检验在整个分析里的作用LEfSe的第一步是Kruskal-Wallis检验这是一个典型的非参数检验方法主要用来判断多个组之间某个特征的中位数是否存在显著差异。这里为什么不直接用ANOVA呢因为微生物丰度数据极少满足正态分布和方差齐性大量物种在多数样本里都是0数据分布是极度偏态的。Kruskal-Wallis检验的思路是把所有样本的特征值混合后按大小排序用秩次代替原始数值做比较。它不太受极端值影响适合处理带有大量0和少数高丰度值的数据。LEfSe默认要求这个检验的p值小于0.05才认为这个物种在多个组之间的分布有差异才有资格进入下一轮筛选。实际分析时要注意这一步筛选出来的候选特征可能仍然非常多尤其当样本分组差异很大时可能会有数百个物种通过检验。这很正常后续两步就是进一步压缩这些候选。2.2 第二关Wilcoxon 两两比较解决“谁和谁不同”Kruskal-Wallis检验只回答“几个组之间是否存在差异”但它不会告诉你是A组和B组不同还是A组和C组不同。所以LEfSe第二步会进行Wilcoxon秩和检验做两两组间比较。默认情况下LEfSe会对通过第一步的特征做Wilcoxon检验任意两组之间的差异都需要达到显著水平这个特征才会被保留。对于多组比较还提供了one-against-all或者更严格的策略具体可以在运行参数里控制。这一步有点像是“组间差异一致性检查”避免某些特征只是被单一极端分组拉高了整体p值。我遇到不少同学只记住了LEfSe会跑Wilcoxon却忽略了两两比较的策略差异。如果分组超过三组建议在跑之前想清楚你是想找某一组相对于其他所有组的特异marker还是想找任意两组之间都能分开的marker这个选择会直接影响最终结果。2.3 第三关LDA score 不是 p 值是效应量通过了前两轮检验的特征接下来会进入线性判别分析阶段。线性判别分析的基本思想是寻找一个线性组合让组间差异尽可能大、组内差异尽可能小。在LEfSe里每个特征都会被计算出一个判别系数这个系数经过标准化后就是LDA score通常取log10后的值默认以2作为显著阈值。很多人第一次看LDA条形图会以为横轴大小就是丰度差异大小其实不是。LDA score反映的是这个特征在区分不同组别时的重要性/贡献度含义更接近“效应量”。所以假设某物种在两组中丰度差异明显但组内波动也很大它的LDA score未必很高。反过来一个丰度不是最高、但组间规律性很强的物种LDA score可能会很高。这个特点决定了LEfSe找出来的物种通常具备“组间差异一致性好”的特征这也是为什么审稿人比较认可这套结果。2.4 阈值用多少alpha 和 LDA score 怎么调LEfSe里有两个重要参数统计检验的显著性阈值alpha默认0.05和LDA score阈值默认2.0。前者控制Kruskal-Wallis和Wilcoxon检验的p值标准后者控制最终展示多少个marker。如果样本量不是很大比如每组只有五六例跑出来一个显著marker都没有可以尝试把LDA阈值降到2.0以下比如1.5但不能降得太狠否则结果会变得很“虚”。反过来如果差异物种特别多图纸一片混乱可以把LDA阈值提高到3.0甚至3.5减少图形上的信息量只保留最有区分度的物种。关于多重检验校正LEfSe默认流程其实没有做特别严格的FDR校正所以它的p值只能作为筛选参考不建议当作确定性结论。如果论文需要更严格的统计支撑最好再用ANCOM-BC或ALDEx2做交叉验证。3. 分析前的数据准备3.1 特征表和分组信息长什么样LEfSe分析需要两个核心输入一个是特征丰度表另一个是样本分组信息表。特征丰度表一般长这样第一行是样本名第一列是特征名中间是丰度值Taxon Sample_CTR1 Sample_CTR2 Sample_CTR3 Sample_DZ1 Sample_DZ2 Sample_DZ3 k__Bacteria|p__Firmicutes 120 150 110 80 90 70 k__Bacteria|p__Bacteroidota 200 180 210 320 340 310 k__Bacteria|p__Proteobacteria 50 40 45 160 150 170分组信息表也必须是tab分隔第一列是样本名第二列是分组。如果有二级分组可以在第三列写上subclassSample_CTR1 Control Batch1 Sample_CTR2 Control Batch1 Sample_CTR3 Control Batch2 Sample_DZ1 Disease Batch1 Sample_DZ2 Disease Batch1 Sample_DZ3 Disease Batch2这里必须强调分组信息中的样本名要和特征表的列名完全一致顺序不要求一致但名字不能有多余空格不要加引号。3.2 输入格式常见雷区LEfSe对输入格式比较挑剔最常见的坑有三个。第一缺失值处理不当。特征表里出现空单元格、NA或者NaN后面跑format_input.py时很容易直接报错。处理方式很简单所有缺失值都填成0。第二空格和特殊字符。Excel里做表格后直接复制到txt经常会把tab变成空格或者让样本名多一个看不见的空格。建议在Linux下用head -5和cat -A检查一下分隔符如果列间是^I才是tab。另外特征名里不要带空格、括号、引号全部用下划线替代。第三分类层级前缀不完整。如果后面想画Cladogram图特征名必须按照k__、p__、c__、o__、f__、g__、s__这种层级格式并用|连接。如果只是输入一个OTU编号或ASV编号没有分类注释LEfSe也能跑但Cladogram画不出来或者画出来只有一个孤零零的点没法看。3.3 要不要先做预处理LEfSe本身没有做丰度归一化它直接使用你输入的特征值进行计算。因此如果你的输入是原始测序得到的reads count建议提前转成相对丰度或者用总丰度缩放的方式归一化。很多流程会输出相对丰度表直接拿来跑就可以。另一个重要操作是过滤低丰度特征。微生物组数据可能有上万行特征但绝大多数物种在少数样本里出现一次这种特征进了LEfSe只会拖慢速度对结果贡献也很小。我一般会先过滤一遍在所有样本中相对丰度都低于0.01%的特征直接去掉如果有明显是叶绿体、线粒体、未分类的序列也在这一步剔除。值得注意的是过滤不能太激进。如果样本之间差异比较细微你把低丰度物种全部删掉可能把真正有生物学意义的稀有物种也丢了。建议同时保留原始未过滤表方便后续对比。3.4 在线版还是本地版LEfSe最早提供在线Galaxy版用户上传两个文件、选好参数就能跑适合快速出图。在线版的优点是省去安装环境对于不熟悉命令行的同学非常友好。缺点是大数据量上传比较费时间而且每次都要手动操作重复性不够好。本地版则需要安装LEfSe程序用命令行运行好处是脚本化、可复现批量处理多个数据集的时候效率高。我个人更推荐本地版尤其是当你需要跑几十个样本、要多次调整阈值或者要在服务器上做批量分析的时候。本地版也没有样本量限制只要机器内存扛得住。如果只是两三个小数据集在线版也能满足需求。总之建议两条路都试一下本地跑出来结果不理想时再去在线平台核对一遍这样心里更有底。4. 实操从命令行跑通 LEfSe4.1 安装LEfSe本地安装时建议用conda创建独立环境避免和项目里其他Python环境冲突。操作步骤如下conda create -n lefse python3.9 -y conda activate lefse conda install -c bioconda lefse -y也可以尝试用pip安装pip install lefse安装完成后验证一下脚本是否可用run_lefse.py --help format_input.py --help如果能看到帮助信息说明安装成功。需要提醒的是LEfSe不同小版本之间参数命名略有差异实际使用前先用--help看一眼以你本机版本为准。老教程里的Python 2命令在新环境中不一定能直接跑。4.2 准备两个数据文件我在服务器上习惯建一个工作目录把所有文件放进去比如LEfSe_test/。目录下放两个文件feature_table.txt特征丰度表第一行样本名第一列特征名。group_info.txt分组信息表第一列样本名第二列分组。写完后检查文件格式head -5 feature_table.txt head -5 group_info.txt这一步很关键很多报错都是因为分隔符不是tab。4.3 format_input.py把原始表格转换成 LEfSe 内部格式正式分析前需要先用format_input.py将特征表和分组信息合并成一个LEfSe内部格式文件。基本命令如下python format_input.py feature_table.txt group_info.txt lefse_input.txt -c 1 -u 2其中-c 1表示分组信息文件里class所在的列是第1列-u 2表示可选的subclass在第2列。如果没有subclass就不要加-u。这一步执行完会生成lefse_input.txt。可以用head看一下内容文件格式已经变成LEfSe能识别的结构。如果这里报错说明输入文件格式有问题需要回到上一步检查分隔符和缺失值。4.4 run_lefse.py核心统计检验生成内部格式文件之后就可以跑LEfSe核心分析python run_lefse.py lefse_input.txt lefse_result.txt -l 2 -a 0.05 -w 0.05参数含义-l 2LDA score阈值默认2。-a 0.05Kruskal-Wallis检验的显著性阈值。-w 0.05Wilcoxon检验的显著性阈值。-o 1多组比较时使用one-against-all策略可选。运行结束后会生成lefse_result.txt里面每一行是一个通过筛选的特征包含它对应的分组、LDA score、p值等信息。这个文件是后续画图的输入也是你自己做数据整理时的重要依据。我觉得第一次跑的时候最好把过程输出的日志保留万一结果异常翻日志能看出是在哪一步筛选掉的。4.5 画图LDA条形图和Cladogram图LEfSe自带几个画图脚本最常用的是python plot_res.py lefse_result.txt lda_bar.pdf python plot_cladogram.py lefse_result.txt cladogram.pdfplot_res.py生成的是LDA条形图展示每个显著物种的LDA scoreplot_cladogram.py生成的是Cladogram分支图展示差异物种的分类层级关系。如果需要查看每个显著特征的丰度分布还可以用python plot_features.py lefse_result.txt feature_plots/这会在feature_plots/目录下生成每一个显著特征在不同分组中的丰度分布图适合用来快速核查结果是否合理。画图时也可调整参数比如plot_res.py可以用--format pdf指定输出格式或者通过-f调整图片宽度。实际需求不同可以自己探索。4.6 一个可复现的完整示例这里给一个从零开始的完整命令序列假设已在LEfSe_test/目录下准备好两个输入文件# 第一步合并数据 python format_input.py feature_table.txt group_info.txt lefse_input.txt -c 1 # 第二步跑LEfSeLDA阈值2.0 python run_lefse.py lefse_input.txt lefse_result.txt -l 2 -a 0.05 -w 0.05 # 第三步画图 python plot_res.py lefse_result.txt lda_bar.pdf python plot_cladogram.py lefse_result.txt cladogram.pdf # 第四步查看结果表格 head -20 lefse_result.txt这样一个流程跑下来几十个样本的数据通常几分钟就能完成。如果特征表特别大耗时可能稍长但一般不会超过半个小时。5. 结果怎么看、怎么写进文章5.1 LDA 条形图先看效应量再看颜色LDA条形图是LEfSe最经典的结果。图中每个条形代表一个显著差异物种横轴是LDA score条形越长说明该物种对分组的判别贡献越大。条形颜色对应富集分组红色、绿色、蓝色通常分别代表不同组。看这张图我一般分三步。第一步看整体数量显著物种是不是太多或太少如果太多可以调高LDA阈值再跑一次第二步看颜色分配是否合理各组的marker是否有生物学逻辑第三步挑几个LDA score最高的物种回到相对丰度表里验证一下它们在组间的差异是否真实。颜色相同的条形聚集在一起通常说明这些物种属于同一个高的分类单元协同变化。5.2 Cladogram 分支图从门到属的富集路线Cladogram图看起来像一圈圈的年轮内部是门往外依次是纲、目、科、属、种。图上用不同颜色标记在不同组中显著富集的物种节点大小代表丰度大小圆圈越靠近外圈分类越具体。读懂Cladogram的关键是理解“由内到外”的分类层级。比如某个节点的分类学注释是p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Lactobacillaceae|g__Lactobacillus在图上会显示为一串从内到外的节点如果整条路径都被标记为红色说明从门到属这个完整分支在红色组里都有显著富集。论文里常说的“Firmicutes门及其下属分类在疾病组显著富集”就是根据这张图得出的。需要提醒的是Cladogram图很好看但信息密度高阅读时容易眼花。汇报时最好把关键分类层级单独列出来配合LDA条形图一起说明。5.3 特征丰度图和热图除了两张主图LEfSe还可以输出每个显著特征在分组间的丰度分布图这就是plot_features.py的作用。这类图能直观展示某个物种在不同样本中的高低情况适合放到补充材料。如果想做热图可以把lefse_result.txt里显著的物种提取出来再结合原始丰度表用R的pheatmap画。一般做法是先选出LDA score排名靠前的物种然后按样本做z-score标准化再画热图。这样能把所有显著marker和样本之间的关系一次性展示出来。5.4 论文里怎么描述 LEfSe 结果写论文时LEfSe方法部分可以写成这样“Linear discriminant analysis Effect Size (LEfSe) was applied to identify taxa with significant abundance differences between groups, with the threshold of LDA score set to 2.0 and the significance level set to 0.05.”结果部分则避免只说“某些物种有差异”最好结合具体效应量。比如“LEfSe analysis identified 15 taxa significantly enriched in the disease group, including g_Lactobacillus (LDA score 4.2, P 0.003) and f_Enterobacteriaceae (LDA score 3.8, P 0.01).”我见过不少同学在结果部分放一张LDA图和一张Cladogram图但正文只字不解释这样读者很难抓住重点。图的展示需要文字引导至少要说明图里最大的特征是哪些物种。6. 常见问题与排查记录6.1 报错与解决办法速查表常见报错或异常可能原因处理办法format_input.py 报 “input file” 错误分隔符不是tab或存在空缺值用cat -A检查文件补0或重导出tab分隔文本样本名对不上特征表和分组信息中样本名不一致用 cut -f1 group_info.txtCladogram 画不出来特征名缺少分类层级前缀或 连接结果文件为空Kruskal-Wallis 或 Wilcoxon 无显著特征降低 alpha或检查分组是否合理LDA图条形数量太多差异物种过多提高-l阈值如3.0或3.5运行时间太长特征过多、样本过多先过滤低丰度特征再跑LEfSe结果和普通t检验不一致LEfSe采用多步筛选更保守结合其他差异分析方法交叉验证6.2 分析结果“没有显著 marker”这是最常遇到的问题之一。有时候分组差异确实存在但LEfSe一个显著结果都没跑出来。原因可能有三个第一样本量太少。每组只有两个样本统计检验几乎没有效能很难通过Kruskal-Wallis检验。这种情况可以尝试把LDA阈值降低到1或者接受当前数据不适合LEfSe的事实改用效果更直接的描述性分析。第二组间差异太弱。微生物群落整体差异不大个别物种的丰度变化也不稳定LEfSe的两次秩和检验会把大部分特征过滤掉。可以考虑换用ANCOM-BC等更灵敏的方法或者增加样本量。第三阈值设得太严。默认LDA大于2确实是一个比较高的标准如果数据集噪声大可以先用默认参数跑一遍再降到1.5左右看看趋势。但不能无限降否则就失去LEfSe的意义了。6.3 分组样本量不均或特别小时怎么处理LEfSe对每组最少样本数没有强制规定但实际经验是每组少于3个就会很不可靠少于5个时结果波动很大。如果你手里就是小样本建议同时做以下操作先看各组样本分布是否重叠降低LDA阈值后观察marker是否仍然一致再用随机森林等方法做特征重要性分析看是否和LEfSe结果有交集。如果有的组样本特别多、有的组特别少比如疾病组20个、对照组4个跑出来的差异物种往往会偏向样本量大的那一组。这种情况下可以对大组做随机抽稀比如从20个样本里随机抽4个重复多次看哪些物种在大样本和小样本组合情况下都能稳定出现。这个做法比较耗时间但能明显提高结论可信度。6.4 结果与之前做的差异分析不一致不少人在跑LEfSe之前已经用过DESeq2或者Wilcoxon做了差异物种筛选结果发现两套方法给出的物种列表对不上。这很正常因为LEfSe更关注“效应量”步骤也更严格。DESeq2这类方法针对的是原始计数矩阵基于负二项分布假设LEfSe直接处理相对丰度或标准化数据不依赖这个假设。两种方法得到交集之外的结果不代表谁错了只能说明它们对“差异物种”的定义不同。写论文时可以把两套结果并排展示取交集作为核心结论把只有一种方法检出的物种放到补充材料。这样做比只放一套结果更有说服力。7. 个人心得与避坑经验7.1 LEfSe 不是万能的别把结果当真理我见过一种比较危险的做法把LEfSe跑出来的所有显著物种直接拿去做机制解释甚至当成“因果证据”。实际上LEfSe只能说明在给定分组条件下这些物种的丰度分布存在差异不能说明是因果关系、也不能排除混杂因素。比如疾病组和对照组年龄不匹配跑出来的很多物种可能只是反映年龄差异而不是疾病本身。因此每次跑完LEfSe之后我都会先做一个“反常识检查”这些显著富集的物种是否有文献支持如果一个物种在以往多个研究中都被认为和疾病相关那这次结果就比较可信如果完全不符合已有认知先别急着高兴很可能数据里存在某个未控制的混杂因素。7.2 实际运行中的几个小习惯我自己的操作习惯是所有原始数据先做标准化和过滤得到相对丰度表后再做一次分组检查确认每组样本数量、样本名是否匹配。然后先用默认参数跑一遍记录结果再微调LDA阈值跑第二遍比较两次结果的重合情况。如果两次差异很大说明数据本身不太稳定需要小心。此外lefse_result.txt这个结果文件一定要保留。很多后续操作都要用到它比如提取显著物种做热图、做韦恩图、或者给其他工具做输入。不要只留PDF图丢了中间结果后面会很被动。7.3 最后一个实用小技巧如果你想在论文里把LDA条形图画得更精致可以直接读取lefse_result.txt用R的ggplot2重新绘制。LEfSe自带的图虽然方便但配色和排版不一定符合期刊要求。自己画图时把物种名按照LDA score降序排列用不同颜色代表分组再调整字体和坐标轴效果会好很多。另外多组比较时LEfSe的Cladogram图可能非常拥挤节点标签重叠严重。这时候与其硬调图片尺寸不如拆开画按每个显著的门或者纲单独画一个子图读者看起来更清楚审稿人也更舒服。希望这些实测经验能帮你少走一点弯路。
LEfSe分析全指南：原理、实操、图表解读与避坑经验
做微生物组研究的人大概率都绕不开 LEfSe 这四个字母。不管是16S rRNA测序还是宏基因组测序在“找差异物种”这一步几乎都会看到 LEfSe 的三件套结果一张五颜六色的LDA条形图、一圈一圈往外扩的Cladogram分支图以及若干张差异物种的丰度分布图。LEfSe全称 Linear discriminant analysis Effect Size直译是线性判别分析效应大小最早由哈佛公共卫生学院的Huttenhower团队在2011年发布。它解决的问题非常具体在样本量不大、特征维度却高达几千上万且大多数物种丰度都是零的微生物数据里怎么找出真正能区分两组或多组样本的物种。这篇文章我就从分析原理、数据准备、命令行实操到结果解读和常见报错把LEfSe整个流程讲透适合正在做16S、ITS、宏基因组差异分析的同学参考。1. 项目概述LEfSe 到底是什么为什么大家都在用1.1 一个不算新但一直能打的差异分析工具LEfSe 的定位很明确它是一套“先筛选、再计算效应量、最后可视化”的生物标志物发现流程主要用于发现不同分组之间丰度存在显著差异的微生物类群。在微生物组论文里LEfSe的结果几乎等于“差异物种”的代名词很多期刊审稿人也默认你要给出LDA score和Cladogram图。从技术实现上看LEfSe并不是一个单一的统计方法而是一套流水线内部串联了Kruskal-Wallis秩和检验、Wilcoxon秩和检验和线性判别分析。它接受一个高维的物种丰度表格输出哪些物种在哪个分组里显著富集以及富集程度有多大。整个过程对数据分布没有严格要求也不需要样品量达到几十个因此在微生物组这类“特征多、样本少、零值多”的数据上特别适用。为什么它十多年了还在被大量使用我个人的看法是它的结果非常直观LDA条形图一眼就能看出哪些物种在哪个组更占优势Cladogram又能把不同分类层级的富集情况展示在一张图里很适合放在论文结果部分。1.2 LEfSe 和普通差异分析有什么不一样很多人一开始会跑DESeq2、edgeR或者直接做Wilcoxon检验然后把显著性p值筛出来就算完事。LEfSe跟这些方法最大的区别在于它不只看p值还要算一个“效应量”。p值只能说明“这个物种在组间有没有统计显著差异”但没法回答“这个差异到底有多大、方向是什么”。LEfSe通过LDA计算出每个物种对区分不同组别的贡献大小用LDA score来表示这样就不容易陷入“p值很小但实际差异微弱”的尴尬。另外LEfSe自带了一套“阶梯式过滤”逻辑先用非参数检验排除大部分无关特征再通过两两比较和LDA筛选进一步压缩候选物种。相比之下单独做几十上百次Wilcoxon检验往往需要自己做多重校正而且结果很零散还要自己画图。LEfSe把这些环节整合起来输出格式相对规范。不过也要实话实说LEfSe不是唯一的选择ANCOM、ALDEx2、ANCOM-BC这些工具在某些场景下可能更严谨。但从实用主义的角度看LEfSe依然是目前“性价比”最高的一套流程尤其适合做探索性分析和论文插图。1.3 什么场景适合用 LEfSe16S/ITS扩增子测序得到的OTU/ASV丰度表或者宏基因组物种注释丰度表。样本有明确分组比如疾病组和对照组、处理组和空白组、不同时间点或不同环境梯度。需要比较两组或多组并且想找出具体富集的物种。想要输出LDA条形图、Cladogram图这类论文常用图。不适合的场景也很明显没有分组的样本、样本量每组只有1到2个、数据质量特别差且分组混杂因素巨大这些情况下LEfSe结果不一定可靠。结论要结合生物学背景不能只靠一个工具拍板。2. 原理不搞懂结果容易误读2.1 第一关Kruskal-Wallis 检验在整个分析里的作用LEfSe的第一步是Kruskal-Wallis检验这是一个典型的非参数检验方法主要用来判断多个组之间某个特征的中位数是否存在显著差异。这里为什么不直接用ANOVA呢因为微生物丰度数据极少满足正态分布和方差齐性大量物种在多数样本里都是0数据分布是极度偏态的。Kruskal-Wallis检验的思路是把所有样本的特征值混合后按大小排序用秩次代替原始数值做比较。它不太受极端值影响适合处理带有大量0和少数高丰度值的数据。LEfSe默认要求这个检验的p值小于0.05才认为这个物种在多个组之间的分布有差异才有资格进入下一轮筛选。实际分析时要注意这一步筛选出来的候选特征可能仍然非常多尤其当样本分组差异很大时可能会有数百个物种通过检验。这很正常后续两步就是进一步压缩这些候选。2.2 第二关Wilcoxon 两两比较解决“谁和谁不同”Kruskal-Wallis检验只回答“几个组之间是否存在差异”但它不会告诉你是A组和B组不同还是A组和C组不同。所以LEfSe第二步会进行Wilcoxon秩和检验做两两组间比较。默认情况下LEfSe会对通过第一步的特征做Wilcoxon检验任意两组之间的差异都需要达到显著水平这个特征才会被保留。对于多组比较还提供了one-against-all或者更严格的策略具体可以在运行参数里控制。这一步有点像是“组间差异一致性检查”避免某些特征只是被单一极端分组拉高了整体p值。我遇到不少同学只记住了LEfSe会跑Wilcoxon却忽略了两两比较的策略差异。如果分组超过三组建议在跑之前想清楚你是想找某一组相对于其他所有组的特异marker还是想找任意两组之间都能分开的marker这个选择会直接影响最终结果。2.3 第三关LDA score 不是 p 值是效应量通过了前两轮检验的特征接下来会进入线性判别分析阶段。线性判别分析的基本思想是寻找一个线性组合让组间差异尽可能大、组内差异尽可能小。在LEfSe里每个特征都会被计算出一个判别系数这个系数经过标准化后就是LDA score通常取log10后的值默认以2作为显著阈值。很多人第一次看LDA条形图会以为横轴大小就是丰度差异大小其实不是。LDA score反映的是这个特征在区分不同组别时的重要性/贡献度含义更接近“效应量”。所以假设某物种在两组中丰度差异明显但组内波动也很大它的LDA score未必很高。反过来一个丰度不是最高、但组间规律性很强的物种LDA score可能会很高。这个特点决定了LEfSe找出来的物种通常具备“组间差异一致性好”的特征这也是为什么审稿人比较认可这套结果。2.4 阈值用多少alpha 和 LDA score 怎么调LEfSe里有两个重要参数统计检验的显著性阈值alpha默认0.05和LDA score阈值默认2.0。前者控制Kruskal-Wallis和Wilcoxon检验的p值标准后者控制最终展示多少个marker。如果样本量不是很大比如每组只有五六例跑出来一个显著marker都没有可以尝试把LDA阈值降到2.0以下比如1.5但不能降得太狠否则结果会变得很“虚”。反过来如果差异物种特别多图纸一片混乱可以把LDA阈值提高到3.0甚至3.5减少图形上的信息量只保留最有区分度的物种。关于多重检验校正LEfSe默认流程其实没有做特别严格的FDR校正所以它的p值只能作为筛选参考不建议当作确定性结论。如果论文需要更严格的统计支撑最好再用ANCOM-BC或ALDEx2做交叉验证。3. 分析前的数据准备3.1 特征表和分组信息长什么样LEfSe分析需要两个核心输入一个是特征丰度表另一个是样本分组信息表。特征丰度表一般长这样第一行是样本名第一列是特征名中间是丰度值Taxon Sample_CTR1 Sample_CTR2 Sample_CTR3 Sample_DZ1 Sample_DZ2 Sample_DZ3 k__Bacteria|p__Firmicutes 120 150 110 80 90 70 k__Bacteria|p__Bacteroidota 200 180 210 320 340 310 k__Bacteria|p__Proteobacteria 50 40 45 160 150 170分组信息表也必须是tab分隔第一列是样本名第二列是分组。如果有二级分组可以在第三列写上subclassSample_CTR1 Control Batch1 Sample_CTR2 Control Batch1 Sample_CTR3 Control Batch2 Sample_DZ1 Disease Batch1 Sample_DZ2 Disease Batch1 Sample_DZ3 Disease Batch2这里必须强调分组信息中的样本名要和特征表的列名完全一致顺序不要求一致但名字不能有多余空格不要加引号。3.2 输入格式常见雷区LEfSe对输入格式比较挑剔最常见的坑有三个。第一缺失值处理不当。特征表里出现空单元格、NA或者NaN后面跑format_input.py时很容易直接报错。处理方式很简单所有缺失值都填成0。第二空格和特殊字符。Excel里做表格后直接复制到txt经常会把tab变成空格或者让样本名多一个看不见的空格。建议在Linux下用head -5和cat -A检查一下分隔符如果列间是^I才是tab。另外特征名里不要带空格、括号、引号全部用下划线替代。第三分类层级前缀不完整。如果后面想画Cladogram图特征名必须按照k__、p__、c__、o__、f__、g__、s__这种层级格式并用|连接。如果只是输入一个OTU编号或ASV编号没有分类注释LEfSe也能跑但Cladogram画不出来或者画出来只有一个孤零零的点没法看。3.3 要不要先做预处理LEfSe本身没有做丰度归一化它直接使用你输入的特征值进行计算。因此如果你的输入是原始测序得到的reads count建议提前转成相对丰度或者用总丰度缩放的方式归一化。很多流程会输出相对丰度表直接拿来跑就可以。另一个重要操作是过滤低丰度特征。微生物组数据可能有上万行特征但绝大多数物种在少数样本里出现一次这种特征进了LEfSe只会拖慢速度对结果贡献也很小。我一般会先过滤一遍在所有样本中相对丰度都低于0.01%的特征直接去掉如果有明显是叶绿体、线粒体、未分类的序列也在这一步剔除。值得注意的是过滤不能太激进。如果样本之间差异比较细微你把低丰度物种全部删掉可能把真正有生物学意义的稀有物种也丢了。建议同时保留原始未过滤表方便后续对比。3.4 在线版还是本地版LEfSe最早提供在线Galaxy版用户上传两个文件、选好参数就能跑适合快速出图。在线版的优点是省去安装环境对于不熟悉命令行的同学非常友好。缺点是大数据量上传比较费时间而且每次都要手动操作重复性不够好。本地版则需要安装LEfSe程序用命令行运行好处是脚本化、可复现批量处理多个数据集的时候效率高。我个人更推荐本地版尤其是当你需要跑几十个样本、要多次调整阈值或者要在服务器上做批量分析的时候。本地版也没有样本量限制只要机器内存扛得住。如果只是两三个小数据集在线版也能满足需求。总之建议两条路都试一下本地跑出来结果不理想时再去在线平台核对一遍这样心里更有底。4. 实操从命令行跑通 LEfSe4.1 安装LEfSe本地安装时建议用conda创建独立环境避免和项目里其他Python环境冲突。操作步骤如下conda create -n lefse python3.9 -y conda activate lefse conda install -c bioconda lefse -y也可以尝试用pip安装pip install lefse安装完成后验证一下脚本是否可用run_lefse.py --help format_input.py --help如果能看到帮助信息说明安装成功。需要提醒的是LEfSe不同小版本之间参数命名略有差异实际使用前先用--help看一眼以你本机版本为准。老教程里的Python 2命令在新环境中不一定能直接跑。4.2 准备两个数据文件我在服务器上习惯建一个工作目录把所有文件放进去比如LEfSe_test/。目录下放两个文件feature_table.txt特征丰度表第一行样本名第一列特征名。group_info.txt分组信息表第一列样本名第二列分组。写完后检查文件格式head -5 feature_table.txt head -5 group_info.txt这一步很关键很多报错都是因为分隔符不是tab。4.3 format_input.py把原始表格转换成 LEfSe 内部格式正式分析前需要先用format_input.py将特征表和分组信息合并成一个LEfSe内部格式文件。基本命令如下python format_input.py feature_table.txt group_info.txt lefse_input.txt -c 1 -u 2其中-c 1表示分组信息文件里class所在的列是第1列-u 2表示可选的subclass在第2列。如果没有subclass就不要加-u。这一步执行完会生成lefse_input.txt。可以用head看一下内容文件格式已经变成LEfSe能识别的结构。如果这里报错说明输入文件格式有问题需要回到上一步检查分隔符和缺失值。4.4 run_lefse.py核心统计检验生成内部格式文件之后就可以跑LEfSe核心分析python run_lefse.py lefse_input.txt lefse_result.txt -l 2 -a 0.05 -w 0.05参数含义-l 2LDA score阈值默认2。-a 0.05Kruskal-Wallis检验的显著性阈值。-w 0.05Wilcoxon检验的显著性阈值。-o 1多组比较时使用one-against-all策略可选。运行结束后会生成lefse_result.txt里面每一行是一个通过筛选的特征包含它对应的分组、LDA score、p值等信息。这个文件是后续画图的输入也是你自己做数据整理时的重要依据。我觉得第一次跑的时候最好把过程输出的日志保留万一结果异常翻日志能看出是在哪一步筛选掉的。4.5 画图LDA条形图和Cladogram图LEfSe自带几个画图脚本最常用的是python plot_res.py lefse_result.txt lda_bar.pdf python plot_cladogram.py lefse_result.txt cladogram.pdfplot_res.py生成的是LDA条形图展示每个显著物种的LDA scoreplot_cladogram.py生成的是Cladogram分支图展示差异物种的分类层级关系。如果需要查看每个显著特征的丰度分布还可以用python plot_features.py lefse_result.txt feature_plots/这会在feature_plots/目录下生成每一个显著特征在不同分组中的丰度分布图适合用来快速核查结果是否合理。画图时也可调整参数比如plot_res.py可以用--format pdf指定输出格式或者通过-f调整图片宽度。实际需求不同可以自己探索。4.6 一个可复现的完整示例这里给一个从零开始的完整命令序列假设已在LEfSe_test/目录下准备好两个输入文件# 第一步合并数据 python format_input.py feature_table.txt group_info.txt lefse_input.txt -c 1 # 第二步跑LEfSeLDA阈值2.0 python run_lefse.py lefse_input.txt lefse_result.txt -l 2 -a 0.05 -w 0.05 # 第三步画图 python plot_res.py lefse_result.txt lda_bar.pdf python plot_cladogram.py lefse_result.txt cladogram.pdf # 第四步查看结果表格 head -20 lefse_result.txt这样一个流程跑下来几十个样本的数据通常几分钟就能完成。如果特征表特别大耗时可能稍长但一般不会超过半个小时。5. 结果怎么看、怎么写进文章5.1 LDA 条形图先看效应量再看颜色LDA条形图是LEfSe最经典的结果。图中每个条形代表一个显著差异物种横轴是LDA score条形越长说明该物种对分组的判别贡献越大。条形颜色对应富集分组红色、绿色、蓝色通常分别代表不同组。看这张图我一般分三步。第一步看整体数量显著物种是不是太多或太少如果太多可以调高LDA阈值再跑一次第二步看颜色分配是否合理各组的marker是否有生物学逻辑第三步挑几个LDA score最高的物种回到相对丰度表里验证一下它们在组间的差异是否真实。颜色相同的条形聚集在一起通常说明这些物种属于同一个高的分类单元协同变化。5.2 Cladogram 分支图从门到属的富集路线Cladogram图看起来像一圈圈的年轮内部是门往外依次是纲、目、科、属、种。图上用不同颜色标记在不同组中显著富集的物种节点大小代表丰度大小圆圈越靠近外圈分类越具体。读懂Cladogram的关键是理解“由内到外”的分类层级。比如某个节点的分类学注释是p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Lactobacillaceae|g__Lactobacillus在图上会显示为一串从内到外的节点如果整条路径都被标记为红色说明从门到属这个完整分支在红色组里都有显著富集。论文里常说的“Firmicutes门及其下属分类在疾病组显著富集”就是根据这张图得出的。需要提醒的是Cladogram图很好看但信息密度高阅读时容易眼花。汇报时最好把关键分类层级单独列出来配合LDA条形图一起说明。5.3 特征丰度图和热图除了两张主图LEfSe还可以输出每个显著特征在分组间的丰度分布图这就是plot_features.py的作用。这类图能直观展示某个物种在不同样本中的高低情况适合放到补充材料。如果想做热图可以把lefse_result.txt里显著的物种提取出来再结合原始丰度表用R的pheatmap画。一般做法是先选出LDA score排名靠前的物种然后按样本做z-score标准化再画热图。这样能把所有显著marker和样本之间的关系一次性展示出来。5.4 论文里怎么描述 LEfSe 结果写论文时LEfSe方法部分可以写成这样“Linear discriminant analysis Effect Size (LEfSe) was applied to identify taxa with significant abundance differences between groups, with the threshold of LDA score set to 2.0 and the significance level set to 0.05.”结果部分则避免只说“某些物种有差异”最好结合具体效应量。比如“LEfSe analysis identified 15 taxa significantly enriched in the disease group, including g_Lactobacillus (LDA score 4.2, P 0.003) and f_Enterobacteriaceae (LDA score 3.8, P 0.01).”我见过不少同学在结果部分放一张LDA图和一张Cladogram图但正文只字不解释这样读者很难抓住重点。图的展示需要文字引导至少要说明图里最大的特征是哪些物种。6. 常见问题与排查记录6.1 报错与解决办法速查表常见报错或异常可能原因处理办法format_input.py 报 “input file” 错误分隔符不是tab或存在空缺值用cat -A检查文件补0或重导出tab分隔文本样本名对不上特征表和分组信息中样本名不一致用 cut -f1 group_info.txtCladogram 画不出来特征名缺少分类层级前缀或 连接结果文件为空Kruskal-Wallis 或 Wilcoxon 无显著特征降低 alpha或检查分组是否合理LDA图条形数量太多差异物种过多提高-l阈值如3.0或3.5运行时间太长特征过多、样本过多先过滤低丰度特征再跑LEfSe结果和普通t检验不一致LEfSe采用多步筛选更保守结合其他差异分析方法交叉验证6.2 分析结果“没有显著 marker”这是最常遇到的问题之一。有时候分组差异确实存在但LEfSe一个显著结果都没跑出来。原因可能有三个第一样本量太少。每组只有两个样本统计检验几乎没有效能很难通过Kruskal-Wallis检验。这种情况可以尝试把LDA阈值降低到1或者接受当前数据不适合LEfSe的事实改用效果更直接的描述性分析。第二组间差异太弱。微生物群落整体差异不大个别物种的丰度变化也不稳定LEfSe的两次秩和检验会把大部分特征过滤掉。可以考虑换用ANCOM-BC等更灵敏的方法或者增加样本量。第三阈值设得太严。默认LDA大于2确实是一个比较高的标准如果数据集噪声大可以先用默认参数跑一遍再降到1.5左右看看趋势。但不能无限降否则就失去LEfSe的意义了。6.3 分组样本量不均或特别小时怎么处理LEfSe对每组最少样本数没有强制规定但实际经验是每组少于3个就会很不可靠少于5个时结果波动很大。如果你手里就是小样本建议同时做以下操作先看各组样本分布是否重叠降低LDA阈值后观察marker是否仍然一致再用随机森林等方法做特征重要性分析看是否和LEfSe结果有交集。如果有的组样本特别多、有的组特别少比如疾病组20个、对照组4个跑出来的差异物种往往会偏向样本量大的那一组。这种情况下可以对大组做随机抽稀比如从20个样本里随机抽4个重复多次看哪些物种在大样本和小样本组合情况下都能稳定出现。这个做法比较耗时间但能明显提高结论可信度。6.4 结果与之前做的差异分析不一致不少人在跑LEfSe之前已经用过DESeq2或者Wilcoxon做了差异物种筛选结果发现两套方法给出的物种列表对不上。这很正常因为LEfSe更关注“效应量”步骤也更严格。DESeq2这类方法针对的是原始计数矩阵基于负二项分布假设LEfSe直接处理相对丰度或标准化数据不依赖这个假设。两种方法得到交集之外的结果不代表谁错了只能说明它们对“差异物种”的定义不同。写论文时可以把两套结果并排展示取交集作为核心结论把只有一种方法检出的物种放到补充材料。这样做比只放一套结果更有说服力。7. 个人心得与避坑经验7.1 LEfSe 不是万能的别把结果当真理我见过一种比较危险的做法把LEfSe跑出来的所有显著物种直接拿去做机制解释甚至当成“因果证据”。实际上LEfSe只能说明在给定分组条件下这些物种的丰度分布存在差异不能说明是因果关系、也不能排除混杂因素。比如疾病组和对照组年龄不匹配跑出来的很多物种可能只是反映年龄差异而不是疾病本身。因此每次跑完LEfSe之后我都会先做一个“反常识检查”这些显著富集的物种是否有文献支持如果一个物种在以往多个研究中都被认为和疾病相关那这次结果就比较可信如果完全不符合已有认知先别急着高兴很可能数据里存在某个未控制的混杂因素。7.2 实际运行中的几个小习惯我自己的操作习惯是所有原始数据先做标准化和过滤得到相对丰度表后再做一次分组检查确认每组样本数量、样本名是否匹配。然后先用默认参数跑一遍记录结果再微调LDA阈值跑第二遍比较两次结果的重合情况。如果两次差异很大说明数据本身不太稳定需要小心。此外lefse_result.txt这个结果文件一定要保留。很多后续操作都要用到它比如提取显著物种做热图、做韦恩图、或者给其他工具做输入。不要只留PDF图丢了中间结果后面会很被动。7.3 最后一个实用小技巧如果你想在论文里把LDA条形图画得更精致可以直接读取lefse_result.txt用R的ggplot2重新绘制。LEfSe自带的图虽然方便但配色和排版不一定符合期刊要求。自己画图时把物种名按照LDA score降序排列用不同颜色代表分组再调整字体和坐标轴效果会好很多。另外多组比较时LEfSe的Cladogram图可能非常拥挤节点标签重叠严重。这时候与其硬调图片尺寸不如拆开画按每个显著的门或者纲单独画一个子图读者看起来更清楚审稿人也更舒服。希望这些实测经验能帮你少走一点弯路。
