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等人，2022），还设计了几个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/vallerlab/MicrobiomeR）和Rhea（Lagkouvardos等人，2022）。在软件设计方面，phyloseq软件包在S4类的基础上发展良好。microbiome、microbiomeSeq、MicrobiomeR和ampvis2包都依赖于phyloseq包的phyloseq类，涵盖了微生物群落生态学研究中的一系列功能（表S1，信息）。然而，它们都缺乏与分析方法相关的封装，缺乏一些重要的、前沿的方法（表S1，信息）。因此，这些封装框架对于目前微生物群落研究的数据挖掘管道来说是不够的。在采用多种统计和绘图方法进行数据分析的情况下，迫切需要一个模块化的软件来减少用户学习和使用的时间成本。Rhea软件包有一个简短而清晰的工作流程框架，但它缺乏许多重要的方法，也不是一个传统的基于函数的软件包。 在此，我们基于R6类系统（Chang 2022）创建了一个名为microeco的包，它提供了封装的面向对象编程范式。为分析方法的每一部分创建了类，以使包的框架模块化、清晰和简短（图1）。每个类都封装了一系列的功能和不同的算法。我们整合了一些常用的方法，使软件包广泛的微生物群落分析，如LEfSe（Segata等人，2011）、冗余分析（RDA）、共现网络分析、功能预测和空模型分析（Stegen等人，2013）。我们还开发了一些独特的高级功能，如维恩图中的分类群组成分析以及原核生物和真菌的功能冗余计算。 microtable class microtable class是一个基本类，用于创建存储输入文件的对象和后续的数据预处理。与phyloseq包中的phyloseq类相比，microtable类很简短。例如，创建对象时需要的文件格式是用户友好的data.frame 类，用于所有OTU表、分类表和样品信息表（见补充材料例1）。用于数据预处理的内置函数包括tidy_dataset(), filter_pollution(), rarefy_samples(), merge_samples(), merge_taxa()。通过应用函数tidy_dataset()，对象中所有的基本文件都可以被修剪，并具有一致的信息。如果各样本的序列数差别很大，建议使用rarefy_samples()使每个样本的序列数相等（小编注释：rarefy_sample 函数功能是对测序数据OTU进行抽平处理），以减少测序深度对α-和β-多样性计算的影响。运行cal_abund()可以自动计算每个分类等级的分类群丰度，并返回一个包含所有表格的列表（taxa_abund）。函数cal_alphap()用于计算α多样性，包括Chao1、Shannon-Wiener指数、Simpson指数和系统发育的α多样性（Faith 1992）。此外，函数cal_betap()可用于获得β多样性的距离矩阵，包括Bray-Curtis、Jaccard和UniFrac（Lozupone和Knight 2005；陈俊2022）。 trans_abund和trans_venn class（小编注释：该函数用于Overlap 类群分析） 为了可视化分类群的丰度变化，我们创建了trans_abund类来实现柱状图、boxplot、热图和饼图的绘制方法。trans_venn类是为文氏图设计的，用于解读组间共享的或特殊的类群，包括对5个以上组的分析。与其他维恩分析包，如VennDiagram（Hanbo Chen 2022）相比，我们还增加了每个特定或共享部分的OTU数量或序列数量的百分比。一般来说，研究人员对文氏分析结果中的特定和共同部分的分类群组成感兴趣。因此，我们开发了独特的函数trans_venn_com()来分析每个部分的分类群组成（Mes等人，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等人，2022）。函数cal_group_distance()和plot_group_distance()可以分别进行群落距离的变换和绘制。聚类分析和绘图也在函数 plot_clustering()中实现。此外，为了简化包罗万象的多元方差分析（perMANOVA）的使用，我们开发了函数cal_manova()，用于perMANOVA的总体比较、配对组比较和多因子比较。 trans_diff class（小编注释：该函数用于筛选不同分组间差异的物种或类群） 指标类群的鉴定对于解释不同群体间群落结构差异的生物学机制具有重要意义。随着微生物群落生态学中测序数据复杂性的增加，评估各组间有显著差异的分类群是一个挑战。因此，一些工具，如LEfSe（Segata等人，2011），结合了监督机器学习和差异丰度测试的优势，以确定区分一个群体与其他群体的重要指标类群。trans_diff类目前实现了三种著名的方法。LEfSe（Segata等人，2011），随机森林（An等人，2022）和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等人，2022）。betaMNTD和betaMPD的算法经过优化，比picante包（Kembel等人，2010）的算法更快。RCbray和betaNTI（或betaNRI）之间的组合可用于推断特定假设下主导群落集合的每个生态过程的强度（Stegen等人，2013；Liu等人，2022）。这可以通过函数cal_process()解析每个推断过程的百分比来实现。 trans_network class 微生物生态学中的共现模式分析是一个热门话题，通常是应用网络分析的方式进行分析。在trans_network类中，提供了三种网络构建方法：相关网络（Deng等人，2012）、SPIEC-EASI（SParse InversE Covariance Estimation for Ecological Association Inference）网络（Kurtz等人，2022）和基于概率图形模型（PGM）的网络（Tackmann, Matias Rodrigues and Mering 2022）。对于基于相关的网络（图2中的步骤1），三种相关计算方法是可选的，包括cor.test函数、WGCNA包（Langfelder和Horvath 2008）和SparCC（Friedman和Alm 2012）。相关系数的优化可以选择使用基于随机矩阵理论的方法（Deng等人，2012）。另一种网络构建方法是SpiecEasi R包中的SPIEC-EASI（Kurtz等人，2022）。我们实施的第三种网络方法是PGM网络，它可以自动调用系统中的Julia语言和FlashWeave包（Tackmann, Matias Rodrigues and Mering 2022）。用函数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 2022）的分类信息或与真菌的Funguild数据库（Nguyen et al. 2022）的分类信息自动匹配（图3的步骤1）。FAPROTAX数据库是根据Bergey's Manual of Systematic Bacteriology和相关文献的信息建立的，所以原核生物性状鉴定的结果是可靠和保守的。随后，可以通过函数cal_spe_func_perc()来计算群落或网络模块中与每个性状相关的类群（丰度未加权）或个体（丰度加权）的百分比（图3中第2或3步）。我们认为一个功能性状的百分比越高，说明该性状在群落或网络模块中的功能冗余度越高。 为了方便预测和分析群落的功能潜力，我们整合了Tax4Fun R软件包（Aßhauer等人，2022）用于元基因组代谢途径预测（图3中的第6步）和FAPROTAX软件（http://www.loucalab.com/archive/FAPROTAX）来预测群落对生物地球化学循环的功能潜力。在获得群落的功能潜力后，通过使用microeco软件包很容易进行以下统计分析（图3中的第7-9步）。 结果和讨论 microeco软件包中的案例研究涉及中国湿地土壤的16S rRNA基因扩增子测序数据集（An等人，2022）和根瘤微生物组研究的ITS扩增子测序数据集（Gao等人，2022）。详细的样品信息在样品信息表（sample_info_16S和sample_info_ITS）中进行了描述。16S rRNA基因扩增子测序数据集被用来演示本文中的例子。数据集中的样本可分为中国内陆湿地（IW）、沿海湿地（CW）和西藏高原湿地（TW）。我们提供了一个详细的在线教程来展示microeco软件包的使用（https://chiliubio.github.io/microeco/）。在此，我们重点介绍了几个重要的例子和相应的代码（见补充材料），可以应用这些代码来实现图4-6的结果。这里分析的这个例子数据的子集涉及90个样本，每组30个（CW、IW和TW）。我们首先创建了microeco表格对象，并按照补充材料例1中的代码进行了基本的预处理操作。在去除未分配到 "k__Archaea "或 "k__Bacteria "界中的分类群或具有 "线粒体 "或 "叶绿体 "分配信息的分类群后，总共留下了13 296个OTU。对OTU表进行了抽平化处理，使各样品的序列号相同（每个样品有10 000个序列）。分类群丰度、α多样性和β多样性的距离矩阵都是为下游分析计算的（补充材料例1）。 首先，在图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）。 为了进一步研究物种共现模式，我们将trans_network类应用于相关网络分析（补充材料例3）。用特定的阈值（spearman相关：R=0.7，P=0.01）构建相关网络，并用gephi可视化（图5A）。我们用和弦图来显示门间内部和之间的联系数（图5B）。为了探索网络中节点的重要性，我们根据模块内和模块间的连接性来分配节点的角色（图5C）。这种方法已被广泛应用于网络分析中，并被验证为对网络中节点作用的综合评价（Deng等人，2012；Shi等人，2022）。然后计算模块的eigengenes，并与非生物因素相关（图5D），揭示了模块的OTU组成可能由不同的环境因素主导。需要指出的是，在非生物因素过滤所带来的强系统发育聚类的条件下，相关网络可能由大量与生物相互作用无关的边组成。因此，最好使用相关网络来解释物种生态位的相似性和其背后的机制，而不是生物相互作用（Röttjers和Faust 2022）（小编注释：共现网络中的边数只能表征微生物之间的潜在相互作用的强和弱，不能直接视为为微生物之间的真实相互作用）。 功能分析对于揭示物种共现模式和群落组合的结构机制具有重要意义。以前的研究几乎没有破译过微生物共现网络模块的功能概况，主要是由于方法上的困难。我们使用FAPROTAX数据库来识别物种特征，并分析了网络模块的功能冗余（图6A）。结果显示，M2（第二大模块）有很高比例的种群（OTU）具有好氧化合功能。在M4中没有发现N-循环和S-循环相关的种群。然后，我们计算了群落的功能冗余度，并将功能冗余度与非生物因素相关联，以探索非生物因素如何影响功能冗余度（图6B）。结果表明，电导率和MAT（年平均温度）都与亚硫酸盐呼吸和光合作用相关的功能冗余有极其显著的正相关关系。最后，我们用综合函数预测代谢途径的丰度，并在CW、IW和TW之间进行LEfSe（图6C）。CW、IW和TW中LDA得分最高的分别是遗传信息处理、新陈代谢和环境信息处理，这意味着每个区域的富集和指示性途径反映了与群落结构相关的独特的生物和基因模式。这一预测性的、有意义的结果可能为元基因组测序研究的设计和数据挖掘带来更多的希望和更深的思考。 有必要指出的是，随着研究的复杂性和深度的增加，尤其是在最近十年，统计和绘图方法在微生物生态学中极为多样。尽管我们在microeco包中打包了许多常用的和独特的方法，但不可能考虑所有的方法，也不可能满足这个领域的所有要求。我们设计microeco包的主要理念是简单而不失力量。因此，关键是如何让研究人员快速、方便地在该领域进行数据挖掘。我们选择了先进的和有代表性的方法来提高功率，并将这些方法分为几类，以使框架清晰。 创建microeco项目的目标之一是帮助研究人员缩小在学习分析方法的难度和对快速数据挖掘的需求之间一直存在的差距，特别是对新手程序员。目前，使用microeco软件包可以实现对许多热点和重大科学问题的数据探索，如物种共现模式（图2）和功能概况（图3）。使用microeco软件包的另一个好处是，它对研究者来说具有高度的可重复性，便于在出版物中描述复杂的方法。 结论 microeco软件包的特点是：(i)高度模块化，每个类涵盖了分析的各个部分，软件包的结构和功能易于理解、记忆、搜索和使用；(ii)高度灵活性，每个部分都实现了不同的算法或方法，使分析具有灵活性；(iii)快速行动，一些算法已被优化；(iv)功能强大，实现了一些前沿的方法，如LEfSe、网络分析、空模型和功能冗余计算。