微生物组学研究手段概览——扩增子测序
林二狗 宇宙实验媛 2018-10-23
微生物组学不依赖于微生物的分离培养，克服了传统的纯培养方法的技术限制，为研究和开发利用不可培养的微生物（占微生物种类99%以上）提供了一种新的途径和良好的策略。另外，通过微生物组学可以得到环境中丰度较低的，甚至是痕量微生物的信息。同时结合宏观生态的研究理念，可以解释环境中微生物菌群的多样性、功能活性等宏观特征，相比于对于单一种属微生物的研究，能更准确地反应出微生物生存的真实状态。
本次系列专题旨在概述微生物组学的主要研究手段，包括扩增子测序、宏基因组、宏转录组、宏蛋白质组、宏代谢组、功能组和培养组，每篇文章除了基础的介绍内容，还结合了刚刚发表的相应文章进行实例解析，敬请关注。
微生物扩增子测序：样本中微生物种类及构成、不同样本间微生物组成差异
微生物扩增子测序主要通过直接扩增环境总 DNA 的特定区域，解释某一特定样本环境微生物分布、丰度变化和群落组成情况。随着高通量测序技术的不断发展，微生物扩增子测序技术已成为环境微生物群落比较和差异分析的主流研究手段之一。
根据不同的研究需求，一般选择扩增16S高变区来区分环境样品的细菌和古菌群落，ITS区域扩增用以评估真菌群落，原生动物等为主的真核微生物群落研究可通过18S高变区测序来实现。此外，也可通过特定的功能基因测序等来揭示特定功能相关的环境微生物群落分布。由于在下主攻细菌和古细菌，所以下面多以此为例。
首先我们来看看为什么选择扩增16S片段来表征菌群。细菌核糖体RNA（rRNA）有三种类型：5S rRNA（120bp）、16S rRNA（约1540bp）和23S rRNA（约2900bp）。其中，5S rRNA基因序列较短，包含的遗传信息较少，不适于细菌种属的分析鉴定；23S rRNA基因的序列太长，且其碱基的突变率较高，不适于鉴定亲缘关系较远的细菌种类；而16S rRNA普遍存在于原核细胞中，种属内具有高度的保守性，且含量较高、拷贝数较多，遗传信息量适中，既能体现不同菌属之间的差异，又能利用测序技术较容易地得到其序列，非常适合作为细菌多样性分析的标准。
图1. 16S rRNA基因区域：包含9个高变区（V1-V9）和10个保守区。通过对某一段高变区序列或几个高变区的组合（比如，V4区或V3+V4区）进行PCR扩增后进行测序。
图2. 16S rRNA基因不同区段和组合对菌群分类的解释效果（doi:10.1038/nrmicro3330）
这里介绍一个概念——OTU。OTU（operational taxonomic units），即操作分类单元，通过一定的距离度量方法计算两两不同序列之间的距离度量或相似性，继而设置特定的分类阈值，获得同一阈值下的距离矩阵，进行聚类操作，形成不同的分类单元。（隐约听到有人让我说人话……）对于样本的16S序列进行测序，往往会得到成千上万条reads，如果直接对每条序列进行物种注释的话，工作量大、耗时长。而且16S序列的扩增、测序等过程中出现的随机错误会降低结果的准确性。所以，在16S分析中引入OTU，首先对相似性（阈值一般为97%）序列进行聚类，分成数量较少的分类单元，基于分类单元进行物种注释。对应关系大概是酱婶儿的：
但是需要说明的是：（1）因为注释基于的是扩增子序列数据库（如RDP、Greengenes、SILVA等），所以有很多未知的微生物种类不能被鉴定出来，尤其是土壤等环境样本中，这些在文章中都被归类为“unclassified”；（2）目前微生物扩增子测序普遍使用二代测序仪（三代测序测16S全长了解一下），因为受到读长限制，注释信息往往不能精确到种水平。
表1. 测序公司提供的基于扩增子测序数据的标准分析内容
以Ferrocino和Ponzo等人2018年8月发表于Scientific Reports的文章为例。研究人员招募了41例妊娠糖尿病患者，分别于孕中期（24-28周胎龄）和孕晚期（38周胎龄）收集粪便，提取肠道微生物组，进行16S扩增子测序，解析这些患者在怀孕过程中的肠道菌群变化。
（1）α多样性分析：
α多样性主要关注局域均匀生境下的物种数目，因此也被称为生境内的多样性（within-habitat diversity）
shannon：菌群多样性指数
H=－∑（Pi）（㏑Pi），Pi=样品中属于第i种的个体的比例，如样品总个体数为N，第i种个体数为ni，则Pi=ni/N。
各种之间，个体分配越均匀，H值就越大。换言之，如果每一个体都属于不同的种，多样性指数就最大；如果每一个体都属于同一种，则其多样性指数就最小。
observed_species：菌种丰富度指标，观测到的OTU数
chao1：菌种丰富度指数，估计群落中的OTU数目
chao1是度量物种丰富度的指标，它和丰度、均匀度无关，但是它对稀有的物种很敏感。
Schao1=Sobs+n1(n1-1)/2(n2+1)，其中Schao1为估计的OTU数，Sobs为观测到的OTU数，n1为只有一条序列的OTU数目，n2为只有两条序列的OTU数目。
图3. α多样性结果展示：绿色（Enrolment）代表24-28周胎龄样本，蓝色（Study end）代表38周胎龄样本（下同）。由图中，可以看出相比于Enrolment组，Study end组的α多样性有所增加。
（2）β多样性分析：
β多样性指沿环境梯度不同生境间群落的物种组成的相异性或物种沿环境梯度的更替速率，也被称为生境间的多样性（between-habitat diversity）。
β多样性分析常用的方式为PCA（Principal Components Analysis）和PcoA（Principal Co-ordinates Analysis），二者都属于排序分析（Ordination analysis）。排序(ordination)的过程就是在一个可视化的低维空间或平面重新排列这些样本,使得样本之间的距离最大程度地反映出平面散点图内样本之间的关系信息（反映样本间菌群结构的相似性和差异性）。每一个点代表一个样本，相同颜色的点来自同一个分组。通过比较样本点的距离来衡量样本/组间的差异程度，样本/组间的距离越近表示两样本/组的组成相似性越高，差异越小。
PCA即主成分分析，也称主分量分析或主成分回归分析法。首先利用线性变换，将数据变换到一个新的坐标系统中；然后再利用降维的思想，使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上，第二大方差在第二个坐标(第二主成分)上。这种降维的思想首先减少数据集的维数，同时还保持数据集的对方差贡献最大的特征，最终使数据直观呈现在二维坐标系。
PCoA即主坐标分析，它与PCA类似，通过一系列的特征值和特征向量进行排序后，选择主要排在前几位的特征值，找到距离矩阵中最主要的坐标，结果是数据矩阵的一个旋转，它没有改变样本点之间的相互位置关系，只是改变了坐标系统。
两者的在于PCA分析是基于原始的物种组成矩阵所做的排序分析，而PCoA分析则是基于由物种组成计算得到的距离矩阵得出的。当样本数量比较多而物种的数量比较少的时候，推荐选择PCA；反之，样本数量少而物种数量多的时候，选择PCoA更好。当然，很多时候我们会PCA和PCoA都做，根据结果再来选择。
图4. 主成分分析结果展示。如图所示，x轴标签PCo 1 (22.9%)代表能最大区分所有样品的第一主坐标轴，可以解释样品中所有差异的22.9%；y轴标签PCo 2 (23.5%)代表能最大区分所有样品的第二主坐标轴，可以解释样品中所有差异的23.5%。这两轴形成的第一个平面展示了样品间近50%以上的差异，不同颜色表示两组间可以很好的区分开，即妊娠糖尿病患者的肠道微生物组成与其所处的不同孕期具有显著相关性。
（3）预测宏基因组的功能基因组成：
PICRUSt（Phylogenetic Investigationof Communities by Reconstruction of Unobserved States）是基于16S rRNA基因序列预测微生物群落功能的工具。其原理是首先对基因内容进行预测（Gene Content Prediction；对Greengenes数据库中的“closed reference”序列划分OTU后构建进化树，因为基因内容在进化距离相近的物种间相对保守，所以可以通过Ancestral state reconstruction算法并结合IMG/M数据库预测出树中未进行全基因组测序的OTU的基因组信息）。然后进行宏基因组预测（Metagenome Prediction；a. 因为细菌中可能含有1-15个16SrRNA基因拷贝，按照预测其他基因方式预测16SrRNA基因拷贝数；b. 将每个OTU对应序列数除以其16SrRNA基因拷贝数来进行标准化；c. 将标准化的数据乘以其各自对应的基因组中基因数从而实现宏基因组预测的目的）。最后，获得的预测结果可以通过KEGG Orthology或COGs等对基因家族进行分类。这样就可以在不进行宏基因组测序（详见后文）的情况下，对微生物组的功能进行一定程度地表征。
图5. PICRUSt工作流程（doi: 10.1038/nbt.2676，下同）
图6. 使用PICRUSt预测不同环境的微生物组的准确性：人体>土壤>其他哺乳动物>高盐
图7. 随测序深度增加PICRUSt预测结果与宏基因组测序结果的准确度比较。图中所示使用PICRUSt，即使低16S测序量（105条序列）也可获得准确度较高的功能预测结果。
为什么说是“一定程度”呢？因为细菌间普遍存在水平转移，即使16S序列相同，不同菌株所含有的功能基因也可能有所差异。所以，PICRUSt预测的结果仅供描述。同时，如果在Greengenes数据库中找不到所测序列的同源物种的参考序列，则该序列的功能基因将无法被预测，而且PICRUSt只能对已知微生物基因的已知功能进行预测，故其并不能完全代替宏基因组研究。但是，宏基因组测序费用相对昂贵，在经费有限的情况下，这种预测方式不失为一种替代的选择。
图8. 使用PICRUSt预测的Study end组菌群代谢通路与OTU相关性分析。由图，糖酵解/糖异生、果糖/甘露糖代谢、半乳糖代谢、淀粉/蔗糖代谢等通路相对富集，脂肪酸代谢等通路则降低。同时，脂多糖（LPS）的生物合成通路与Sutterella和Bacteroides等菌属显示出正相关。
微生物组学研究手段概览2——宏基因组和宏转录组
宏基因组
宏基因组测序是将环境总DNA提取出来，随机打断成300/500bp的小片段，然后在片段两端加入通用引物进行PCR扩增测序，然后对测序数据进行质控，再将高质量序列拼接，根据数据库参考信息，对基因序列进行预测和功能注释，最终获得重要的宏基因组信息，如序列组成（GC含量、基因组大小等）、物种组成、功能组成和群落特征等。
相比于16S扩增子测序，宏基因组测序能够使物种鉴定深度达到“种”，而前者往往只能达到“属”的级别。另外，基于16S扩增子测序的基因预测结果（预测方法请参见往期内容“微生物组学研究手段概览——扩增子测序”）是依据数据库中的参考序列得到的，宏基因组测序则提供了菌群实际的基因信息。所以，如果不考虑经费的限制，宏基因组测序是能够更准确地研究微生物组及其功能的方法。
图1. 宏基因组学研究的生物信息工作内容（doi: 10.1371/journal.pcbi.1002808）。
图2. 宏基因组数据处理的基本流程和需用软件（有具体操作需求的同学可以通过下面网址下载英文相关教程
https://github.com/TGAC/361Division/tree/master/Metagenomics%202015)
表1. 测序公司提供的基于宏基因组测序数据的标准分析内容
目前很多研究都聚焦在人体肠道微生物组的宏基因组学研究，期望从中挖掘出多种疾病的因果关系。已经采集到的样本来自不同国家、遗传背景、生活习惯、身体状况的人群。同时，环境微生物的宏基因组学研究也在如火如荼地开展中。相比于宿主样本，环境样本（如水体、土壤等）中存在大量未知或不可培养的微生物，这些微生物有很多是尚未被鉴定或者深入研究的。宏基因组测序直接研究环境中的总DNA，为开发新的生物活性物质、发现新的基因和物种、研究特定环境中微生物群落结构与功能的关系（功能网络与互作）、微生物对环境变化的响应与反馈（比如，农耕土壤质量的恶化与恢复）、微生物群落的演替与进化、微生物区域分布与生物地理学等开辟了一条新的途径，为解释和解决一些重大农业和环境问题提供重要依据。
图3. NCBI Sequence Read Archive鸟枪法宏基因组测序数据的增长（doi: 10.1016/j.cell.2016.08.007）
Nature文章实例
以Crits-Christoph和Diamond等人2018年6月发表于Nature的文章为例。研究人员用宏基因组测序的方式，从草地土壤样本中获得了376株细菌基因组，的基因组序列，发现了完全新颖的生物合成基因簇，并描述了它们的基因组学、系统发育学和生态学状况。他们共鉴定出了1599个生物合成基因簇（biosynthetic gene clusters），推测这些基因簇可能合成非核糖体多肽、聚酮化合物等，很多非核糖体多肽合成酶（NPRS）和聚酮合成酶（PKS）都能够合成抗生素、抗真菌素、嗜铁素或免疫抑制剂，所以该研究中着重关注这两类基因簇（图4），以期发现新的抗生素和药用化合物。
图4. 样本中提取到的基因组和其生物合成基因簇的情况差异展示。a，几个菌门的平均相对丰度。b，生物合成基因簇在几个菌门之间的分布，不同颜色代表由antiSMASH（用于生物合成基因簇的分析；http://antismash.secondarymetabolites.org）推测的产物类型c，研究鉴定到了240个NRPS、PKS（根据酶结构域的组成差异分为Ⅰ型、Ⅱ型和Ⅲ型）和NRPS-PKS混合型基因簇，以及86个可能不完整的基因簇。尽管这些基因含量水平各异，但是由于这些酶具有保守的结构域，因此他们的生物合成途径是可识辨的。d，生物合成基因簇的网络图，连线代表两端具有共有基因，共有基因所占比例越大连线越粗、颜色越深。由该图可知，Verrucomicrobia、Acidobacteria和Rokubacteria中存在不同且稀疏的NRPS和PKS系统，而多数稀有型NRPS基因簇之间的距离较为疏远。在Rokubacteria和Acidobacteria的一支中保守的Ⅲ型PKS基因簇形成了一个密集的网络集群，高度的保守性可能意味着一类新的代谢物的广泛分布。
图5. Acidobacteria中的生物合成NRPS/PKS基因簇。Acidobacteria中发现了两个几乎含有完整NRP和PKS基因簇的基因组，它们分别被命名为Candidatus Eelbacter (Eelbacter_gp4_AA13)和Candidatus Angelobacter (Angelobacter_gp1_AA117)。通过对Acidobacteria基因组中核糖体蛋白序列进行系统发育学分析，发现两个基因组均进化出独特的生物合成操纵子(5a)。Candidatus Angelobacter基因组中含有多种抗生素合成蛋白、一种细菌素合成簇、Ⅵ型和Ⅱ型分泌系统的多基因操纵子元件和一些含有RHS重复结构的大蛋白；Candidatus Eelbacter基因组包含6个长度超过45kb的复合Ⅰ型NRPS-PKS混合型基因簇(5b)。
宏基因组测序有点儿贵，老板看了价格有点儿心碎。所以为了节约成本，很多同学就只能走上自己建库的道路了。然而建库质量的优劣直接影响测序数据的信息完整性和正确性，所以选择合适的建库方法和试剂就显得尤为重要了。
这里，小编推荐VAHTS™ Universal DNA Library Prep Kit for Illumina® V3快速通用型文库构建试剂盒。该试剂盒可用于100 pg - 4 μg Input DNA建库，并且通过对经典DNA建库流程的末端修复模块、连接模块和文库扩增模块等环节的整体改进，使文库转化率和扩增文库产出得到大幅提升，广泛适用于各种样本的PCR或PCR-Free文库构建。试剂盒中提供的所有试剂都经过了严格的质量控制和功能验证，最大程度上保证了文库构建的稳定性和重复性。
以上试剂盒可免费申请试用，数量有限先到先得哦。如需申请，请撩小编。
宏转录组
宏转录组测序是指从整体水平上研究某一特定环境、特定时期群体生命全部基因组转录情况以及转录调控规律的研究手段。它以生态环境中的全部RNA为研究对象，提取环境微生物群落中的全部转录本，进行高通量测序和生物信息学分析，与宏基因组研究相辅相成，能够很好地揭示复杂微生物群落的变化，有效地扩展微生物资源的利用空间。
图6. 宏转录组分析简化流程（doi: 10.1007/s00253-018-8976-7）
表2. 测序公司提供的基于宏转录组测序数据的标准分析内容
Nature文章实例
还是同一篇文章，结合宏转录组测序和Kallisto量化分析发现，在133个差异表达的NRPS和PKS基因簇中，共检测到198个NRPS/PKS基因的表达。在4个微生物门类中均检测到了NRPS和PKS的表达，并且酸杆菌门有84个活性基因簇产生了表达。本文还检测了Candidatus Eelbacter的10个生物合成基因簇（含11个NRPS/PKS结构域基因）和Candidatus Angelobacter的25个生物合成基因簇（含25个NRPS/PKS结构域基因），发现7个基因组中有10个NRPS/PKS基因簇在24h修正试验过程中，表现出时间依赖性。
Candidatus Angelobacter中有5个生物合成基因簇与多种感知、响应环境的基因具有共表达模型，如ToB(嗜铁素摄取受体)、MacB(大环内酯转运物)、pbp(青霉素结合蛋白)、16s rRNA MT(16s rRNA甲基转移酶)和gvp(气泡蛋白)等(图7c)。Angelobacter多个基因簇的生物合成基因表达是同步的，表明这些基因能够协同响应生态竞争。此外，Acidobactera_gp22_AA4和Gemmatimonadetes_AG49也具有多个重要的共表达基因，如rsb X/R/S(应激响应调节操纵子)、vgb(维及霉素同源酶B)和ToB等(图7c)。
图7. 生物合成基因的宏转录组。在添加了底物后12-24h，Candidatus Angelobacter基因组中几个基因簇的基因表达水平显著升高（7a），并且一些生物合成基因的表达与核心核糖体基因的表达不一致（7b）。表明Candidatus Angelobacter能够对水和基质的添加产生响应，并且在核心代谢基因的表达量上升数小时后单独调控次级代谢产物合成相关基因的表达。从基因共表达模型来看，7个基因组中有4个发生了显著的次级代谢基因共表达富集（p＜0.05）（7c）。
简单来说，宏基因组能够说明菌群“能做什么”，而宏转录组则是“要做什么”。尤其是在不同取样时间点之间和实验条件处理前后，宏转录组往往会发生明显变化，据此可以建立菌群功能与影响因子之间的关联性。后面我们还会介绍表征菌群“做了什么”和“做出了什么”的研究方法，敬请关注。
转载出自林二狗宇宙实验媛 ，如有侵权请及时联系。
