microeco：一个用于微生物群落生态学数据挖掘的R软件包

microeco: an R package for data mining in microbial community ecology

发表时间：2020年12月17

摘要

在使用高通量测序技术的微生物群落生态学研究中产生了大量的测序数据，尤其是基于扩增子测序的群落数据。在对扩增子测序数据进行初步的生物信息学分析后，根据操作分类单元和分类分配表进行后续的统计和数据挖掘仍然是复杂和耗时的。为了解决这个问题，我们提出了一个集成的R包–‘microeco’作为处理微生物群落和环境数据的分析管道。这个包是基于R6类系统开发的，结合了微生物群落生态学研究中的一系列常用和先进的方法。该软件包包括数据预处理、分类群丰度绘制、维恩图、α多样性分析、β多样性分析、差异丰度测试和指标分类群分析、环境数据分析、null模型分析、 络分析和功能分析等类。每个类的设计都是为了提供一套可以让用户方便使用的方法。与微生物生态学领域的其他R包相比，microeco包使用起来快速、灵活、模块化，为研究人员提供了强大而便捷的工具。microeco包可以从CRAN（The Comprehensive R Archive Network）或github（
https://github.com/ChiLiubio/microeco）安装。

引言

高通量测序技术在微生物生态学中的应用产生了大量的测序数据。现在，微生物群落生态学中的测序数据分析可以任意分为生物信息学分析和随后的统计分析。生物信息学分析是一个典型的计算密集型工作，unix/linux服务器适合这一操作。在根据扩增子测序数据获得必要的操作分类单元（OTU）表、分类分配和系统发育树后，下游的数据分析往往是多样化、复杂化和耗时的。虽然在生物信息学平台上有一些统计和绘图方法，如QIIME（Caporaso等人，2010）、mothur（Schloss等人，2009）、RDP（http://rdp.cme.msu.edu/）和SiLVAngs（
https://ngs.arb-silva.de/silvangs/），但要有效地进行数据挖掘，需要更强大、灵活和全面的工具。

R编程语言及其用于统计计算的软件包系统在科学应用中因其强大和灵活而脱颖而出。除了经典的生态学软件包，如vegan（Oksanen等人，2019），还设计了几个R软件包来进行基于测序的微生物群落数据的复杂分析，如phyloseq（Mcmurdie和Holmes 2013）、microbiome（
https://github.com/microbiome/microbiome）、microbiomeSeq（
http://www.github.com/umerijaz/microbiomeSeq）、ampvis2（
https://madsalbertsen.github.io/ampvis2/reference/index.html）、MicrobiomeR（
https://github.com/vallenderlab/MicrobiomeR）和Rhea（Lagkouvardos等人，2017）。在软件设计方面，phyloseq软件包在S4类的基础上发展良好。microbiome、microbiomeSeq、MicrobiomeR和ampvis2包都依赖于phyloseq包的phyloseq类，涵盖了微生物群落生态学研究中的一系列功能（表S1，支持信息）。然而，它们都缺乏与分析方法相关的封装，缺乏一些重要的、前沿的方法（表S1，支持信息）。因此，这些封装框架对于目前微生物群落研究的数据挖掘管道来说是不够的。在采用多种统计和绘图方法进行数据分析的情况下，迫切需要一个模块化的软件来减少用户学习和使用的时间成本。Rhea软件包有一个简短而清晰的工作流程框架，但它缺乏许多重要的方法，也不是一个传统的基于函数的软件包。

在此，我们基于R6类系统（Chang 2019）创建了一个名为microeco的包，它提供了封装的面向对象编程范式。为分析方法的每一部分创建了类，以使包的框架模块化、清晰和简短（图1）。每个类都封装了一系列的功能和不同的算法。我们整合了一些常用的方法，使软件包支持广泛的微生物群落分析，如LEfSe（Segata等人，2011）、冗余分析（RDA）、共现 络分析、功能预测和空模型分析（Stegen等人，2013）。我们还开发了一些独特的高级功能，如维恩图中的分类群组成分析以及原核生物和真菌的功能冗余计算。

microtable class

trans_abund和trans_venn class（小编注释：该函数用于Overlap 类群分析）

为了可视化分类群的丰度变化，我们创建了trans_abund类来实现柱状图、boxplot、热图和饼图的绘制方法。trans_venn类是为文氏图设计的，用于解读组间共享的或特殊的类群，包括对5个以上组的分析。与其他维恩分析包，如VennDiagram（Hanbo Chen 2018）相比，我们还增加了每个特定或共享部分的OTU数量或序列数量的百分比。一般来说，研究人员对文氏分析结果中的特定和共同部分的分类群组成感兴趣。因此，我们开发了独特的函数trans_venn_com()来分析每个部分的分类群组成（Mendes等人，2011；Andrew等人，2012）。该函数可以将每个馏分中的OTU组成转化为群落格式表，并返回一个新的微表对象，用于快速绘制每个馏分的分类群组成（见补充材料例2）。

trans_alpha class（小编注释：该函数用于alpha多样性数据格式转换、绘图和差异性统计）

trans_alpha类被设计用来将microtable对象的alpha多样性表转化为其他格式的统计和绘图。函数cal_diff()可以通过Kruskal-Wallis秩和检验或ANOVA（方差分析）对所有α多样性指数进行组间差异检验，并进行多重比较。下面的函数plot_alpha()用于展示α多样性数据，并通过多重比较或成对比较直接添加显著性。

trans_beta class

微生物群落数据的生态学分析中的一个关键部分是β多样性，这可以用trans_beta类来进行。trans_beta类目前实现了几种常用的无约束排序方法，如主成分分析、主坐标分析（PCoA）和非计量多维缩放（NMDS）等。此外，同一组或不同组的样本之间的距离可以表明组内或组间的群落异同（Lim等人，2015）。函数cal_group_distance()和plot_group_distance()可以分别进行群落距离的变换和绘制。聚类分析和绘图也在函数 plot_clustering()中实现。此外，为了简化包罗万象的多元方差分析（perMANOVA）的使用，我们开发了函数cal_manova()，用于perMANOVA的总体比较、配对组比较和多因子比较。

trans_diff class（小编注释：该函数用于筛选不同分组间差异的物种或类群）

指标类群的鉴定对于解释不同群体间群落结构差异的生物学机制具有重要意义。随着微生物群落生态学中测序数据复杂性的增加，评估各组间有显著差异的分类群是一个挑战。因此，一些工具，如LEfSe（Segata等人，2011），结合了监督机器学习和差异丰度测试的优势，以确定区分一个群体与其他群体的重要指标类群。trans_diff类目前实现了三种著名的方法。LEfSe（Segata等人，2011），随机森林（An等人，2019）和metastat（White, Nagarajan and Pop，2009）。LEfSe结合了线性判别分析（LDA）和差异丰度测试来寻找指标类群。与Python版本相比，trans_diff类中的LEfSe方法是用R代码重新实现的，以减少对软件的依赖。同时，结果的绘制（LDA得分的柱状图和支系图）也易于使用该软件包进行调整。我们整合了随机森林分析和差异丰度测试，与randomForest软件包相比，增加了随机森林方法的便利性和力量。metastat方法在两组之间的差异丰度测试中特别有用（White, Nagarajan and Pop 2009）。我们实现了metastat方法以减少文件格式转换的工作量，并增加了相应的绘图功能。

trans_env class

评估环境因素对微生物群落结构的影响，对于推断支配群落组装的基本机制至关重要。在trans_env类中，冗余分析（RDA）和基于距离的RDA（db-RDA）在cal_rda()中基于素描函数rda()实现。为了方便RDA结果的可视化，我们在trans_rda()和plot_rda()中分别加强了转换和绘图的方法。环境数据和β-多样性距离矩阵之间的mantel检验可以通过cal_mantel()对所有环境变量和所有距离矩阵进行方便的计算。相关热图是展示分类群丰度和环境因素之间相关性的重要方法。这可以通过使用trans_env类的两个步骤完成。首先，分类群和环境因素之间的相关性可以用cal_cor()计算所有样品或在不同组内进行。然后，可以使用ggplot2包风格的相关热图或带有聚类图的pheatmap包风格的相关热图进行plot_corr()，即使相关分析是针对不同组的。

trans_nullmodel class

近几十年来，系统发育分析和null模型的整合，通过增加系统发育的维度，更有力地促进了生态位和中性（niche and neutral）对群落组合影响的推断（Stegen等人，2013）。trans_nullmodel类提供了一个封装，包括系统发育信 、β平均成对系统发育距离（βMPD）、β平均最近分类群距离（βMNTD）、β最近分类群指数（βNTI）、β净相关度指数（βNRI）和基于Bray-Curtis的Raup-Crick（RCbray）的计算。系统发育信 分析的方法是mantel correlogram（Liu等人，2017）。betaMNTD和betaMPD的算法经过优化，比picante包（Kembel等人，2010）的算法更快。RCbray和betaNTI（或betaNRI）之间的组合可用于推断特定假设下主导群落集合的每个生态过程的强度（Stegen等人，2013；Liu等人，2017）。这可以通过函数cal_process()解析每个推断过程的百分比来实现。

trans_network class

微生物生态学中的共现模式分析是一个热门话题，通常是应用 络分析的方式进行分析。在trans_network类中，提供了三种 络构建方法：相关 络（Deng等人，2012）、SPIEC-EASI（SParse InversE Covariance Estimation for Ecological Association Inference） 络（Kurtz等人，2015）和基于概率图形模型（PGM）的 络（Tackmann, Matias Rodrigues and Mering 2019）。对于基于相关的 络（图2中的步骤1），三种相关计算方法是可选的，包括cor.test函数、WGCNA包（Langfelder和Horvath 2008）和SparCC（Friedman和Alm 2012）。相关系数的优化可以选择使用基于随机矩阵理论的方法（Deng等人，2012）。另一种 络构建方法是SpiecEasi R包中的SPIEC-EASI（Kurtz等人，2015）。我们实施的第三种 络方法是PGM 络，它可以自动调用系统中的Julia语言和FlashWeave包（Tackmann, Matias Rodrigues and Mering 2019）。用函数save_network()和默认设置保存 络可以生成network.gexf文件，该文件可以用Gephi软件（https://gephi.org）打开，用于 络可视化。模块分类、 络拓扑属性和节点拓扑属性的计算都是基于igraph软件包（Csardi and Nepusz 2006）开发的。快速贪婪的模块优化算法通过igraph中的cluster_fast_greedy()函数实现了模块的分类。在PCA分析的基础上，提取每个模块的OTU丰度矩阵的第一个主成分作为模块的eigengenes，以揭示高阶组织信息（Deng等人，2012）。模块的eigengene分析由函数cal_eigen()提供。此外，我们还开发了函数subset_network()用于提取 络的一部分，函数cal_sum_links()用于总结任何分类等级的分类群之间的联系。

trans_func class

功能分析是微生物群落数据分析中一个有吸引力的部分，主要来自其对生物学问题的可解释性。在微生物群落生态学中，探索功能冗余是一个挑战，因为很难识别OTU或物种的功能特征。为了使这种分析在trans_func类中可行，我们开发了cal_spe_func()函数，将OTU的分类分配与原核生物的FAPROTAX数据库（Louca, Parfrey and Doebeli 2016）的分类信息或与真菌的Funguild数据库（Nguyen et al. 2016）的分类信息自动匹配（图3的步骤1）。FAPROTAX数据库是根据Bergey’s Manual of Systematic Bacteriology和相关文献的信息建立的，所以原核生物性状鉴定的结果是可靠和保守的。随后，可以通过函数cal_spe_func_perc()来计算群落或 络模块中与每个性状相关的类群（丰度未加权）或个体（丰度加权）的百分比（图3中第2或3步）。我们认为一个功能性状的百分比越高，说明该性状在群落或 络模块中的功能冗余度越高。

为了方便预测和分析群落的功能潜力，我们整合了Tax4Fun R软件包（A?hauer等人，2015）用于元基因组代谢途径预测（图3中的第6步）和FAPROTAX软件（
http://www.loucalab.com/archive/FAPROTAX）来预测群落对生物地球化学循环的功能潜力。在获得群落的功能潜力后，通过使用microeco软件包很容易进行以下统计分析（图3中的第7-9步）。

结果和讨论

首先，在图4和补充材料例2中显示了几个操作。我们进行了韦恩图分析，并分析了每个共享或特定部分的OTU组成（图4A）。在维恩图中加入OTU的序列 百分比来说明每个馏分中OTU的丰度信息（图4A）。结果表明，CW、IW和TW之间共享的OTU以及这些共享的OTU的丰度占总OTU数量和丰度的比例最大。为了探索每个馏分中的分类群组成，我们对结果进行了转换，得到了每个馏分的属级的OTU组成（图4B）。这种方法可以帮助解释一些分类群（这里是指属）是由特定的还是共享的OTU组成的。例如，在共享部分CW-IW-TW中，Planctomyces属的比例比其他组相对较高。接下来的例子是组间α-多样性的配对比较，这表明IW中的Chao1明显高于CW（图4C）。为了进一步显示各组间的β多样性是否存在较大差异，如补充材料例2所示，进行了基于Bray-Curtis距离的PCoA（图4D）。我们还用trans_nullmodel来计算βNRI，以评估不同组内是否存在不同的系统发育周转模式。结果表明，TW的基本系统发育周转的强度明显高于CW和IW（图4E）。用LEfSe（图4F和G）来识别指示性类群，表明指示性类群在不同组中是不同的。例如，Proteobacteria and Firmicutes 在CW中明显富集（图4F和G）。为了进一步探讨环境因素如何影响这些指示性类群，我们采用trans_env类将非生物因素与这些类群的丰度相关联，并进行相关热图（图4H）。我们还评估了分类群和环境因素的聚类情况，以揭示一些分类群比其他分类群具有更相似的生态位偏好的程度。结果清楚地表明，相关热图可以提供有用的帮助，使不同的指标类群可能受到环境因素的影响（图4H）。