转自玉鑫之家，编辑: 刘驰 责编: 孟凡正 审核: 邵瑞鑫 统筹：杨青华
【文献速递】随着微生物组领域实验设计的日益复杂、数据量的激增，以及组学数据类型的多样化，高效的统计分析与可视化面临巨大挑战。福建农林大学环境基因组团队在权威期刊《Nature Protocols》发表了题为“A workflow for statistical analysis and visualization of microbiome omics data using the R microeco package”的方法类文章。福建农林大学资源与环境学院博士毕业生刘驰为第一作者（现为河南农业大学农学院玉鑫实践团成员），李香真教授和姚敏杰教授为通讯作者。该方法文章提供了一个基于R语言microeco包，用于微生物多组学数据统计分析和可视化方法全面的、可重复的参考流程。详细的中文翻译文稿见：https://zhuanlan.zhihu.com/p/1950934368441831618
【处理与方法】本方案要求最低R语言版本为4.4.0。各个所需的R包的详细安装步骤提供在了GitHub仓库的'Step1_Install_packages.R'文件中（见"Data and code availability"章节）。采用一项涉及种植与施肥双因素处理的田间作物试验数据来演示本流程。种植处理包含大豆-玉米轮作（RC）和玉米连作（CC）两种模式；施肥处理设CK（不施肥）、NPK（施用无机肥）、NPKS（NPK+秸秆还田）三种方式，共形成六个处理组：RC-CK、RC-NPK、RC-NPKS、CC-CK、CC-NPK、CC-NPKS，每组设5个生物学重复。扩增子测序数据涵盖大田土、根际土和根系（用于分析内生菌）3个部分；宏基因组测序与非靶向代谢组学数据则包含大田土和根际土2个类型。
图1 | 本方案的整体框架。在完成R语言环境配置（步骤1）后，首先进行扩增子测序数据生物信息学分析结果的预处理与标准化（步骤2-4）；随后依次开展核心微生物组分析（步骤5-6）、α多样性分析（步骤7-9）、β多样性分析（步骤10-11）、差异丰度检验分析（步骤12-18）以及机器学习分析（步骤19-22）。宏基因组测序数据分析分别针对HUMAnN3（步骤23-24）和Kraken2/Bracken（步骤25-26）流程展开。非靶向代谢组学数据分析包含数据预处理（步骤27）、机器学习分析（步骤28）、差异丰度检验（步骤29）以及与微生物类群的关联分析（步骤30）。
【流程整体概况】本方案描述了一套针对扩增子测序数据、宏基因组测序数据和代谢组学数据进行统计分析及可视化的模块化方法（图1）。工作流程的各章节根据分析内容类型进行组织，特别将扩增子测序数据分析划分为多个阶段以便于理解与操作，这些阶段及其对应步骤将在后续章节详细阐述。我们选取了一个具有多处理因子和环境元数据的代表性案例研究，该案例足以展示方法的各种应用场景。在设计整个流程时，我们特别关注了方法的组合使用——这正是微生物群落数据分析中的常见场景。方法组合凸显了软件包的灵活性与高效性，从而为用户的分析过程提供更多选择。对于扩增子测序数据，由于样本间测序深度存在差异，需通过数据标准化避免测序深度对样本比较的影响，因此每个分析步骤都根据经验知识应用了相应的数据标准化方法。例如核心ASV分析、α和β多样性计算使用重抽样方法，而差异丰度检验则采用尺度因子标准化方法。
【主要内容】
第一阶段:安装R包●耗时约30分钟
第二阶段:产生用于扩增子数据后续分析使用的microtable对象●耗时约10分钟
2.导入生物信息分析数据并生成microtable对象。
3.microtable对象预处理。
4.标准化处理。
第三阶段:核心微生物分析●耗时约10分钟
5.识别各生态位下的核心ASV。
6. 各生态位下核心ASV的共享与独有分析。
第四阶段:α多样性分析●耗时约5分钟
7. Alpha多样性的计算。
8. Alpha多样性的差异检验和可视化。
9. Alpha多样性与环境因素的相关性分析。
图2 | 核心物种与α多样性分析。
第五阶段：Beta多样性分析●耗时约10分钟
10. 不涉及环境因子的β多样性分析。
11. 涉及环境因子的β多样性分析。
图3 | Beta多样性分析。
第六阶段:差异丰度分析●耗时约3.5小时
12. ASV水平差异丰度检验。
13. ASV水平不同方法结果比较。
14. 高分类层级差异丰度检验。
15. 差异类群与环境因子相关性分析。
16. 生成模拟ASV丰度用于定制化基准测试。
17. 模拟数据差异检验。
18. 模拟数据差异检验结果可视化。
图4 | 差异丰度检验。
第七阶段:机器学习模型●耗时约15分钟
19. ASV水平下的分析。
20. 属水平下的分析。
21. 多模型分类性能比较。
22. 属水平下的回归分析。
图5 | 机器学习分析结果。
第八阶段：宏基因组数据分析（HUMAnN3软件）● 耗时约10分钟
23. 导入HUMAnN分析结果。
24. 代谢通路的统计分析与可视化。
第九阶段：宏基因组数据分析（Kraken2/Bracken软件）● 耗时约10分钟
25. 导入Bracken软件结果。
26. 对导入丰度数据的统计分析与可视化。
图6 | 宏基因组数据分析。
步骤10：代谢组学数据分析 ● 耗时约20分钟
27.导入代谢组数据。
28. 代谢组数据的机器学习分析。
29. 代谢组数据差异检验。
30. 代谢物与菌属相关性分析。
图7 | 代谢组学数据分析。
本研究得到国家自然科学基金项目（42077206、U22A20608、42477117）、国家重点研发计划（2018YFE0107000）和福建农林大学科技创新专项基金（KFb22072XA）的资助。
