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

微生物群落多样性测序与功能分析微生物群落多样性测序与功能分析 微生物群落测序是指对微生物群体进行高通量测序 通过分析测序序列的 构成分析特定环境中微生物群体的构成情况或基因的组成以及功能 借助不同 环境下微生物群落的构成差异分析我们可以分析微生物与环境因素或宿主之间 的关系 寻找标志性菌群或特定功能的基因 对微生物群落进行测序包括两类 一类是通过 16s rDNA 18s rDNA ITS 区域进行扩增测序分析微生物的群体构 成和多样性 还有一类是宏基因组测序 是不经过分离培养微生物 而对所有 微生物 DNA 进行测序 从而分析微生物群落构成 基因构成 挖掘有应用价值 的基因资源 以 16s rDNA 扩增进行测序分析主要用于微生物群落多样性和构成的分析 目前的生物信息学分析也可以基于 16s rDNA 的测序对微生物群落的基因构成和 代谢途径进行预测分析 大大拓展了我们对于环境微生物的微生态认知 目前我们根据 16s 的测序数据可以将微生物群落分类到种 species 一 般只能对部分菌进行种的鉴定 甚至对亚种级别进行分析 几个概念 16S rDNA 或 16S rRNA 16S rRNA 基因是编码原核生物核糖体小亚基的 基因 长度约为 1542bp 其分子大小适中 突变率小 是细菌系统分类学研究 中最常用和最有用的标志 16S rRNA 基因序列包括 9 个可变区和 10 个保守区 保守区序列反映了物种间的亲缘关系 而可变区序列则能体现物种间的差异 16S rRNA 基因测序以细菌 16S rRNA 基因测序为主 核心是研究样品中的物种分 类 物种丰度以及系统进化 OTU operational taxonomic units OTUs 在微生物的免培养分析中经常 用到 通过提取样品的总基因组 DNA 利用 16S rRNA 或 ITS 的通用引物进行 PCR 扩增 通过测序以后就可以分析样品中的微生物多样性 那怎么区分这些 不同的序列呢 这个时候就需要引入 operational taxonomic units 一般情 况下 如果序列之间 比如不同的 16S rRNA 序列的相似性高于 97 就可以把它 定义为一个 OTU 每个 OTU 对应于一个不同的 16S rRNA 序列 也就是每个 OTU 对应于一个不同的细菌 微生物 种 通过 OTU 分析 就可以知道样品中的微 生物多样性和不同微生物的丰度 测序区段 由于 16s rDNA 较长 1 5kb 我们只能对其中经常变化的区 域也就是可变区进行测序 16s rDNA 包含有 9 个可变区 分别是 v1 v9 一般 我们对 v3 v4 双可变区域进行扩增和测序 也有对 v1 v3 区进行扩增测序 工具 原料 16s rDNA 测序首先需要提取环境样品的 DNA 这些 DNA 可以来自土壤 粪便 空气或水体等任何来源 提取 DNA 后需要经过质检和纯化 一般 16s rDNA 测序扩增对 DNA 的总量 要求并不高 总量大于 100ng 浓度大于 10ng ul 一般都可以满足要求 如果 是来自和寄主共生的环境如昆虫的肠道微生物 提取时可能包括了寄主本身的 大量 DNA 对 DNA 的总量要求会提高 微生物菌群多样性测序受 DNA 提取和扩 增影响很大 不同的扩增区段和扩增引物甚至 PCR 循环数的差异都会对结果有 所影响 因而建议同一项目不同样品的都采用相同的条件和测序方法 这样相 互之间才存在可比性 完成 PCR 之后的产物一般可以直接上测序仪测序 在上机测序前我们需 要对所有样本进行定量和均一化 通常要进行荧光定量 PCR 完成定量的样品 混合后就可以上机测序 16s rDNA 测序目前可以采用多种不同的测序仪进行测序 包括罗氏的 454 Illumina 的 MiSeq Life 的 PGM 或 Pacbio 的 RSII 三代测序仪 不同的 仪器各有优缺点 目前最主流的是 Illumina 公司的 MiSeq 因为其在通量 长 度和价格三者之间最为平衡 MiSeq 测序仪可以产生 2x300bp 的测序读长 一 次可以产生 15Gb 的测序数据远远大于其他测序仪的测序通量 方法 步骤 1 1 16s rDNA 分析基本流程 2 2 原始数据处理 原始测序数据需要去除接头序列 并将双端测序序列进行拼接成单条序 列 根据测序 barcode 序列区分不同的样本序列 过滤低质量序列和无法比对到 16s rDNA 数据库的序列 3 3 OTU 分类和统计 OTU operational taxonomic units 是在系统发生学研究或群体遗传学 研究中 为了便于进行分析 人为给某一个分类单元 品系 种 属 分组 等 设置的同一标志 通常按照 97 的相似性阈值将序列划分为不同的 OTU 每一个 OTU 通常被视为一个微生物物种 相似性小于 97 就可以认为 属于不同的种 相似性小于 93 95 可以认为属于不同的属 样品中的微 生物多样性和不同微生物的丰度都是基于对 OTU 的分析 使用 QIIME version 1 8 0 工具包进行统计注释 使用 QIIME version 1 9 0 ucluster 方法根据 97 的序列相似度将所有序列进行同源比对并聚类成 operational taxonomic units OTUs 然后与数据库 GreenGenes version gg 13 8 http greengenes lbl gov cgi bin JD Tutorial nph 16S cgi 进行比对 比对方法 uclust identity 0 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 16S rRNA 总长约 1540 bp 包含 9 个可变区 由于高通量测序的测 序长度的限制 不可能将 16S rRNA 的 9 个可变区全部测序 所以在 PCR 扩增 时往往只能选择 1 3 个可变区作为扩增片段 Kozich 等评估了 Miseq 测序仪 分析的不同 16S rRNA 可变区的准确性发现 测定 V4 区效果最佳 根据我们 的测序长度 v3 v4 区是最佳选择 5 5 我们还需要对样本之间或分组之间的 OTU 进行比较获得韦恩图 注意 韦恩图目前一般最多只能显示 5 个样本或分组 过多的样本无法 无法进行韦恩图绘制 6 6 样品构成丰度 稀释曲线 微生物多样性分析中需要验证测序数据量是否足以反映样品中的物种多 样性 稀释曲线 丰富度曲线 可以用来检验这一指标 稀释曲线是用来评价测序量是否足以覆盖所有类群 并间接反映样品中 物种的丰富程度 稀释曲线是利用已测得 16S rDNA 序列中已知的各种 OTU 的 相对比例 来计算抽取 n 个 n 小于测得 reads 序列总数 reads 时出现 OTU 数量的期望值 然后根据一组 n 值 一般为一组小于总序列数的等差数列 与其相对应的 OTU 数量的期望值做出曲线来 当曲线趋于平缓或者达到平台 期时也就可以认为测序深度已经基本覆盖到样品中所有的物种 反之 则表 示样品中物种多样性较高 还存在较多未被测序检测到的物种 下图中的稀释曲线 横坐标代表随机抽取的序列数量 纵坐标代表观测到的 OTU 数量 样本 曲线的延伸终点的横坐标位置为该样本的测序数量 如果曲线趋于平坦表明 测序已趋于饱和 增加测序数据无法再找到更多的 OTU 反之表明不饱和 增加数据量可以发现更多 OTU 7 7 Shannon Winner 曲线 Shannon Wiener 曲线 是利用 shannon 指数来进行绘制的 反映样品中 微生物多样性的指数 利用各样品的测序量在不同测序深度时的微生物多样 性指数构建曲线 以此反映各样本在不同测序数量时的微生物多样性 当曲 线趋向平坦时 说明测序数据量足够大 可以反映样品中绝大多数的微生物 物种信息 与上图一样 横坐标代表随机抽取的序列数量 纵坐标代表的是反映物 种多样性的 Shannon 指数 样本曲线的延伸终点的横坐标位置为该样本的测序数量 如果曲线趋于 平坦表明测序已趋于饱和 增加测序数据无法再找到更多的 OTU 反之表明 不饱和 增加数据量可以发现更多 OTU 其中曲线的最高点也就是该样本的 Shannon 指数 指数越高表明样品的 物种多样性越高 Q Shannon 指数怎么算的 A Shannon 指数公式 其中 Sobs 实际测量出的 OTU 数目 ni 含有 i 条序列的 OTU 数目 N 所有的序列数 8 8 Rank Abundance 曲线 用于同时解释样品多样性的两个方面 即样品所含物种的丰富程度和均 匀程度 物种的丰富程度由曲线在横轴上的长度来反映 曲线越宽 表示物种的 组成越丰富 物种组成的均匀程度由曲线的形状来反映 曲线越平坦 表示物种组成 的均匀程度越高 一般超过 20 个样本图就会变得非常复杂而且不美观 所以一般 20 个样 本以下会做该图 图片保存为结果目录中 rank pdf 横坐标代表物种排序的数量 纵坐标代表观测到的相对丰度 样本曲线的延伸终点的横坐标位置为该样本的物种数量 如果曲线越平 滑下降表明样本的物种多样性越高 而曲线快速陡然下降表明样本中的优势 菌群所占比例很高 多样性较低 9 9 Alpha 多样性 样本内多样性 Alpha 多样性是指一个特定区域或者生态系统内的多样性 常用的度量 指标有 Chao1 丰富度估计量 Chao1 richness estimator 香农 威纳 多样性指数 Shannon wiener diversity index 辛普森多样性指数 Simpson diversity index 等 计算菌群丰度 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 数目 n2 为只有两条序列的 OTU 数目 Shannon 用来估算样品中微生物的多样性指数之一 它与 Simpson 多 样性指数均为常用的反映 alpha 多样性的指数 Shannon 值越大 说明群落 多样性越高 Ace 用来估计群落中含有 OTU 数目的指数 由 Chao 提出 是生态学中 估计物种总数的常用指数之一 与 Chao1 的算法不同 Simpson 用来估算样品中微生物的多样性指数之一 由 Edward Hugh Simpson 1949 提出 在生态学中常用来定量的描述一个区域的生物多样 性 Simpson 指数值越大 说明群落多样性越高 辛普森多样性指数 随机取样的两个个体属于不同种的概率 1 随机取样的两个个体属于同种的概率 10 10 Beta 多样性分析 样品间差异分析 Beta 多样性度量时空尺度上物种组成的变化 是生物多样性的重要组成 部分 与许多生态学和进化生物学问题密切相关 因此在最近 10 年间成为生 物多样性研究的热点问题之一 PCoA 分析 PCoA principal co ordinates analysis 是一种研究数据相似性或差 异性的可视化方法 通过一系列的特征值和特征向量进行排序后 选择主要 排在前几位的特征值 PCoA 可以找到距离矩阵中最主要的坐标 结果是数据 矩阵的一个旋转 它没有改变样品点之间的相互位置关系 只是改变了坐标 系统 通过 PCoA 可以观察个体或群体间的差异 每一个点代表一个样本 相同颜色的点来自同一个分组 两点之间距离 越近表明两者的群落构成差异越小 PCoA 有多张图 分别代表的 PCoA1 2 2 3 3 1 11 11 NMDS 分析 非度量多维尺度分析 NMDS Nonmetric Multidimensional Scaling 常用于比对样本组之间 的差异 可以基于进化关系或数量距离矩阵 横轴和纵轴 表示基于进化或者数量距离矩阵的数值 在二维表中成图 与 PCA 分析的主要差异在于考量了进化上的信息 每一个点代表一个样本 相同颜色的点来自同一个分组 两点之间距离 越近表明两者的群落构成差异越小 12 12 PCA 分析 主成分分析 PCA Principal component analysis 是一种研究数据相 似性或差异性的可视化方法 通过一系列的特征值和特征向量进行排序后 选择主要的前几位特征值 采取降维的思想 PCA 可以找到距离矩阵中最主 要的坐标 结果是数据矩阵的一个旋转 它没有改变样品点之间的相互位置 关系 只是改变了坐标系统 详细关于主成分分析的解释推荐大家看一篇文 章 通过 PCA 可以观察个体或群体间的差异 每一个点代表一个样本 相同颜色的点来自同一个分组 两点之间距离 越近表明两者的群落构成差异越小 以上三个图可能遇到的问题 1 PCA PcoA NMDS 分析分别是基于什么数据画的 回答 PCA PcoA NMDS 分析均是基于 OTU 分类 taxon 数据所画 用的 是 R 语言 Vegan 包中的相关函数画成 其中 PcoA 与 NMDS 还要基于样本之间 的距离矩阵才能画成 2 PCA 分析如果图中大部分点集中在一起 少数点在很远的外围 是什 么原因造成的 回答 是因为样本 OTU 分类时候 少数样本某些菌含量特别高所造成 导致这些样本偏离正常范围 建议单独拿出这些样本观察 看是否是实验错 误 3 PCA 分析时 不是有 PC1 PC2 PC3 三个坐标吗 是给出三张图吗 还是三维立体图 回答 PCA 作图时 会得出 PC1 PC2 PC3 三个坐标 可以根据 PC12 PC13 PC23 分别作图 一般给出的是 PC12 的图 当 PC12 图质量不好 看不出明显的样本分类效果时 可以看 PC13 或 PC23 的图分类是否清晰 也 可以用 R 语言 rgl 包做出 PC123 三维图 QIIME 本身结果中有提供 PCA 的三维图结果 可以通过网页打开 13 13 LDA 差异贡献分析 PCA 和 LDA 的差别在于 PCA 它所作的只是将整组数据整体映射到最方 便表示这组数据的坐标轴上 映射时没有利用任何数据内部的分类信息 是 无监督的 而 LDA 是由监督的 增加了种属之间的信息关系后 结合显著性 差异标准测试 克鲁斯卡尔 沃利斯检验和两两 Wilcoxon 测试 和线性判别分 析的方法进行特征选择 除了可以检测重要特征 他还可以根据效应值进行 功能特性排序 这些功能特性可以解释顶部的大部分生物学差异 详细说明 可以参考这篇文章 不同颜色代表不同样本或组之间的显著差异物种 使用 LefSe 软件分析 获得 其中显著差异的 logarithmic LDA score 设为 2 问题 LDA 分析有什么用 回答 组间差异显著物种又可以称作生物标记物 biomarkers 该分 析主要是想找到组间在丰度上有显著差异的物种 14 物种进化树的样本群落分布图 是将不同样本的群落构成及分布以物种分类树的形式在一个环图中展示 数据经过分析后 将物种分类树和分类丰度信息通过软件 GraPhlAn http huttenhower sph harvard edu GraPhlAn 进行绘制 其 目的是将物种之间的进化关系以及不同样本的物种分布丰度和最高分布样本 的信息在一个视觉集中的环图中一次展示 其提供的信息量较其他图最为丰 富 中间为物种进化分类树 不同颜色的分支代表不同的纲 具体的代表颜 色见右上角的图例 之后外圈的灰色标示字母的环表示的是本次研究中比 例最高的 15 个科 字母代表的科参见左上角的图例 之后的外圈提供的是 热力图 如果样本数 10 个则绘制样本 如果样本数超过 10 个则按照分组 绘制 每一环为一个样本 根据其丰度绘制的热力图 最外圈为柱状图 绘 制的是该属所占比例最高的样本的丰度和样本颜色 样本颜色见环最下方的 样本名字的颜色 其中热力图和柱状图取值均为原比例值 x10000 后进行 log2 转换后的值 参考文献 1 Vazquez Baeza Y Pirrung M Gonzalez A Knight R 2013 Emperor A tool for visualizing high throughput microbial community data Gigascience 2 1 16 2 Legendre P and Legendre L 1998 Numerical Ecology Second English Edition Developments in Environmental Modelling 20 Elsevier Amsterdam 3 Segata N Izard J Waldron L et al Metagenomic biomarker discovery and explanation J Genome Biol 2011 12 6 R60 4 Langille MGI Zaneveld J Caporaso JG McDonald D Knights D Reyes JA et al 2013 Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences Nat Biotechnol 31 814 821 15 物种相关性分析 根据各个物种在各个样品中的丰度以及变化情况 计算物种之间的相关 性 包括正相关和负相关 相关性分析使用 CCREPE 算法 首先对原始 16s 测序数据的种属数量进行 标准化 然后进行 Spearman 和 Pearson 秩相关分析并进行统计检验 计算出 各个物种之间的相关性 之后在所有物种中根据 simscore 绝对值的大小 挑 选出相关性最高的前 100 组数据 基于 Cytoscap 绘制共表达分析网络图 网 络图采用两种不同的形式表现出来 物种相关性网络图 A 图中每一个点代表一个物种 存在相关性的物种 用连线连接 其中 红色的连线代表负相关 绿色的先代表正相关 连线颜 色的深浅代表相关性的高低 物种相关性网络图 B 图中每一个点代表一个物种 点的大小表示与其 他物种的关联关系的多少 其中与之有相关性的物种数越多 点的半径和字 体越大 连线的粗细代表两物种之间相关性的大小 连线越粗 相关性越高 参考文献 Schwager E Weingart G Bielski C et al CCREPE Compositionality Corrected by Permutation and Renormalization J 2014 16 聚类分析 根据 OUT 数据进行标准化处理 1wlog10 之后 选取数目最多的前 60 个物种 基于 R heatmap 进行作图 热图中的每一个色块代表一个样品的一 个属的丰度 样品横向排列 属纵向排列 两个热图 差异是是否对样品进 行聚类 从聚类中可以了解样品之间的相似性以及属水平上的群落构成相似 性 如果聚类结果中出现大面积的白或黑是因为大量的菌含量非常低 导致 都没有数值 可以在绘制之前进行标准化操作 对每一类菌单独自身进行 Z 标准化 17 群落功能差异分析 通过对已有测序微生物基因组的基因功能的构成进行分析后 我们可以 通过 16s 测序获得的物种构成推测样本中的功能基因的构成 从而分析不同 样本和分组之间在功能上的差异 PICRUSt Nature Biotechnology 1 10 8 2013 通过对宏基因组测序数据功能分析和对应 16s 预测功能分析结果的比较 发现 此方法的准确性在 84 95 对肠道微生物菌群和土壤菌群的功能分 析接近 95 能非常好的反映样品中的功能基因构成 为了能够通过 16s 测序数据来准确的预测出功能构成 首先需要对原始 16s 测序数据的种属数量进行标准化 因为不同的种属菌包含的 16s 拷贝数 不相同 然后将 16s 的种属构成信息通过构建好的已测序基因组的种属功能 基因构成表映射获得预测的功能结果 根据属这个水平 对不同样本间的 物种丰度进行显著性差异两两检验 我们这里的检验方法使用 STAMP 中的 two sample 中 T TEST 方法 Pvalue 值过滤为 0 05 作 Extent error bar 图 此处提供 COG KO 基因预测以及 KEGG 代谢途径预测 用户也可自行使用 我们提供的文件和软件 STAMP 对不同层级以及不同分组之间进行统计分析 和制图 以及选择不同的统计方法和显著性水平 参考文献 Donovan H Parks1 Gene W Tyson STAMP statistical analysis of taxonomic and functional profiles Bioinformatics 2014 30 21 3123 3124 doi 10 1093 18 COG 构成差异分析图 图中不同颜色代表不同的分组 列出了 COG 构成在组间存在显著差异的 功能分类以及在各组的比例 此外右侧还给出了差异的比例和置信区间以及 P value 19 KEGG 代谢途径差异分析图 通过 KEGG 代谢途径的预测差异分析 我们可以了解到不同分组的样品之 间在微生物群落的功能基因在代谢途径上的差异 以及变化的高低 为我们 了解群落样本的环境适应变化的代谢过程提供一种简便快捷的方法 图解读 图中不同颜色代表不同的分组 列出了在第三层级的构成在组 间存在显著差异的 KEGG 代谢途径第三层分类以及在各组的比例 此外右侧还 给出了差异的比例和置信区间以及 P value 本例图所显示的是第三层级的 KEGG 代谢途径的差异分析 也可以针对第 二或第一层的分级进行分析 20 基因的差异分析图 除了能对大的基因功能分类和代谢途径进行预测外 我们还能提供精细 的功能基因的数量和构成的预测 以及进行样本间以及组间的差异分析 并 给出具有统计意义和置信区间的分析结果 这一分析将我们对于样本群落的差异进一步深入到了每一类基因的层面 图解读 图中不同颜色代表不同的分组 列出了在组间 样本间存在显著 差异的每一个功能基因 酶 以及在各组的比例 此外右侧还给出了差异的 比例和置信区间以及 P value 21 在获得标准报告后如果希望单独修改分组或对某些组之间进行显著 性差异分析 可以使用 STAMP 软件在自己的电脑上进行数据分析 STAMP 提 供了丰富的统计检验方法和图形化结果的输出 在使用 STAMP 之前需要首先准备需要的 spf 格式文件和样品分组信息表 在我们的报告中已经将 KEGG 和 KO 以及 COG 的结果文件后经过转换生成了适 用于 STAMP 软件打开的 spf 格式文件 还有对应的分组信息表文件 groupfile txt 以下是使用 STAMP 时的一些相关问题 详细的 STAMP 使用教程可以参考 我们提供的 STAMP 使用教程 1 stamp 作图用的原始数据的来源 STAMP 可以直接使用来自 QIIME 的 biom 文件和 PICUST 的 KEGG 和 ko 文 件 groupfile txt 文件的格式为 tab saperated value tab 键隔开的数据 2 分组问题 导入数据之后 view group legend 在窗口右侧会出现 分组栏 根据需要进行分组 3 Unclassiffied 选项中 remain Unclassiffied reads remove Unclassiffied reads 和 use only for calculating frequency profiles 方法的区别 remain Unclassiffied reads 和 use only for calculating frequency profiles 方法会保留所有的数据 而 r