人体肠道宏基因组生物信息分析方法

肠道是含有微生物物种数量最多、功能最全面的微生态群落器官，据报道，人类肠道中有多达39万亿个细菌细胞，是人类细胞数量的1.3倍，而且这些细菌编码的基因数量是人类基因的100多倍[1]。人体肠道微生物的定殖显著影响人类的生理和营养，对人们生活至关重要。先前研究已经表明肠道微生物有助于从食物中获取能量，调节免疫力、体内激素平衡、中枢神经系统成熟以及行为特征。此外，肠道微生物组的变化也与人类多种疾病相关，如肥胖[2-5]、代谢综合征[6]、类风湿性关节炎[7]、炎症性肠病[8-9]、阿尔茨海默病[10]、糖尿病[11]、流感病毒感染[12]、症状性动脉粥样硬化[13]和肝脏疾病[14-17]。

要了解肠道微生物对人类健康的影响，需阐明肠道微环境的组成、多样性和功能[1]。基于DNA的测序分析克服了传统微生物培养技术的局限性，并取得了新的进展，极大地扩大了目前对肠道微生物多样性的认识[18]。由美国国立卫生研究院启动的人类微生物组计划从300个美国人的不同人体部位产生了大量16S rRNA基因数据集，其中显示人类肠道微生物群主要由2个门主导，拟杆菌门和厚壁菌门。然而，16S rRNA基因基于扩增片段的分析只能检测到属和以上级别的微生物，这意味着分类学多样性与功能将在种及其以下水平(即亚种和株)被忽略。宏基因组测序的目标是整个基因组DNA序列，体现了替代16S rRNA基因测序的强大策略，用于分析复杂的微生物群落。宏基因组学研究已应用于人类肠道微生物群，以产生大量可用于表征微生物群落的组成和功能的数据[19]。更重要的是，基于全基因组DNA的测序可以检测物种水平的微生物，并且比基于16S rRNA基因的方法更适合于肠道微生物组和人类疾病之间的关联研究。

宏基因组物种注释也叫宏基因组物种比对，是将经过拆分、质控、去宿主后的高质量clean reads，与参考基因组或marker基因做比对。宏基因组物种注释的参考微生物包括细菌、古细菌、病毒和真菌四大类，它们都来源于当前流行的各大数据库，包括NCBI数据库(http://www.ncbi.nlm.nih.gov/genome)、人体微生物计划(http://www.hmpdacc.org)、GOLD数据库(https://gold.jgi.doe.gov)、酵母基因组数据库(http://www.yeastgenome.org)、真菌基因组行动(http://www.broadinstitute.org)、人类口腔微生物数据库(http://www.homd.org/)、真菌数据库(http://fungidb.org)。随着生物信息分析行业的发展，越来越多好的软件被开发出来，比对软件只是其中一小类，常用短序列比对软件包括：SOAP[20]、BWA[21]、Bowtie[22]、Martin[23]等。有人使用了RefCov (http://gmt.genome.wustl.edu/gmt-refcov)软件，对6种常见短序列比对软件做了评估，主要通过统计它们的覆盖广度和覆盖深度。这项研究表明，这些比对软件都可以在肠道宏基因组数据中发挥各自的价值，其中CLC aligner软件的计算时间最短，使用的内存最少。宏基因组物种注释还需要基于比对结果，获得相应的物种分类信息。虽然现在可用的数据库很多，但是环境中还有许多微生物分类信息没有收录在参考数据库中。目前有少量软件既可以比对微生物基因组，又可以对微生物分类，还可以统计微生物丰度，例如MetaPhlAn2[24]软件、Kraken[10]软件、Genometa[25]软件。其中Kraken软件是基于Kmer的策略，能非常快速、准确找出复杂环境中的微生物基因组，并且物种丰度计算的速度比Metaphlan2软件快11倍多，其核心是Kraken有一种特殊数据库，用以预先计算基因组中包含的特殊Kmer序列，而Genometa有很好的交互性。

2011年，肠型(Enterotypes)的概念首次在《自然》杂志上由Arumugam等[26]提出，该研究发现可以将人类肠道微生物组分成稳定的3种类型，因为这3种类型不受年龄、性别、体重以及地域限制，表现出较高的稳定性，与血型具有很高的相似性，所以将其定义为肠型。2012年，在其他灵长类动物中也发现了肠型的存在[27]。但是由于肠型并不是完全离散的分布，有研究发现在人肠道微生物组中存在4种可互相转换的肠型[28]，另一项针对美国人群的研究则报道称发现了2种肠型，这些结果也引起了国际上对于肠型概念的广泛讨论与争议。

肠型是肠道微生物组在高维数据空间中呈现出的客观存在的聚集效应。目前流行的肠型计算方法有2种，一种是基于样品间的Jensen-Shannon距离，利用围绕中心点划分算法(PAM)进行聚类，最佳分类数目通过Calinski-Harabasz (CH)指数确定[29]；另一种则是直接基于物种丰度，利用狄利克雷多项混合模型(DMM)进行肠型分类[30]。之后，肠型的可视化通过使用R软件的‘ade4’包的类间分析、主成分分析或主坐标分析构建图形来实现[31]。由于在样品中某些属或种的相对丰度很低，所以在分析之前，可以用适当的阈值来过滤相对丰度较低的物种以降低噪声。

根据物种平均相对丰度大于0.1%过滤后的肠型分类可视化，如图 1所示。

生物信息学工具已被广泛用于肠型的研究，相比与血型一样的离散类型，对肠型分布更好的理解应该是在物种特征的高维空间中的密集分布(就如http://enterotype.embl.de/所描述的：Enterotypes are not discrete types，like blood types，but rather densely populated regions in a higher dimensional space of microbiome features)，因此将肠型作为根据肠道微生物组对个体进行分层的标准或许将更为合适。

最新的《自然∙微生物学》杂志报道了学界对于肠型的最新认识，调解了之前关于肠型的一些争论[32]。首先，肯定了肠道菌群大数据在人群中成层现象(population stratification)的客观存在性[32]；其次，对于之前出现或2种肠型、或3种肠型、或4种肠型的争论进行了系统分析，发现随着大数据的积累，肠型必然是连续的而非离散的，并且研究了各种不同方法对结果的影响，发现造成差异的原因是由于不同方法抽取不同维度进行观测，观测的其实都是复杂肠道菌群大数据整体的某个子集的不同切面[32]。然而不论如何观测，以后更多数据积累发现新类型，都无法否认目前人群中常见的几种肠型，如Bacteroides型、Prevotella型、Firmicutes型等的客观存在性。

在宏基因组学中，分类多样性通常表现为阿尔法(α)多样性和贝塔(β)多样性，这是生态学研究中常用的术语。

阿尔法多样性：单个样品中的分类群数量和这些群体的平均丰度称为阿尔法多样性，可以使用统计学方法确定，如Shannon-Wiener指数[33]、系统发育多样性全树分析[34]和Simpson指数[35]。通常，这些指数是逐步计算和递归计算的，以确定采样数据量是否具有足够的代表性以形成稀疏曲线。

贝塔多样性：阿尔法多样性描述了单个样本的特性，而贝塔多样性描述了样本之间的相似性。对于贝塔多样性分析，样品之间的距离/相似性，往往使用例如Bray-Curtis测量[36-37]方法和加权/未加权的UniFrac[1, 24, 38-40]方法。样本之间的相似性可以使用精准算法在二维图中可视化[41]。样本组成可以看作是一个向量，因此可以看作是存在于多维坐标系中的一个点。协调系统目录的简单和转换表明样本之间的距离，这是主成分分析的精确性质。

因为二代测序技术得到的reads仍旧较短，且其中的微生物组装复杂，因而基于组装得到相对更长的contigs进行的下游分析将更有效。由于在宏基因组中超过90%的微生物是未知的[42]，不依赖于参考序列的从头组装就显得很有必要了。现有的宏基因组组装算法有很多，绝大部分都是基于de Brujin图构建的方法，例如，MetaVelvet[43]和Meta-IBDA[44]通过reads构建de Brujin图，并根据de Brujin图的特性识别代表基因组特异性的子图(sub-graphs)；Genovo[45]根据组装概率模型识别最大可能的序列重构，从而得到具有最大似然值的contigs集合。Vollmers[46]对这些组装软件进行了评测，并提出了最优组装软件的选择依赖于样品的类型、可获得的计算资源，以及更重要的是研究目的。

通过基因预测可以识别宏基因组中的基因编码序列，从而了解宏基因组样品中微生物的基因功能。常用的基因预测软件有MetaGeneMark[47]、Glimmer-MG[48]、MetaGene[49]、Orphelia[50]、FragGeneScan[51]和MetaGun[52]。Trimble等[53]通过统计模拟数据对这些预测方法进行了比较，发现这些方法的性能会根据reads的性质(例如测序错误率和reads长度等)而有所变化，并且不同的方法将根据不同性质的阈值产生最优准确性，因此建议研究者需要基于自身研究数据的特性来选择合适的方法以及对应的参数。

当对多个宏基因组样品进行比较分析的时候，我们需要先根据每个样品预测得到的基因构建一个去冗余的基因集。我们通常使用cd-hit[54]进行去冗余基因组的构建，首先，将预测的基因序列进行两两比对，然后根据比对结果将能比对到另外一条较长序列(相似度大于90%且覆盖度大于95%)的基因序列列为冗余序列去除，最终得到一个去冗余的基因集。考虑到由于测序数据有限会导致一些低丰度微生物不能被检测到，所以通过整合已经公布的基因集，例如MetaHIT基因集[55]、HMP基因集合[56]和LC基因集[14]等，将有利于优化自身数据分析得到的基因集合。

差异基因分析(或差异物种分析)就是在基因丰度表(或物种丰度表)中寻找针对不同组之间有显著性差异的一种分析，一般使用针对非正态分布的非参检验方法，具体而言，有Wilcoxon秩和检验、Mann-Whitney U检验和Wilcoxon signed-rank检验等。其中，对于2个独立样本(两组来源于不同总体)的差异分析方法可以使用Wilcoxon秩和检验与Mann-Whitney U检验，Wilcoxon秩和检验与Mann-Whitney U检验并没有实质性的差别。对于2个配对样本(两组来源于同一个总体)，应该使用Wilcoxon signed-rank检验来代替Mann-Whitney U检验(或Wilcoxon秩和检验)。当样本量很小时，我们采用Wilcoxon rank-sum permutation检验，在这种情况下，P值是基于置换(排列)的估计分布来代替一个固定的正态分布来计算的。对于3组或以上的差异基因检验方法有Kruskal-Wallis检验、Friendman检验等。

通常的方法是，在每组中过滤掉平均丰度低于10-8的基因，然后用差异检验方法来鉴定与实验变量相关的基因。因为做差异基因的筛选，有多个基因要筛选，需要做多重假设检验，通过检验所获得的P值直接与设定的显著性水平进行比较会有很大的误差，所以我们通过多重假设检验，将所获得的P值通过FDR (false discovery rate)来进行校正，并通常设置校正后的阈值(Q值)为0.05。

物种丰度是指物种数量的多少。肠道的物种丰度表为肠型分析、肠道多样性计算以及寻找疾病的微生物标记物奠定了基础。

目前，有两种较为常用的计算物种多样性的方法。一种是将测序所得到的宏基因组序列与全基因组序列进行比对，经计算后得到物种丰度表。该计算方法将唯一比对到的序列丰度与多重比对到的序列丰度求和得到物种丰度[14]，按照公式(1-4)计算。

式中，Ab(S)表示S物种的物种丰度，U表示唯一比对到S物种上的序列数量，M表示比对到多个物种的序列数量，l表示参考基因组的长度，N表示物种的数量，若物种S的Ab(U)为0，那么，Co也为0。在计算完成所有物种丰度之后，将物种S的物种丰度除以所有物种丰度之和，就可以得到某样本的物种丰度表。

另一种方法是用MetaPhlAn2[24]或mOTU[39]软件将测序得到的宏基因组序列比对到标记序列数据库上，计算后就可以得到物种丰度表。mOTU和MetaPhlAn2比对的数据库有所不同，MetaPhlAn2是将序列比对到一个已构建好的肠道数据库上，该数据库包含了13500多个细菌古菌，3500个病毒和110个真菌的序列标记物，将比对到物种S的标记序列上的序列数除以比对到标记序列上的序列总数，就得到了物种S的物种丰度，计算所有的物种丰度就构成了物种丰度表，mOTU计算物种丰度需要依赖于40个通用单拷贝系统发育标记基因序列，这些基因在绝大多数已知生物的基因组中以单拷贝形式出现，将测序序列比对到这些标记基因可以生成物种丰度表，这40个标记基因都是从3445个完整的细菌基因组和来自263个公开的宏基因组的未知物种中提取得到的。

基因丰度表是寻找差异基因、研究疾病与健康差异、挖掘深层生物机理的基础。在计算基因丰度之前，需要先构建非冗余基因集，构建好基因集之后就可以将原始序列用SOAPaligner[20]、BLAST[57]、BLAT[58]、BWA[21]或Bowtie2[32]等比对软件比对到基因集上，通过计算可以得到基因丰度。计算方法与计算物种丰度类似，假设现有基因G，基因G的丰度就等于唯一比对到的序列丰度与多重比对到的序列丰度之和。按照公式(5-8)计算。

式中，Ab(G)表示基因G的基因丰度，U表示唯一比对到某基因上的序列数量，M表示比对到多个基因的序列数量，l表示基因序列的长度，N表示基因集中的数量，若基因G的Ab(U)为0，那么，Co也为0。在计算完成所有基因丰度之后，将基因G的基因丰度除以所有的基因丰度之和，就得到某样本的基因丰度表。

目前，功能分析的应用非常普遍，以代谢信号通路为例，与其相关的文献逐年增加[59-61]，所以宏基因组功能分析将会是一个非常普遍而且重要的分析内容。宏基因组功能分析可以获得基因的注释信息，了解相关的代谢信号通路，有助于研究者从系统水平去了解基因的生物学功能，进一步揭示宿主与菌群之间的关系[60]。

宏基因组功能分析的主要参考数据库包括京都基因与基因组百科全书(KEGG)数据库[62]、EggNOG数据库[63]、GO数据库[64]、耐药基因(ARDB)数据库[65]、碳水化合物酶(CAZY)数据库[66]、群体感应(QSDB)数据库、病原与宿主互作(PHI)数据库[67]、病原菌毒力因子(VFDB)数据库[68]、转运蛋白分类(TCDB)数据库[69]等。

代谢信号通路(pathway)的主要参考数据库有KEGG和EggNOG。KEGG是基因组破译方面的数据库，是一个整合了基因组、化学和系统功能信息的数据库[62]，包括完整和部分测序的基因组序列，更高级的功能信息存储在PATHWAY数据库里，包括图解的细胞生化过程如代谢、膜转运、信号传递、细胞周期，还包括同系保守的子通路等信息[70]。KEGG通过图形方式呈现众多的代谢途径以及各途径之间的关系，这样可以使研究者能够对其所要研究的代谢途径有一个直观全面的了解。

EggNOG数据库是利用Smith-Waterman比对算法对构建的基因直系同源类群进行了功能描述和功能分类的注释，从另一个角度提供基因功能分析。

与KEGG和EggNOG存在一定差异的是，GO数据库主要用于研究基因功能，而非研究基因和蛋白功能。

其他数据库更多专注于具有特殊功能的基因，如CAZy数据库包含了碳水化合物酶类的物种来源、酶功能EC分类、基因序列、蛋白质序列及其结构等信息。

ARDB数据库包含细菌病原菌的多种抗性基因数据，通过该数据库的注释可以找到耐药性相关的基因名称以及所耐受的抗生素的种类。

PHI数据库包含了被实验证实对病原体与宿主相互作用有影响的基因，主要来源于真菌、卵菌和细菌病原体，宿主包括动物、植物和真菌。

TCDB数据库是一个用于运输蛋白质研究的参考数据，包含了各种生物体运输系统有关蛋白的序列、结构、功能以及进化信息。

VFDB数据库包含细菌病原体的毒力基因序列信息，被广泛应用于毒力因子基因鉴定。

通过序列比对，将基因集比对到相应的参考数据库，获得基因的注释信息，对注释信息分析，了解基因参与的生物学功能。通常使用BLAST[57]软件进行基因集注释的序列比对。

基因之间的关联分析与β多样性分析相似，通常采用斯皮尔曼等级相关系数和皮尔逊相关系数计算相关性。当2个物种之间的相关系数高时，它可能表明这些物种是相关的。同样地，当宏基因组分析中2个基因之间的相关系数高时，可以表明这些基因来自相同的物种。

在宏基因组分析中，可以根据它们的相关系数将基因分组在一起。当选择合适的分组标准时，来自相同物种的基因很可能被聚集在一起，这样的基因群被称为MGS，并且这些基因群中的基因可以与目前可用的参考基因组比对以验证“物种”检测的分组方法，一些未比对上的基因群也许来源于未知物种。从MGS数据中可以鉴定和提取属于特定物种的基因。因此，微生物基因组可以使用宏基因组数据进行组装，MGS可以提供可用于计算分类组成的潜在参考基因组。

CAG方法是2014年在《自然∙生物技术》杂志上报道出来的人体肠道宏基因组数据分析方法[13]。这种方法的思想和MGS的思想有一定相同之处，都利用了这样一种思想或理论：在同一个株或种中的基因，这些基因的丰度在群体中具有丰度一致性。为了解决使用MGS方法时仅关注基因标记的局限性，canopy算法被用来鉴定共丰富基因簇(CAG)，CAG分析基于基因丰度鉴定物种，并且每个CAG可以被看作是部分微生物或完整微生物。

当基因丰度表包含足够数量的样本时，可以应用canopy算法来识别CAG。canopy算法使用Pearson相关系数和Spearman相关系数作为阈值，在对只有1个基因的canopies进行第一轮聚类和过滤之后，根据簇的平均丰度，使用Pearson相关系数作为阈值再次被应用到canopy算法中，第二轮簇可以包含重叠基因。因此，对于出现在多于1个簇中的基因，基因及其相关簇之间的距离能够被确定，对于每个重叠基因，选择最近的簇。最后，选择含有超过700个基因的簇最可能是含有细菌基因组部分的CAG。

基于CAG分析的结果，可以在人肠道微生物群中鉴定出可能显著影响宿主的未分类或难以理解的微生物(也称为生物暗物质)。在过去的研究中，生物暗物质一般被忽略，因此，许多疾病相关的肠道微生物可能未被检测到。CAG中的重叠群和基因可以产生无法从当前的细菌基因组数据库获得的大量新信息。

医学分类器是通过特征工程中的特征选择的方法挑选出若干不同级别的bio-markers，并根据bio-markers来构建的一种分类模型，它采用分类评价指标来评估模型效果的好坏，目的是用于识别疾病与健康人群。

医学分类器是来自机器学习中的分类算法，包括Logistictic回归、Support Vector Machine (SVM)、Random Forest、人工神经网络、朴素贝叶斯等。除了分类算法，特征选择对分类性能也会有直接的影响，特征选择就是在我们获得的给定特征集中选择出与分类相关的特征子集的过程。mRMR (maximal relevance and minimal redundancy)是一种常用的特征选择的评估准则。评估方法一般包括Filter、Wrapper、Hybrid、Embeded，其中Wrapper是较为传统和常用的方法，它是直接把最终将要使用的学习器的性能作为特征子集的评价准则。通过分类算法得到的分类器的性能最终是由分类评价指标来确定，分类评价指标是用来评价分类器分类性能好坏的一种指标，常见的评价指标有正确率、错误率、灵敏度、特异度、精度、Matthews correlation coefficient (MCC)、Receiver operating characteristic (ROC)。其中MCC和ROC不受分组样本数量影响；ROC是根据一系列不同的二分类方式，以真阳性率(灵敏度)为纵坐标，假阳性率(1-特异度)为横坐标绘制的曲线，直观描绘了分类器在TP和FP间的trade-off，AUC (area under curve)是ROC曲线下方面积之和，它以数值的形式来评估分类器的好坏，AUC的取值范围在(0.5, 1.0)范围内，AUC的值越大，分类的性能就越好。若分类评价指标最终显示分类器性能不佳，则需重新构建分类器直至分类性能达到稳定为止。

医学分类器的引入对于构建用于临床诊断的基于肠道微生物的非接触检测方法具有重要意义。

随着生物信息行业的迅猛发展，越来越多的组装软件被开发出来，然而还有许多问题有待解决。从人体肠道宏基因组测序数据进行全基因组组装，为下一步基因预测、基因集构建、基因丰度统计提供了重要的基础数据。肠道菌群非常复杂，许多不同菌株来自相同物种，以及许多不同物种来自相同属，这些高度相似DNA序列会混淆组装结果，导致不正确的基因组序列出现，因此组装算法还有很大的提升空间。二代测序的短读长限制了contig的长度，随着三代测序技术成熟，基于三代PacBio的SMRT单分子测序在宏基因组方面也已开展了一些应用，PacBio产生的序列平均长度在5-50 kb，其中50%左右的序列长度大于14 kb，宏基因组拼接组装效果有了很大程度的提高，能组装出更长的contig，提高了短reads的使用率，增加了完整基因的数量，这表明宏基因组的组装又上升了一个台阶。随着测序读长不断增加，测序质量的不断提高，将三代测序数据应用到宏基因组中将是未来研究的一大方向。

随着技术的发展，亚种或株水平的研究必然对人体肠道宏基因组领域产生深刻影响。大量(60%-70%)的肠道微生物目前还无法从现有国际数据库中找到参考基因组，虽然人体肠道宏基因组领域有CAG等方法对亚种或株通过生物信息学方法进行分析，但对这个方向的根本性突破还是需要准确度更高、读长更长的测序技术如三代测序技术去推动。当测序准确度提高、测序读长接近微生物自身基因组长度，大量未知的肠道微生物基因组将通过最简单、直接的方式整体性地被呈现出来，所谓的生物暗物质将向我们揭开神秘面纱。

2017年，美国克利夫兰医学中心将“利用微生物组预防、诊断和治疗疾病”列为2017年十大医疗科技创新之首。微生物组技术本身就是医疗前沿创新技术，随着近年来各项新技术的发展，微生物组与多个热门领域新技术相结合，更体现出了微生物组技术的活力和未来的广阔应用价值。从2017年到2018年，在肿瘤免疫治疗领域，高水平学术期刊《科学》连续报道了肠道菌群显著影响PD-1 (程序性死亡受体1)药物治疗的效果[71-72]，并且文章通过宏基因组技术找到与PD-1免疫疗法相关的10种微生物。在未来肿瘤免疫治疗用药方案选择前，通过肠道菌群检测及宏基因组分析或是一种常规的手段。而另一个炙手可热的领域CRISPR/Cas9基因编辑技术，在未来改造菌群对抗日益严峻的抗生素滥用、细菌耐药性乃至通过菌群改造治疗疾病，也有与宏基因组技术相结合应用的广阔空间。

感谢吴春燕、范芳芳、宋文峰、李心如、马圣、熊潇、马玮在资料收集和文章撰写中提供了帮助。