1.Abstract

本文作者从实验设计，分子分析技术的选择，数据分析方法以及多种组学数据的解析等方面，对如何实现最优的微生物组学研究进行探讨。比如对近期快速发展的精确序列变异(exact sequence variants/ESV，详者注：目前更多使用ASV的名称)的方法替代传统基于聚类的OTU分析，整合宏基因组学和代谢组学的方法，组成型数据分析问题等方面的近期突破性的进展开展探讨。但在研究中还是应当关注实验设计和与研究可重复性相关的经典问题。

2.Introduction

和其他研究一样，在微生物组研究中，扎实的实验是至关重要的，实验方法、环境因素和分析都会影响最终结果。本领域当前研究获得了很多引人注目的成果，但仍然缺少数据收集和分析方法的标准。因此本领域目前最重要的挑战是，整合微生物组研究中独有的新方法，同时记得采用广泛应用于科学研究的标准方法。本综述旨在为微生物组实验设计和分析数据结果提供直接的指导标准，特别关注人类、模式生物以及环境微生物组。

3.Experimental design

设计可以获得有意义数据的实验是分析的第一步。典型的科学问题，例如疾病-对照(case–control)和,纵向干预(longitudinal interventions)研究等都可以放在微生物组的背景下研究。研究者可以分析在不同群落之间或时间序列下，微生物群落之间结构组成、遗传学或功能的潜在差异。需要考虑的重要问题是实验设计和样品收集**。严密的实验设计中隔离和询问感兴趣的可变因素等都是至关重要的。

微生物组实验设计中的注意事项

• 混杂因子对照：年龄、性别、饮食和生活方式。疾病组 vs 对照组
• 纵向取样：取样地点，季节变化-春、夏、秋、冬
• 实验技术引入的偏差：引物、空白对照、试剂差异或污染

首先必须确定实验范围，然后为感兴趣的问题选取适合的实验设计，例如，横向研究(cross-sectional studies)适用于发现不同人群(如健康和疾病)或生活在不同区域人群之间的微生物群落差异。基于分析计划和特定的科学问题来设计实验可以帮助我们确定样本量。（(推荐阅读：样品生物学重复数据选择 1必要性 2需要多少重复？).然而统计能力和效应量分析(effect size analysis)仍然是微生物组研究中的一大挑战。目前用于分析统计能力和效应量分析的方法大多基于置换多元回归方差分析(permutational multivariate analysis of
variance，PERMANOVA，原理详见：Adonis和ANOSIM方法组间整体差异评估原理)、狄利克雷-多项式分布(Dirichlet Multinomial)或者随机森林分析，详见的分类、回归实例。

3.Box 1 Good working practices

相似的微生物组研究常常产生相反的结果，如果没有详细的样品采集方法、实验设计、数据处理和分析过程的详细记录，就很难检查和解释出现差异的原因。
在收集样品时，采样的详细过程应当完整记录，并且应当考虑到更多的影响因素。另外，实验中要遵循基因组标准联盟提出的标记基因(marker gene)和宏基因组的基因组最小信息标准(minimum information standards, MIxS)：MIMARKS和MIMS(此外还有不可培养细菌、古菌的单扩增基因组Single Amplified Genome (MISAG)和宏基因组装基因组Metagenome-Assembled Genome (MIMAG)、末培养病毒基因Uncultivated Virus Genome (MIUViG)等标准可参考)。研究人员应该跟踪它们运行的所有命令和软件版本，并且将原始数据和样本元数据储存在公共数据库中。我们推荐使用Jupyter Notebooks，R Markdown等工具来实现这个目的，然后将其储存在GitHub等版本控制管理系统中。

4.Box 2 Considerations for different microbiomes

实验设计和方法的选择还是需要认真全面的考虑不同的样品类型。对于被非微生物DNA严重污染的样品，如植物、动物组织(通常宿主DNA占样本的90-99%，想要获得6 Gb微生物数据，理论上需要测序60 - 600 GB原始数据)等如果不排除掉宿主的DNA，鸟枪法宏基因组测序是不太可行的，如NBT封面：纳米孔宏基因组6小时识别下呼吸道病原体中采用皂苷去除99.99%宿主使病原体可以准确检测、Science：病原菌激活植物内生菌群的抑病功能中采用离心等方法富集根内生菌的方法。如果样品被死亡微生物等DNA遗迹严重污染(如土壤样品)，则需要在提取DNA之前使用物理方法来去除遗迹DNA(relic DNA)，例如使用单叠氮化丙锭。研究中收集、保存和储存的方法应该在所有样品中保持一致，以避免混淆因子引起的变异。

5.Defining controls and exclusion criteria

确定明确的选择和排除标准，可有效的限制混淆因子的影响。

6.Technical variation

从DNA提取到测序，不同实验方法之间的技术差异很大。在研究中所有样品必须使用相同的试剂盒，并且在纵向研究中应当收集多个基础样品用来评估时间点间在变异性。在采样、DNA提取、PCR和测序过程中，设计空白(阴性)对照对于监测污染至关重要。总之，使用不同方法产生的微生物组数据一致性依然是一个未能解决的难题。

7.Sequencing targets and methods

标记基因、宏基因组以及宏转录组测序研究微生物组会产生不同的结果。问题、假设、样品类型和分析目标都应该与所选的方法相匹配。
不同分析方法的优缺点。

标记基因分析(扩增子)
1.优点

• 样品制备和分析速度快、简单、成本较低
• 与基因组含量的相关性较高
• 适合于生物量较低、宿主DNA污染程度较高的样品
• 可用于与现有的大量公共数据集比较
  2.缺点
• 不能区分DNA来源中生物体是否有生命
• 受到扩增偏好性的影响较大
• 引物和可变区的选择对结果影响较大
• 要求对微生物群落有一定的先验知识
• 物种鉴定分辨率通常限于属水平
• 需要适当的阴性对照
• 获得的功能信息有限

全宏基因组测序
1.优点

• 可以直接获得微生物功能基因的相对丰度；基于已知物种可鉴定分辨率可达物种、甚至菌株水平
• 不需要微生物群落相关的先验知识(如捕获噬菌体、病毒、质粒以及微小真核生物等)
• 一般不会产生PCR偏好性
• 可以估算有参考基因组微生物的原位生长速率
• 可组装获得群体平均基因组(甚至可以获得其中一些微生物较完整的基因组)
• 可以挖掘新的基因家族
  2.缺点
• 成本相对较高，样品制备和分析较复杂
• 来自宿主和细胞器的DNA污染可能会掩盖微生物的特征
• 病毒和质粒通常无法自动化注释
• 与其他方法相比，通常需要较高的测序通量(几G - 几百G)
• 不能区分DNA来源于有生命或无生命的生物体
• 由于受组装影响，平均群体微生物基因组往往不准确

全宏转录组测序
1.优点

• 当与标记基因分析结合使用时，可以估算群落中哪些微生物正在进行积极的转录过程
• 只能鉴定活动生物，排除休眠、死亡微生物及胞外DNA
• 能够捕捉个体内部的动态变化
• 直接评估微生物的活性，包括对干扰或者暴露等情况的响应
  全宏转录组测序
  2.缺点
• 费用最高，样品制备和分析过程最复杂
• 必须排除宿主的mRNA、和rRNA污染
• 样品的收集和存储要十分小心
• 数据结果对有高转录率的生物体有偏向性
• 需要与DNA测序结果结合，才能获得细菌丰度变化和转录率

16S测序后，我们推荐使用Deblur获得单碱基变异的参考序列(sOTUs),尽管DADA2与Deblur结果类似，但Deblur支持并行处理速度更快且更稳定(在不同样品中获得相同sOTUs)。宏基因组和宏转录组首先要去除宿主DNA或rRNA，和宿主RNA。过滤后的序列可以采用Kraken、MEGAN或HUMAnN等有参方法(read-based)，或De novo组装的方法metaSPAdes和MEGAHIT分析。

8.Marker gene analysis

这个区域通常包含高度可变区，可用于区分研究对象的组成，并且两侧包含可以用作PCR引物结合位点的高度保守区。但是由于这些引物扩增区域的DNA序列不同，可能对DNA序列的亲和力不同产生偏好性，从而影响PCR扩增结果。。标记基因测序中的偏好性来源可能是由于不同的可变区选择、扩增子片段大小和PCR循环次数等。关于扩增子引物选择，可进一步阅读：16S结构 16S单V4区是最佳选择?引物评估等文章。

9.Whole metagenome analysis

宏基因组分析就是对样本内所有微生物基因组进行测序的方法。宏基因组测序在基因水平上获得整个群落功能的能力远超标记基因可分析的范围。关于宏基因组学的全面综述，推荐阅读2017年自然生物技术的综述：宏基因组从取样到分析，详见《Nature综述：2万字带你系统入门鸟枪法宏基因组实验和分析》

10.Metatranscriptome analysis

宏转录组分析是通过使用RNA测序来分析微生物组的转录过程，从而提供关于基因表达和微生物组功能活性等信息。值得注意的是，采用叠氮溴化丙锭(propidium monoazide, PMA)去除遗迹DNA的方法也是获得活性微生物组的可选方法之一。如果你想全面了解宏转录组学分析，请阅读《使用宏转录组进行微生物组研究》的文章(Bashiardes, S., et al. 2016. Bioinform. Biol. Insights.)。

11. Marker gene analyses

分析标记基因数据的第一步是去除序列错误：尽管序列错误率很低，在Illumina测序中，每个核苷酸的错误率仅为 ~ 0.1%，但是很大部分明显的序列多样性来源于测序错误（如1M碱基可能拥有1000个测序错误，造成增长成百上千的多样性；大规模的实验测序量可达 Billion）。OTUs聚类，即将相似的序列(通常具有97％相似性阈值)合并归为单个的特征，然后将序列的变体(包括通过序列错误引入的序列变体)合并成可用于随后分析的单个OTU。但是这种方法会在一定程度上，遗漏一些细微但真实的生物序列变异，例如存在SNP的序列本该为多个独立OTUs。寡聚分型基于16S rRNA基因测序中位置的特异性信息，来鉴定单碱基变异(SNP)从而加以区分非常相似但不同的分类群。诸如Deblur和DADA2等算法，使用测序错误校正的模型来转换测序数据为精确序列(标记基因序列)，也称为亚-OTUs(sOTUs)。这些方法得到的结果是一个DNA序列表，是每个样品中的不同序列数，而不是OTU群组。因此我们推荐，当需要与常见的全长数据参考数据库比对的时候，这些方法替代现有基于OTU的方法，除非需要组合使用不同技术(即Illumina测序和454焦磷酸测序)产生的测序数据或者是引物不同。
一个关键的分析步骤是为微生物序列进行物种分类注释。原则上，与三大参考数据库(三个最具特色且经常使用的是Greengenes，RDP和Silva)精确匹配应当提供更好特异性的分类学分配，但当存在大量未知的分类群时这种方法的敏感性较差。此外，由片段较短的标记基因构建的系统发育树通常结果较差，将标记基因序列插入到基于全长序列的参考序列系统发育树中是一种更好的做法。另外，应当对未分类的微生物进行核糖核酸序列分析是否为细胞器的序列，如叶绿体、线粒体**(宿主非特异扩增序列)。
功能预测分析是一种将标记基因和可用的微生物基因组相联系的技术，用来预测宏基因组，从而推断其生物功能。这种分析通常需要基于参考数据库生成OTU表，然后基于演绎模型(如PICRUSt,其中PICRUSt2：OTU/ASV等16S序列随意预测宏基因组，参考数据库增大10倍)为这些基因含量预测提供置信区间，即在距离参考基因组较远的树置信度低，而在许多参考基因组可用的区域则置信度高。

12. Metagenome and metatranscriptome analyses

将未组装的DNA或mRNA序列与参考数据库比对，可以获得物种和功能基因注释.Cell：20种宏基因组学物种分类工具大比拼,物种分类标记基因方法采用广泛关注的单拷贝基因，MetaPhlAn2一条命令获得宏基因组物种组成。
想要组装获得部分或完善的单菌基因组，长序列通常采用MaxBin2，或CONCOCT进行分箱(binning，或分类)为假定的单菌基因组，分箱主要原理是基于丰度和核酸组成等信息。
此外，有许多工具可以进行更为复杂的标准化方法，如edgeR和DESeq2(采用基于负二项分析的标准化方法，在高通量测序数据领域应用极广泛，edgeR使用实战详见3热图：差异菌、OTU及功能)。

13. Higher-level analyses

微生物组数据经过处理，可以获得特征(features，如物种不同分类级或基因)与样本的丰度矩阵。但这一结果是存在迷惑性(deceptively)的，因为微生物组数据通常是高维数据，包括几千个不同物种，表格存在许多零值的稀疏性特点；因此需要注意的统计处理方法，以挖掘有意义的结果。Alpha和Beta多样性常用于评估微生物组的整体变异。Alpha多样性可以量化样品内的特征多样性，也可以进行样品组间比较。Alpha多样性物种测量的方法有三类：丰富度(richness)的测量常用观测的物种数(Observed OTU / Richness)和Chao1丰度估计(估计真实物种多样性)，进化距离测量采用信任系统发育多样性(Faith’s phylogenetic diversity)，这两类方法受样本测序深度影响很大；此外还有一类即考虑丰富度，又考虑均匀度的Shannon指数，对测序量不敏感。详见箱线图：Alpha多样性解读。Beta多样性比较每对样品间的差异，产生所有样品对间的距离矩阵。度量标准的选择对结果影响较大，需牢记我们在挖掘生物学数据并关注其背景的意义。Bray-Curtis、Canberra、 有权重的UniFrac等定量度量采用特征的丰度信息进行计算，binary-Jaccard、无权重的UniFrac定性方法仅考虑特征的有无。进化方法的Unifrac分析需要进化树文件，可提供更生物学的解析，但缺少树文件时无法使用。(Alpha多样性组间常用ANOVA)，无参数的置换(permutation)检验方法PERMANOVA、ANOSIM用于估计的不同组间beta多样性的显著性，其中PERMANOVA应用于组间变异较大的数据集更好用。计算Alpha和beta多样性，需要研究者掌握抽样技术(即每个样本中抽取相同数量的序列)，不同的抽样数量级可影响结果。目前计算Unifrac最好的方式是稀疏/稀释(rarefracation)，但一些特殊情况下的成对差异丰度比较需要完整的样本数据集。Beta多样性数据可视化采用排序的技术，常用如主坐标轴分析(PCoA)或主成分分析(PCA)。点我读懂PCA和PCoA、PCoA绘制时间序列R代码实例 。这类方法将复杂的距离矩阵，转换为可观察的2或3维空间，代表样品间距离。样品可以按分组信息着色，方便观察组间差异，属于无监督的方法。EMPeror框架提供可交互式的显示PCoA图。微生物组数据具有高维、松散、组成型等特点，鉴定解析微生物群落差异的分类群具有挑战性。例如，己知某个病人的药物只影响一个微生物属，对其它菌无任何影响。尽量其它微生物不受药物影响，但它们相对丰度减少，是由单个微生物属过度生长引起的。这种情况影响许多经典方法的结果，如参数统计检验(如student’ t-test和ANOVA，还有wilcoxon秩和检验)、计算相关性(如Spearman排序相关系数)通常导致完全不可接受的超90%假阳性率。一种方法是在统计检验上强制进行强生物假设：如Lovell’s比例度量方法仅检测正相关。其它一些工具为微生物组数据专门做了优化，假定小部分物种是相关的，大多数的相关系数为0，如SparCC和SPEIC-EASI。机器学习是在微生物组领域非常有效的方法，可基于当前状态区分样品(分类，由己知的分类与结果学习，预测末知分类，如健康和疾病、亚种分类【NBT封面：预测水稻亚种】)。ourceTracker可以估计末知群体微生物来源和组成，最有用的是可根据环境样品来分类微生物的来源。

14.Integrating other omics data

了解微生物群落的组成并不是研究的终点，我们更想知道群体的功能。扩增子测序，宏基因组，宏转录组，宏蛋白组，宏代谢组和其它技术的多组学数据整合，可用于特定微生物群体功能和组成的深入研究。
Nature Protocols：整合宏基因组、代谢组和表型分析的的计算框架综上所述，整合多组学数据可以更全面的理解微生物组——从DNA鉴定到蛋白和代谢物的功能，使用研究结果可有科学意义。

14.Metabolomics and the microbiome

微生物产生代谢物可影响宿主和微生物群体动态变化，并与宿主的疾病和健康状态有关。代谢物有益处(如短链脂肪酸)或毒性(基因毒率大肠杆菌素)影响宿主

15.Conclusions

本综述讨论了微生物组研究各阶段工作的指南，从实验设计、收集储存样品、测序数据的图形结果中挖掘规律等，均对结果与生物学解释有影响。

本文主要参考文献：Best practices for analysing microbiomes
以及参考公众号刘永鑫Adam