微生物群落多样性测序与功能分析

微生物群落多样性测序与功能分析微生物群落测序是指对微生物群体进行高通量测序，通过分析测序序列的构成分析特定环境中微生物群体的构成情况或基因的组成以及功能。借助不同环境下微生物群落的构成差异分析我们可以分析微生物与环境因素或宿主之间的关系，寻找标志性菌群或特定功能的基因。对微生物群落进行测序包括两类，一类是通过16srDNA，18srDNA，ITS区域进行扩增测序分析微生物的群体构成和多样性；还有一类是宏基因组测序，是不经过分离培养微生物，而对所有微生物DNA进行测序，从而分析微生物群落构成，基因构成，挖掘有应用价值的基因资源。以16srDNA扩增进行测序分析主要用于微生物群落多样性和构成的分析，目前的生物信息学分析也可以基于16srDNA的测序对微生物群落的基因构成和代谢途径进行预测分析，大大拓展了我们对于环境微生物的微生态认知。目前我们根据16s的测序数据可以将微生物群落分类到种（species）（一般只能对部分菌进行种的鉴定），甚至对亚种级别进行分析，几个概念：16SrDNA（或16SrRNA）：16SrRNA基因是编码原核生物核糖体小亚基的基因，长度约为1542bp，其分子大小适中，突变率小，是细菌系统分类学研究中最常用和最有用的标志。16SrRNA基因序列包括9个可变区和10个保守区，保守区序列反映了物种间的亲缘关系，而可变区序列则能体现物种间的差异。16SrRNA基因测序以细菌16SrRNA基因测序为主,核心是研究样品中的物种分类、物种丰度以及系统进化。OTU：operationaltaxonomicunits(OTUs)在微生物的免培养分析中经常用到，通过提取样品的总基因组DNA，利用16SrRNA或ITS的通用引物进行PCR扩增，通过测序以后就可以分析样品中的微生物多样性，那怎么区分这些不同的序列呢，这个时候就需要引入operationaltaxonomicunits，一般情况下，如果序列之间，比如不同的16SrRNA序列的相似性高于97%就可以把它定义为一个OTU，每个OTU对应于一个不同的16SrRNA序列，也就是每个OTU对应于一个不同的细菌（微生物）种。通过OTU分析，就可以知道样品中的微生物多样性和不同微生物的丰度。测序区段：由于16srDNA较长（1.5kb），我们只能对其中经常变化的区域也就是可变区进行测序。16srDNA包含有9个可变区，分别是v1-v9。一般我们对v3-v4双可变区域进行扩增和测序，也有对v1-v3区进行扩增测序。工具/原料16srDNA测序首先需要提取环境样品的DNA，这些DNA可以来自土壤、粪便、空气或水体等任何来源。提取DNA后需要经过质检和纯化，一般16srDNA测序扩增对DNA的总量要求并不高，总量大于100ng，浓度大于10ng/ul一般都可以满足要求。如果是来自和寄主共生的环境如昆虫的肠道微生物，提取时可能包括了寄主本身的大量DNA，对DNA的总量要求会提高。微生物菌群多样性测序受DNA提取和扩增影响很大，不同的扩增区段和扩增引物甚至PCR循环数的差异都会对结果有所影响。因而建议同一项目不同样品的都采用相同的条件和测序方法，这样相互之间才存在可比性。完成PCR之后的产物一般可以直接上测序仪测序，在上机测序前我们需要对所有样本进行定量和均一化，通常要进行荧光定量PCR。完成定量的样品混合后就可以上机测序。16srDNA测序目前可以采用多种不同的测序仪进行测序，包括罗氏的454，Illumina的MiSeq，Life的PGM或Pacbio的RSII三代测序仪。不同的仪器各有优缺点，目前最主流的是Illumina公司的MiSeq，因为其在通量、长度和价格三者之间最为平衡。MiSeq测序仪可以产生2x300bp的测序读长，一次可以产生15Gb的测序数据远远大于其他测序仪的测序通量。方法/步骤1.116srDNA分析基本流程：2.2原始数据处理：原始测序数据需要去除接头序列，并将双端测序序列进行拼接成单条序列。根据测序barcode序列区分不同的样本序列。过滤低质量序列和无法比对到16srDNA数据库的序列。3.3OTU分类和统计：OTU(operationaltaxonomicunits)是在系统发生学研究或群体遗传学研究中，为了便于进行分析，人为给某一个分类单元（品系，种，属，分组等）设置的同一标志。通常按照97%的相似性阈值将序列划分为不同的OTU，每一个OTU通常被视为一个微生物物种。相似性小于97%就可以认为属于不同的种，相似性小于93%-95%，可以认为属于不同的属。样品中的微生物多样性和不同微生物的丰度都是基于对OTU的分析。使用QIIME（version1.8.0）工具包进行统计注释。使用QIIME（version1.9.0,）的ucluster方法根据97%的序列相似度将所有序列进行同源比对并聚类成operationaltaxonomicunits(OTUs)。然后与数据库GreenGenes（versiongg_13_8,）进行比对，比对方法uclust，identity0.9。然后对每个OTUs进行reads数目统计。下面的2个表，其中一个表是对每个样本的测序数量和OTU数目进行统计，并且在表栺中列出了测序覆盖的完整度（显示前10个样本）。另一个表是对每个样本在分类字水平上的数量进行统计，并且在表栺中列出了在每个分类字水平上的物种数目（显示前10个样本）。可以看到绝大部分的OTU都分类到了属（Genus），也有很多分类到了种（Species）。但是仍然有很多无法完全分类到种一级，这是由于环境微生物本身存在非常丰富的多样性，还有大量的菌仍然没有被测序和发现。测序数目统计表主要是对每个样本的测序数量和OTU数目进行统计，并且在表格中列出了测序覆盖的完整度（显示前10个样本，如果样本超过10个，请查看结果中otu_stat.txt文件）其中SampleName表示样本名称；SampleSize表示样本序列总数；OTUsNumber表示注释上的OTU数目；OTUsSeq表示注释上OTU的样本序列总数。Coverage是指各样品文库的覆盖率，其数值越高，则样本中序列没有被测出的概率越低。该指数实际反映了本次测序结果是否代表样本的真实情况。计算公式为：C=1-n1/N其中n1=只含有一条序列的OTU的数目；N=抽样中出现的总的序列数目。分类水平统计表主要是对每个样本在分类学水平上的数量进行统计，并且在表格中列出了在每个分类学水平上的物种数目（只显示前10个样本，如果样本超过10个，请查看结果中taxon_all.txt文件）其中SampleName表示样本名称；Phylum表示分类到门的OTU数量；Class表示分类到纲的OTU数量；Order表示分类到目的OTU数量；Family表示分类到科的OTU数量；Genus表示分类到属的OTU数量；Species表示分类到种的OTU数量。4.4我们还可以对这些种属的构成进行柱状图显示：横坐标中每一个条形图代表一个样本，纵坐标代表该分类层级的序列数目或比例。同一种颜色代表相同的分类级别。图中的每根柱子中的颜色表示该样本在不同级别（门、纲、目等）的序列数目，序列数目只计算级别最低的分类，例如在属中计算过了，则在科中则不重复计算。Q:为什么要选择V3-V4区的测序长度？为什么有些文献是V6区，有什么区别？A:16SrRNA总长约1540bp，包含9个可变区。由于高通量测序的测序长度的限制，不可能将16SrRNA的9个可变区全部测序，所以在PCR扩增时往往只能选择1-3个可变区作为扩增片段。Kozich等评估了Miseq测序仪分析的不同16SrRNA可变区的准确性发现，测定V4区效果最佳。根据我们的测序长度，v3-v4区是最佳选择。5.5我们还需要对样本之间或分组之间的OTU进行比较获得韦恩图：注意，韦恩图目前一般最多只能显示5个样本或分组，过多的样本无法无法进行韦恩图绘制6.6样品构成丰度：稀释曲线微生物多样性分析中需要验证测序数据量是否足以反映样品中的物种多样性，稀释曲线（丰富度曲线）可以用来检验这一指标。稀释曲线是用来评价测序量是否足以覆盖所有类群，并间接反映样品中物种的丰富程度。稀释曲线是利用已测得16SrDNA序列中已知的各种OTU的相对比例，来计算抽取n个（n小于测得reads序列总数）reads时出现OTU数量的期望值，然后根据一组n值（一般为一组小于总序列数的等差数列）与其相对应的OTU数量的期望值做出曲线来。当曲线趋于平缓或者达到平台期时也就可以认为测序深度已经基本覆盖到样品中所有的物种；反之，则表示样品中物种多样性较高，还存在较多未被测序检测到的物种。下图中的稀释曲线横坐标代表随机抽取的序列数量；纵坐标代表观测到的OTU数量。样本曲线的延伸终点的横坐标位置为该样本的测序数量，如果曲线趋于平坦表明测序已趋于饱和，增加测序数据无法再找到更多的OTU；反之表明不饱和，增加数据量可以发现更多OTU。7.7Shannon-Winner曲线Shannon-Wiener曲线，是利用shannon指数来进行绘制的，反映样品中微生物多样性的指数，利用各样品的测序量在不同测序深度时的微生物多样性指数构建曲线，以此反映各样本在不同测序数量时的微生物多样性。当曲线趋向平坦时，说明测序数据量足够大，可以反映样品中绝大多数的微生物物种信息。与上图一样，横坐标代表随机抽取的序列数量；纵坐标代表的是反映物种多样性的Shannon指数。样本曲线的延伸终点的横坐标位置为该样本的测序数量，如果曲线趋于平坦表明测序已趋于饱和，增加测序数据无法再找到更多的OTU；反之表明不饱和，增加数据量可以发现更多OTU。其中曲线的最高点也就是该样本的Shannon指数，指数越高表明样品的物种多样性越高。Q:Shannon指数怎么算的？A:Shannon指数公式：其中，Sobs=实际测量出的OTU数目；ni=含有i条序列的OTU数目；N=所有的序列数。8.8Rank-Abundance曲线用于同时解释样品多样性的两个方面，即样品所含物种的丰富程度和均匀程度。物种的丰富程度由曲线在横轴上的长度来反映，曲线越宽，表示物种的组成越丰富；物种组成的均匀程度由曲线的形状来反映，曲线越平坦，表示物种组成的均匀程度越高。一般超过20个样本图就会变得非常复杂而且不美观，所以一般20个样本以下会做该图，图片保存为结果目录中rank.pdf。横坐标代表物种排序的数量；纵坐标代表观测到的相对丰度。样本曲线的延伸终点的横坐标位置为该样本的物种数量，如果曲线越平滑下降表明样本的物种多样性越高，而曲线快速陡然下降表明样本中的优势菌群所占比例很高，多样性较低。9.9Alpha多样性（样本内多样性）Alpha多样性是指一个特定区域或者生态系统内的多样性，常用的度量指标有Chao1丰富度估计量（Chao1richnessestimator）、香农-威纳多样性指数（Shannon-wienerdiversityindex）、辛普森多样性指数（Simpsondiversityindex）等。计算菌群丰度：Chao、ace；计算菌群多样性：Shannon、Simpson。Simpson指数值越大，说明群落多样性越高；Shannon指数越大，说明群落多样性越高。表中显示前10个样本，如果样本大于10个，详见结果目录中的alpha_div.txt。Q:能不能解释下每个指数（如chao1、shannon）？A:Chao1：是用chao1算法估计群落中含OTU数目的指数，chao1在生态学中常用来估计物种总数，由Chao(1984)最早提出。Chao1值越大代表物种总数越多。Schao1=Sobs+n1(n1-1)/2(n2+1)其中Schao1为估计的OTU数，Sobs为观测到的OTU数，n1为只有一条序列的OTU数目