人类微生物组研究设计、样本采集和生物信息分析指南

A guide to human microbiome research: study design, sample collection, and bioinformatics analysis

Chinese Medical Journal [IF: 1.585]

DOI: https://doi.org/10.1097/CM9.0000000000000871

Review: 2020-6-26

钱旭波¹, 陈同², 徐益萍¹, 陈雷³, 孙馥香⁴, 卢美萍¹, 刘永鑫^5,6

1. 浙江大学医学院附属儿童医院风湿、免疫和变态反应科
2. 中国中医科学院中药资源中心
3. 首都医科大学附属复兴医院
4. 易汉博基因科技（北京）有限公司
5. 中国科学院遗传与发育生物学研究所
6. 中国科学院大学，生物互作卓越创新中心

钱旭波和陈同为共同第一作者

通讯作者：卢美萍，浙江大学医学院附属儿童医院风湿、免疫和变态反应科，中国浙江杭州竹竿巷57号，邮编：310003，邮箱：meipinglu@zju.edu.cn

摘要

这篇综述的目的是为医学研究人员，特别是那些没有生物信息学背景的研究者提供简单易懂的微生物组学知识，包括研究中常用的概念、技术和分析方法等。首先，我们介绍了基本概念，例如微生物群（microbiota）、微生物组（microbiome）和宏基因组（metagenome）等。然后，我们讨论了研究设计方案、样本量计算方法以及提高研究可靠性的方法。我们特别强调了阳性和阴性对照的重要性。接下来，我们讨论了微生物组研究中常用的统计分析方法，重点关注多重比较的问题以及组间β多样性分析的方法。最后，我们介绍了生物信息学分析的具体流程。总之，严谨的研究设计是获得有意义结果的关键步骤，而适当的统计方法对于准确解释微生物组数据很重要。通过阅读这篇文章，研究者能获得研究设计、样本采集和生物信息分析等全方位的微生物组学知识。

关键词：微生物组、研究设计、统计分析、样本量、生物信息分析、分析流程

1. 前言

随着测序技术和数据分析方法的发展，近几年医学微生物组研究领域出现了一些令人瞩目的成果[1-3]，比如微生物组与代谢性疾病[4-6]、消化系统疾病[7-10]和心血管系统疾病[11]之间的关系日益明确。这些发展和发现增加了医生在微生物组研究方面的兴趣，进而也涌现出了大量有价值的论文[12]。另外，随着QIIME 2[13]和多组学方法[1, 9]等先进技术和分析流程的出现，微生物组分析方法也不断进步。然而，理解和掌握这些技术和分析流程并非易事，特别对于医生来说更是如此。

本文的目的是为研究者，特别是那些没有生物信息学背景的医生提供易懂的微生物组学知识，这些知识包括详细的微生物组学基本概念、科研设计方法、样本采集和保存方法、统计分析方法以及生物信息分析方法。我们希望医生们通过阅读此文能够快速掌握以上知识和方法，进而有效地挖掘数据背后的生物学意义。

2. 基本概念

2.1 Microbiota、Microbiome等术语

Microbiota(微生物群/微生物组)是指定植在人体特定部位的微生物，包括细菌、古菌、病毒、真菌和原生动物[14, 15]。在医学研究中，如果测序技术采用的是16S rRNA基因（又称为rDNA），则microbiota是指细菌和古菌。Microbiome是指整个微生境，包括微生物、基因组和周围环境[14, 15]。不过，microbiota和microbiome有时存在混用情况。我们建议，如果你的研究仅涉及微生物本身，则应该使用microbiota，否则应该使用microbiome（图1）。例如，如果研究者想探索肠道短链脂肪酸与微生物的关系，使用microbiome更合适。宏基因组（metagenome）是指微生物基因组的集合[14]，一般用鸟枪法宏基因组测序获得，宏基因组学则是研究宏基因组的学科[12, 14]。病毒组（virome）指人体内或表面的病毒集合，包括内源性逆转录病毒、真核生物病毒和噬菌体[16]。研究病毒组的学科就是病毒组学。作者注：Microbiota国内有些学者翻译为“微生物群”，microbiome翻译为“微生物组”。不过中文文献用“微生物组”或“××菌群”即可，多数情况下不需要区分是microbiota或microbiome。

图1：微生物组、微生物群、宏基因组和16S rDNA的概念。

（A）微生物组（microbiome）的概念不仅涵盖微生物，而且涵盖周围的环境条件。微生物群（microbiota）仅指微生物本身。（B）宏基因组是指微生物的所有基因组，而16S rDNA仅涵盖基因组的一部分。（C）α多样性衡量样本中的多样性，而β多样性比较样本之间的物种差异。

2.2 细菌层级分类

细菌分类最常用的层级为门、纲、目、科、属、种、株。例如，临床上十分常见的大肠埃希菌的层级分类见表1。

表 1: 大肠埃希菌细菌层级分类

2.3 操作分类单元和扩增子序列变异

操作分类单元（operational taxonomic units，OTUs）的构建对于标记基因（扩增子）数据分析非常重要[17]。OTU是指一组高度相似的序列，通常将具有97%相似性的一组序列归为一个OTU[18, 19]。不过，这种OTU的方法有显著的缺点，它人为地设置一个相似性阈值，漏掉了细微的和真正的生物学序列差异[20]。最近开发的扩增子序列变异（amplicon sequence variants, ASVs）方法可以解决这些问题，它使用序列变异信息将序列数据解析为准确的序列特征。ASV具有单核苷酸分辨率，并且具有比OTU相似或更好的敏感性和特异性[20]。注意，OTU或ASV不等于物种，一个OTU / ASV可能包括多个物种，反之亦然[21]。

2.4 α-多样性

α-多样性是指样本内的多样性，常见的样本有粪便，唾液或支气管肺泡灌洗液等[15]。医学研究中经常使用3种α多样性指数：Chao 1指数、香农指数和辛普森指数。Chao 1指数主要反映物种数量（richness），它计算时考虑以下三个因素：物种数量、单条序列数量和双条序列数量[22]。这意味着它不能反映微生物组的丰度（abundance）。香农指数结合了丰度和均匀度信息[23]，它赋予稀有物种更多的权重[22]，这意味着当稀有物种的数量增加时，它的值会更大。香农指数的值通常不超过5.0；它的值越高，α多样性就越丰富[22]。辛普森指数也整合了丰度和均匀度，不过与香农指数比较，计算时它对常见物种有更大权重。它的值介于0-1之间，这个值越大，α多样性越丰富[22]。在以上指数中，richness是指一个样本中物种的数量[17, 24]，而abundance（丰度）指物种的原始序列读数[24]。如果原始序列读数被转换成百分比后，它就称为相对丰度。

2.5 β-多样性

β-多样性是指样本或组间的微生物组差异，通常用于了解两组微生物组组成的差异是否显著。在这里，我们关注两个常用的β多样性指数：Bray-Curtis相异性和UniFrac距离。Bray-Curtis相异性是一种用于量化两个样本或组间的物种组成差异的指标，其值的范围是0到1，其中0表示两个样本或组间具有相同物种，而1则表示它们不共享任何物种[25]。此外，它在计算时给予常见物种更大的权重[23]。请注意，Bray-Curtis相异性不是真正的距离度量指标，因此用“Bray-Curtis相异性”的叫法比“Bray-Curtis距离”更恰当[22]。

UniFrac距离可以不加权，也可以加权，它基于系统发育距离估算微生物组样本或组间的差异[26]。未加权的UniFrac距离只考虑了物种是否存在，它对于检测稀有物种的数量变化很敏感，但是在计算中忽略了丰度信息[27]。加权UniFrac距离计算时纳入了丰度信息[28]，并减少了稀有物种的权重[29]。

2.6 排序

排序用于探索数据结构，由降维后的正交轴图形表示。排序图是可视化β多样性的有效方法。排序可以分为2大类：非约束排序和约束排序[30-32]。如果图形上的点不受环境因素（样本元数据）的约束，这种排序叫做非约束排序，否则叫约束排序[32]。常用的非约束标准包括主成分分析（principal component analysis, PCA）、对应分析（correspondence analysis, CA）、主坐标分析（principal coordinate analysis, PCoA）和非度量多维标度（non-metric multidimensional scaling, NMDS）[30, 32]。常用的约束排序有冗余分析（redundancy analysis, RDA）和典范对应分析（canonical correspondence analysis, CCA）[31, 32]。

微生物组信息是高维数据。PCA通过将数据以几何方式投影到较少的维度上来简化复杂性，它在计算中使用欧几里得（Euclidean）距离[30]。通常情况下它并不适用于物种丰度数据的分析，因为PCA分析的数据必须是线性的[30]。但是如果物种数据经过Hellinger转换，则PCA可以用于物种数据分析[30]。相反，CA适合于物种丰度数据分析，而且无需预先转换数据。在CA分析中，所有样本均使用Pearson卡方距离进行排序[30]。但是请注意，稀有物种可能会对CA分析产生过大影响[33]。如果研究人员希望基于相异性指标来对样本或特征进行排序，那么PCoA是一个不错的选择。在微生物组研究中，PCoA分析最常使用Bray-Curtis相异性和UniFrac距离。NMDS用于表示排序图中样本的相对位置。与PCoA相似，NMDS分析可以使用任何距离或相异矩阵。参考文献[30]详细介绍了PCoA和NMDS之间的差异，在大多数情况下PCoA比较常用。

RDA是一种结合了PCA和回归的约束排序，它的响应矩阵是微生物组数据，解释矩阵是临床指标（样本元数据）。RDA对于显示微生物组数据是否受临床指标影响很有用。但是请注意，由于PCA计算过程要求响应矩阵的数据结构必须是线性的，因此可能需要对数据进行预转换。最后，CCA其实就是CA的约束版本，它具有CA的基本特性和缺点[31]。

3. 研究设计

3.1 研究设计方案

严谨的研究设计对于获得准确而有意义的结果很重要。医学微生物组研究中最常使用的研究方法包括横断面研究、病例对照研究、纵向研究和随机对照试验（randomized controlled trial, RCT）。前3种是不应用干预因素的观察性研究，而最后一个是典型的实验性研究。

横断面研究分为描述性横断面研究和分析性横断面研究[34]。前者仅是描述性的，主要用于调查一个或多个人群中的微生物组成，而后者则用于探讨微生物组与健康结果之间的关联。但是，微生物组与健康结果之间的关联可能源于混杂因素，例如性别[35]、年龄[36]、体重指数（body mass index, BMI）[37]、饮食[5, 38]、季节[39]和药物治疗[40, 41]。此外，横断面研究时，微生物组和结果是同时测量的，因此很难确定它们之间的因果关系。通常，横断面研究仅用于探索微生物组的基本特征，并且可以作为后续研究的初步实验。

在大多数情况下，微生物组被视为暴露（exposure），疾病被视为结局（outcome）。在这些假设下，传统的病例对照研究很少用于微生物组研究，因为以前的暴露（微生物组）信息很难获得。但是，如果暴露和结局对调，则可以使用病例对照研究设计方案。

同样，在上述假设下进行前瞻性队列研究也很困难，因为很难知道哪些微生物是潜在的暴露。而且，定义可用作暴露或非暴露因素的特定微生物组并非易事，因此难以将研究对象确定为暴露或非暴露个体。在实践中，有或没有疾病的个体通常归入研究组或对照组，然后在不同时间点前瞻性地收集含有微生物组的样本[17]。也就是说，前瞻性队列研究中的研究对象通常根据临床结局而不是特定的微生物组模式进行分组。

RCT或其他实验研究的目的是评估干预措施的有效性。干预措施可以是药物或微生物组。例如，粪菌移植研究中的干预措施是微生物群[42, 43]。

值得注意的是，对照组的选择应恰当。以上这些研究设计中应注意匹配混杂因素，这部分内容将在下面讨论。有时对照的选择很困难，尤其是在临床研究中干预措施是微生物群本身的情况下。在这种情况下，如果其他研究设计不合适，那么进行有对照组的前后自身对照试验（controlled before-after trial）或历史对照试验将是一个不错的选择[44]。

3.2 定义纳入和排除标准

定义确切的纳入标准和排除标准可以使组间更好地匹配，并且有利于控制混淆因素，比如年龄[36, 45]、性别[35]、BMI[46]、饮食[47]、季节因素[39]、药物治疗[40, 41]、种族[48]、地理区域[45]和共存疾病等[7]。年龄可显著影响微生物组，对于那些小于16岁的人更是如此[36, 45]。因此，对于涉及儿童的研究，年龄必须很好地匹配。饮食是另一个对微生物组改变有影响的因素，所以也要进行匹配[47]。为了增加组间的可比性，地理区域因素在研究设计时也需要考虑在内[45]。由于药物治疗对于微生物组有显著影响，所以入组前数月内接受过药物治疗的患者应该排除在外[41, 49]，这里讲的数月通常指入组前3~6个月[49]。

3.3 微生物组研究的样本量和检验效能计算

在进行研究设计时估计样本量大小非常重要。适当的样本量可使微生物组研究识别出组间的差异，并节省资源和时间。但是，样本量和检验效能计算对于研究者来说仍然是一个挑战[50]。微生物组研究中最常用的样本量和效能计算方法可以用t检验、方差分析、χ2检验和Dirichlet多项式模型[51]。以t检验为例，分3个步骤确定样本大小和效能计算。首先，通过初步实验获得少量扩增子数据。其次，使用R包vegan计算出每个样品的香农指数[52]。最后一步是使用R软件包pwr中的power.t.test()函数计算样本量和效能。当研究者仅关注两组之间物种多样性的差异时，可使用t检验计算样本量和效能。在参考文献[51]中有样本量和效能计算的详细介绍。

3.4 阴性和阳性对照的重要性

微生物组研究的结果可能会受到多种因素的影响，例如DNA提取试剂盒、采样方法、污染和测序方法等[53]，不过可以通过使用阴性和阳性对照来减少这些影响。不幸的是，以前的研究中只有30％报告使用了阴性对照，只有10％报告使用了阳性对照[53]。使用对照对于准确认识微生物组非常重要，尤其是当样本的微生物含量较低时。以前的研究发现，过去被认为是无菌的标本（例如胎盘和关节液）可能会被微生物定植[54]。但是，这些阳性的结果可能是由其他因素导致的，例如污染。有趣的是，这些低生物含量标本在采用阴性和／或阳性对照后已被证明是无菌的[55]。因此，我们建议当样本为低生物含量样本（例如血液、羊水、脑脊液、关节液和胎盘等）时，应考虑使用阴性和阳性对照。值得注意的是，阴性对照和阳性对照在病毒学研究中也很重要，因为病毒和细菌通常是同时进行检测的[16]。此外，R包decontam可用于鉴定和去除扩增子和宏基因组学数据中的污染物序列[56]。

3.5 测序方法的选择

微生物组研究中使用的测序方法包括扩增子测序、宏基因组测序和宏转录组测序。扩增子测序包括适用于细菌和古菌的16S rDNA测序以及适用于真菌的内部转录间隔区（internal transcribed spacer, ITS）测序。每种测序方法的优缺点在这两篇参考文献中有详细讨论[17, 57]。简而言之，扩增子测序很便宜，可应用于受宿主DNA污染的低生物含量标本，但一般仅能注释到“属”层级，并且易受某些固有偏倚来源的影响，例如PCR循环数[58]。宏基因组测序方法对样品中存在的所有DNA进行测序，包括细菌、病毒、真核生物和宿主的DNA。它不仅将其分类学分辨率扩展到“种”或“株”的水平，而且还提供了潜在功能信息[17]。但是，扩增子和宏基因组测序方法都无法区分死微生物或活微生物[17]。转录组测序仅产生群落的活跃功能信息。鉴于这些测序方法的优缺点不同，建议将多种测序方法整合在一起以优化研究设计。简而言之，测序方法的选择主要取决于实验成本和样本质量。扩增子测序通常用于获得微生物群落的概况[59]，并且通常适用于大规模研究[6, 60]。如果您有足够的项目资金，并且想要获得菌株水平的分辨率和潜在功能，甚至想要恢复整个基因组，宏基因组测序是一种首选方法[61-65]。

3.6 提高研究可靠性的方法

简单的横断面研究在微生物组研究中的意义有限。在本小节中我们讨论了提高研究可靠性的方法。首先，首选纵向研究或RCT研究，而不是横断面研究或病例对照研究[17, 66]。其次，应计算样本量[51]。第三，混淆因素应匹配，元数据（即各种临床指标等信息）应仔细收集。第四，应详细定义纳入和排除标准。例如，幼年特发性关节炎有几种亚型，每种亚型可能代表不同的疾病[67]。研究者应确定患者组中是否包括所有亚型。第五，最好考虑使用阴性和/或阳性对照[68]。第六，整合其他组学方法，例如代谢组学、转录组学和蛋白质组学，这对于全面了解微生物群落的结构和功能至关重要[17]。因此，应考虑获取微生物群落代谢物概况和／或其他多组学数据。目前，仅探索微生物群落结构的研究不被视为论证效率强的研究设计[17]。最后，建议在动物模型中验证从临床试验获得的初步结果。

表2列出了设计临床微生物组研究需要考虑的因素，图2展示了典型的工作流程。实验研究需要考虑的因素见参考文献[49]。

表 2: 临床微生物组研究设计需要考虑的要素核对表

RCT：随机对照试验

4. 样本类型、保存和储藏

4.1 样本类型

人类微生物组研究的样本类型包括粪便、结肠灌洗液和腔内刷等（表2）。样本类型的选择取决于感兴趣的研究假设。例如，粪便样本易于收集，可用于大规模和纵向研究。另一方面，活检样本对于探索微生物群与宿主之间的相互作用更有用[69]。注意，在一项研究中应该固定采样位置，因为人体的不同部位定植着不同的微生物群[70, 71]。

图2：人类微生物组研究的典型流程。

4.2 保存和储藏

样品保存和储藏的方法应适合实验方法和样品类型。最通用的方法是直接冷冻样品，它可用于各种测序和实验方法，例如扩增子、宏基因组、转录组测序和代谢组学测定。建议将样品收集后15分钟内保存在-20℃下[72, 73]，然后在收集24小时内用干冰转移到-80℃冰箱中储藏。不过样本通常是在家里而不是在医院收集的，在这种情况下可以使用保存液。保存液中保存的样本可以在环境温度下保存一周以上[74]。请注意，样品的保存和储藏方法应一致，以最大程度地减少潜在的混淆因素干扰。

5. 微生物组研究中的统计分析方法

医学研究者通常熟悉单变量统计方法，例如t检验、方差分析、χ²检验和秩和检验。因此，我们在这里仅讨论与多重比较和其他多元统计方法有关的问题。我们首先讨论多重比较会遇到的问题及其解决方案，包括P值调整和使用错误发现率（FDR）。然后，我们讨论其他多元统计方法，例如置换多元方差分析（permutational multivariate analysis of variance, PERMANOVA）和Mantel检验。

5.1 多重比较的问题及解决方法

由于微生物组数据是高维的，因此多重比较经常在微生物组研究中使用。例如，特征表（feature table）具有成百上千个OTU或ASV，并且每个OTU或ASV都可以进行多次比较。医学研究者经常遇到的另一个例子可能更容易理解。假设一项研究分为3组，例如A组、B组和C组，而研究者想比较这3组之间的差异。在这种情况下就应调整P值，因为每个组都进行了2次比较，即A组与B组，A组与C组，B组与C组。如果有任何组或变量需要进行多次比较则必须进行P值调整，以便减少假阳性率[75]。

调整P值的经典方法是控制family-wise错误率，即Ⅰ类错误或α水平。Bonferroni是校正α水平最常用的方法。校正P值的计算非常容易：单个检验的α值除以检验次数。因此，对于上述具有3个检验次数的例子，调整后的P值为0.05 / 3 = 0.017，即只有P <0.017的检验结果才被认为是有意义的[75]。请注意，Bonferroni校正仅适用于多重比较次数较少的假设检验，否则会导致较高的假阴性率（图3）[75]。

解决多重比较问题的另一种方法是控制错误发现率（false discovery rate, FDR），它是I类错误或假阳性的数量与所有被拒绝的无效假设的预期比例。例如，如果100个阳性假设检验结果中有5个是错误发现，则FDR为5％。在微生物组研究中，通常使用“Benjamini-Hochberg（BH）校正的P值”而不是原始P值。校正后的P = 原始P m/i，其中m是检验次数，i是每个P值从小到大排序的序号[75]。如果校正后的P值小于你选择的所选FDR，则认为该检验是有统计学意义的。与Bonferroni方法相比，BH方法不那么保守（即校正强度不是很大），BH法通常用于微生物组特征的多重比较。Bonferroni和BH是最常用的P值校正方法[76]，这两种P值校正方法的校正强度见*图3所示。

图3：不同P值校正方法的校正强度

该图显示，Benjamini-Hochberg校正强度小于Bonferroni。随着原始P值的增加，Bonferroni校正法生成的校正后P值快速接近1.0。

5.2 PERMANOVA检验

有几种统计方法或模型可以用于组间β多样性比较，比如PERMANOVA、Mantel检验、相似性分析（ANOSIM）和多响应置换程序（multi-response permutation procedures, MRPP）。PERMANOVA最常用，并且被认为是以上检验方法中检验效能最大的一种[77]，它可通过R包vegan中的函数adonis()实现[52]。vegan包可计算4种常用相异性或距离度量：Bray-Curtis相异性、Jaccard距离以及加权和未加权UniFrac距离[29]。如果PERMANOVA检验的P值小于0.05，则表明不同组间的β多样性差异具有统计学意义；该检验的另一个输出结果是R²，它表示总方差可以用分组因素来解释的比例[29]。

5.3 Mantel检验

Mantel检验通常用于分析元数据矩阵和微生物组矩阵之间的关联[77]，它可使用R包vegan中的mantel()函数实现[52, 77]。该检验的输出至少2个主要统计量：P值和r。与其他类型的相关系数类似，r的值范围是-1 ~ +1[29]。例如，假设研究人员想知道元数据种的分组因素（例如吸烟状态）是否对肠道微生物组的组成产生影响。如果P＜0.05并且r＞0，这表明吸烟组和不吸烟组之间肠道微生物组的组成不同，元数据矩阵和微生物组矩阵呈正相关。

6. 生物信息分析

6.1 扩增子数据分析：从原始数据到物种分类表

有几种流行的软件或分析流程（pipeline）可用于扩增子数据分析，例如QIIME 2[13]、USEARCH[78]、VSEARCH[79]和mothur[80]。前两者具有许多优点，并已被许多研究者使用和推荐。每种软件或分析流程的优缺点已在我们先前的论文中详细描述[81]遗传：微生物组数据分析方法与应用 和 Protein Cell：扩增子和宏基因组数据分析实用指南。扩增子分析的主要步骤见图4A。我们通常从fastq格式的原始双端Illumina数据开始，最终输出是一个特征表，也称为OTU表或ASV表。

第一步是从原始数据中恢复纯净的扩增子序列，因为原始数据包括人造序列产物（artifact），例如引物和标签（barcode）。它包括3个主要过程：合并双端序列，通过标签拆分序列和去除引物。由于原始数据没有统一的标准格式，因此我们需要设计适合上述过程的分析流程。另外，我们也可以使用基因测序公司提供的纯净扩增子数据。图4B显示了用于恢复纯净扩增子序列的典型分析流程。

第二步是滤除低质量序列，以便减少背景“噪音”。

第三步是识别非冗余序列并且计数。高质量序列仍然有许多人造序列产物，例如错误序列和嵌合体。非冗余序列的计数是找出可靠序列的关键信息。

第四步是选择代表性序列（特征）。此步骤基于唯一序列，并通过将序列聚类成OTU或降噪生成ASV来实现[18, 82]。此步骤还包括de novo检测和去除嵌合体。

第五步是有参嵌合体检测，这是可选项[83]。通过将序列比对到数据库中，例如rRNA数据库SILVA[84]，可以进一步过滤特征序列。应当指出的是，该步骤可以降低假阳性率但易于导致假阴性结果。

最后，通过将纯净的扩增子数据与特征序列进行比较来生成特征表（图4A）。然后使用基于RDP[85]、SILVA或Greengenes[86]数据库的分类器实现特征序列的物种分类。此外，基于16S rRNA基因谱，使用PICRUSt[87, 88]、FAPROTAX[87, 89]和BugBase[90]等工具可实现功能预测。

6.2 宏基因组分析：从原始数据到物种和功能分类表

扩增子测序仅能获得微生物组部分的分类学信息，而且PCR过程很容易产生偏倚和嵌合体[83]。鸟枪宏基因组测序比扩增子测序提供更详细的基因组信息和更高的分类学分辨率[66]。与扩增子方法相比，宏基因组学分析更为复杂，但是它提供了更准确的物种分类、多维度的功能信息，甚至是末培养微生物的基因组草图。宏基因组分析流程如图4C所示。

第一步是预处理原始序列数据。原始数据包含低质量的污染序列以及与宿主相关序列。我们可以使用FastQC软件（http://www.bioinformatics.babraham.ac.uk/ projects/fastqc/）进行数据质量检查，然后使用KneadData流程进行质量控制[91]并去除宿主DNA[92]。有关更多KneadData的信息，请访问 http://huttenhower.sph.harvard. edu/kneaddata 。

第二步是使用基于序列的方法分析物种分类和功能代谢特征。人类微生物组具有高质量的基因集（gene catalog）和基因组[64, 65]，因此我们建议使用HUMAnN2[93]工具并采用基于序列的方法进行物种分类和代谢通路分析，该方法高效且易于操作。但是，这种方法只使用一小部分序列信息，而且分析结果受到已知数据库的限制[66]。

如果需要发现新物种或基因功能，则需要进行第三步。有几个好的软件工具可以用于将纯净序列组装为重叠群（contigs），例如MEGAHIT[94]和metaSPAdes[95]。 然后通过MetaProdigal[96]或Prokka[97]从长序列中预测基因。另外，其他软件工具也可以用于从短序列中预测编码基因，例如MetaGeneAnnotator[98]、MetaGeneMark[99]、Glimmer-MG[100]、MetaGUN[101]、FragGeneScan[102]和Orphelia[103]。为了减少重复基因，在分析多个样品或批次时需要使用CD-HIT构建非冗余基因集[104]。通过采用Bowtie 2[92]或Salmon[105]工具进行比对的方法可以计算基因丰度。目前至少有20个软件工具可用于宏基因组数据物种分类[106]。我们建议使用超快速分类器Kraken 2，它可以提供快速、准确和“种”级别的分类结果[107]。至于功能注释，许多研究人员都推荐使用DIAMOND[108]，它是一种快速、敏感的蛋白质比对工具[108]。每个数据库都提供了独特的功能视角，例如，京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes, KEGG）[109]、EggNOG（一个提供直系同源关系、功能注释和基因进化历史的数据库）[110]、碳水化合物活性酶数据库（Carbohydrate-Active enZYmes Database, CAZy）[111]、致病菌的毒力因子（Virulence Factors of Pathogenic Bacteria, VFDF）[112]和综合抗生素抗性数据库（Comprehensive Antibiotic Resistance Database, CARD）[113]。宏基因组通常包含100~1000个物种[64]，很难厘清彼此关系。分箱算法可以恢复无法培养的高丰度菌的基因组草图，并重建系统发育和代谢通路。

最后一步是使用metaWRAP[114]或DASTool[115]执行分箱流程（图4C）。这些软件工具有逐步操作教程，并且在其网站上提供了有关人类微生物组的一些样本数据集[81]。另外，几个集成的分析流程，例如MOCAT 2[116]、bioBakery[98]、IMP[117]和微生物组助手（Microbiome Helper）[118]，可以执行上述部分或全部分析步骤。你可以在微信公众号“宏基因组”中找到一些受欢迎软件的中文教程。

现在你已经获得了物种分类和功能信息文件。通过STAMP或LEfSe可以轻松找到你感兴趣的生物标记[119, 120]。使用R语言或ImageGP（http://www.ehbio.com/ ImageGP）可以将所有结果可视化。

图4：人类微生物组研究的生物信息学分析流程

（A）扩增子数据分析的主要步骤。（B）扩增子数据预处理的典型流程图：从原始的双端序列到纯净的扩增子。（C）宏基因组测序数据的分析流程。（a）预处理。它涉及删除低质量序列、接头和宿主序列。输出文件是纯净序列。（b）基于序列的分析。它将序列与数据库比对来推断物种分类和代谢特征。（c）基于组装的分析。它将短序列组装为长序列，预测基因，构建非冗余基因集，并与数据库比对进行物种分类和功能注释。（d）分箱。它涉及恢复未培养微生物的基因组草图，并重建系统发育和代谢通路。KEGG：京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes）；eggNOG：基因进化谱系：非监督直系同源群（Evolutionary genealogy of genes: non-supervised orthologous groups）；CAZy：碳水化合物活性酶数据库（Carbohydrate-active enzymes database）；CARD：抗性基因综合数据库（Comprehensive antibiotic resistance database）；VFDB：毒力因子数据库（Virulence factor database）。

7. 病毒组在人类疾病中的作用

近年来病毒组在人类疾病中的作用吸引了医学研究者的关注[121]。使用病毒组学的方法已发现了许多令人信服的研究成果[122]，其中一些技术已经用于临床[123]。在微生物组研究中，病毒组学与其他多组学方法整合后显示出广阔的应用前景。但是，病毒组学研究仍然面临一些挑战。例如，至少40％的病毒序列无法注释[124]。此外，病毒的测序结果容易受到背景噪音的影响[17]。最后，很难获得用于病毒组研究的商业化阳性对照，即病毒模拟群落[16]。

8. 总结和结论

本文讨论了用于微生物组研究的研究设计、样本收集、统计方法和生物信息学分析方法。在“研究设计”部分，我们强调了研究设计的重要性，特别是设计方案、样本量计算以及用于提高研究可靠性的多种措施。研究设计非常重要，因为不好的研究设计可能会产生无意义的数据。在“统计分析”部分，我们介绍了详细的多重比较P值校正方法。选择合适的统计方法对于准确解释微生物组数据很重要。最后，“生物信息学分析”部分介绍了用于分析微生物组数据分析的方法。本文图中使用的脚本可从 https://github.com/YongxinLiu/Qian2020CMJ 获得。

综上所述，对于微生物组研究而言，严谨的研究设计在获得有意义的结果方面具有举足轻重的作用，而适当的统计方法对于准确解释微生物组数据非常重要。循序渐进的分析流程为研究者掌握最新生物信息学分析方法提供了帮助。

参考文献

略，详见原文

Xu-Bo Qian, Tong Chen, Yi-Ping Xu, Lei Chen, Fu-Xiang Sun, Mei-Ping Lu & Yong-Xin Liu. (2020). A guide to human microbiome research: study design, sample collection, and bioinformatics analysis. Chinese Medical Journal Publish Ahead of Print, doi: https://doi.org/10.1097/cm9.0000000000000871

相关文章

• Protein Cell：扩增子和宏基因组数据分析实用指南
• 遗传：微生物组数据分析方法与应用
• 刘永鑫：想学菌群生物信息分析-21分钟带你入门

猜你喜欢

• 10000+: 菌群分析
  宝宝与猫狗 提DNA发Nature 实验分析谁对结果影响大  Cell微生物专刊 肠道指挥大脑
• 系列教程：微生物组入门 Biostar 微生物组  宏基因组
• 专业技能：生信宝典 学术图表 高分文章 不可或缺的人
• 一文读懂：宏基因组 寄生虫益处 进化树
• 必备技能：提问 搜索  Endnote
• 文献阅读 热心肠 SemanticScholar Geenmedical
• 扩增子分析：图表解读 分析流程 统计绘图
• 16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun
• 在线工具：16S预测培养基 生信绘图
• 科研经验：云笔记  云协作 公众号
• 编程模板: Shell  R Perl
• 生物科普: 肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA