宏基因组测序.pdf

1、 宏 基因组 测序 结题报告 项目名称 : XXXXXX 客户 单位 : XXXXXX 报告单位 : XXXXXX 报告日期 : XXXXXX Metagenome 目 录 一、 项目信息 1 二、 项目流程 2 2.1 项目执行流程 2 2.1.1 DNA提取及检测 2 2.1.2 文库构建及质检 . 2 2.1.3 Hiseq测序 . 2 2.2 生物信息分析流程 3 2.3 项目分析内容 4 三、 项目结果报告 5 3.1 测序数据预处理 5 3.1.1 测序数据过滤与拼接 . 5 3.2 基因预测及丰度分析 6 3.2.1 基因预测 . 6 3.2.2 非冗余基因集的构建 . 7 3.

2、2.3 Core-pan基因分析 8 3.2.4 基因数目 差异分析 . 9 3.3 物种注释分析 .11 3.3.1 物种注释 11 3.3.2 物种相对丰度 . 13 3.3.3 优势物种分析 . 15 3.4 功能注释分析 16 Metagenome 3.4.1 eggNOG功能注释 16 3.4.2 KEGG注释 . 17 3.4.3 基于 CAZy数据库功能 注释（高级分析） . 19 3.4.4 基于 ARDB数据库功能注释（高级分析） . 21 3.5组间差异显著性分析 . 23 3.5.1 PCA分析 . 23 3.5.2 Anosim分析 24 3.5.3 NMDS分析（高级

3、分析） . 25 3.5.4 PCoA分析 . 26 3.5.5 Test . 27 3.5.6 STAMP差异分析（高级分析） . 29 3.5.7 LEfSe分析 30 3.5.8 Random forest分析（高级分析） . 31 3.6 环境因子关联分析 33 3.6.1 RDA/CCA分析 . 33 3.6.2 Cytoscape网络分析（高级分析） 34 四、 参考文献 35 1 Metagenome 一、 项目信息 合同 编号 项目 名称 项目 样品信息 样品 来源 样品 数量 客户信息 单位名称 单位 地址 课题组 电话 邮箱 项目联系人 电话 邮箱 项目 信息 项目 负责人

4、 电话 邮箱 项目 整体说明 2 Metagenome 二、 项目流程 2.1 项目 执行 流程 2.1.1 DNA提取 及 检测 基因组 DNA提取采用 公司 自主研发的 meta提取流程 完成， 其中微生物 细胞破碎采用了 玻璃珠 研磨 结合化学 裂解及酶解等 方法，最大程度的还原微生物的 组成 。提取 完成后利用 1%琼脂糖凝胶电泳 检测抽提 的基因组 DNA 的 完整性，并用Nanodrop或 Qubit对 DNA的 浓 度进行精确定量 。 2.1.2 文库构建 及质检 检测合格的 DNA样品用 Covaris超声波破碎仪随机打断成长度约 为 300-400bp的片段 ，经末端修复、加

5、 A尾、加测序接头、 PCR扩增等步骤完成整个文库制备。 文库构建完成后，先使用 Qubit 进行初步定量，稀释文库至 2ng/ul，随后使用Agilent 2100对文库的 insert size进行检测， insert size符合预期后，使用 qPCR方法对文库的有效浓度进行准确定量，以保证文库质量。 2.1.3 Hiseq测序 库检合格后，把不同文库按照有效浓度及目标下机数据量的需求 pooling 后进行 Illumina HiSeq X Ten测序 。 3 Metagenome 2.2 生物 信息分析流程 首先对原始数据 (Raw Data)进行质控及宿主过滤，得到有效数据 (Cl

6、ean Data)。然后进行 Metagenome组装，并将各样品未被利用上的 reads放在一起进行混合组装，以期发现样品中的低丰度物种信息。继而进行基因预测、物种注释、常用 功能数据库注释、 物种分类分析。另外，还可以基于标准分析结果，进行一系列高级信息分析；同时，结合环境因子、病理指标或特殊表型进行深入关联研究。 原始数据数据质控个体组装 、 混合组装物种注释 功能注释K EG Ge g g N O GC A Zy （ 高级分析 ）A D R B （ 高级分析 ）组间差异显著性分析P C A 分析 S TA M P 分析P C o A 分析A n o s i m 分析N M D S 分析

7、肠型分析Te s tLE f S e 分析R a n d o m fo r e s t组装结果基因预测非冗余基因集数据质控高质量数据组装基因预测M e ta P h l A nM LG环境因子关联分析C C A / R D A物种 因子相关系数分析功能 因子相关系数分析物种 功能相关系数分析相关性网络 （ 高级分析 ）高级分析图 1 生物 信息分析流程 4 Metagenome 2.3 项目 分析 内容 本次 实验生物信息学分析内容见下表： 基本 分析 标准 分析 高级 分析 数据统计 COG功能 注释 CAZy数据库功能注释 数据质控 KEGG功能注释 ARDB数据库 功能注释 序列拼接组装

8、 物种 /功能 PCA分析 物种 /功能 STAMP分析 基因 预测 物种 /功能 PCoA分析 物种 /功能 NMDS分析 非 冗余基因集构建 物种 /功能 Anosim分析 Random forest Core-pan基因分析 Test Cytoscape网络分析 基因 数目差异分析 LEfSe分析 物种相对 丰度 物种 /功能 CCA/RDA分析 优势 物种分析 5 Metagenome 三、 项目 结果报告 3.1 测序 数据预处理 3.1.1 测序 数据过滤与拼接 Illumina Hiseq Paired-end测序 获得的原始序列 3端带有 adapter接头序列，以及一些少量低

9、质量序列和杂质序列，为了提高后续分析质量和可靠性，对原始序列进行去接头、质量剪切以及去除宿主 DNA污染等。 数据过滤 方法： a. 去除 3端 adapter，保留 reads长度大于等于 50bp、平均质量值大于等于 20的 reads； b. 去除含 N碱基的 reads，对序列 3端和 5端进行质量剪切； c. 去除剪切后长度小于 50bp及平均质量值低于 20的 reads，保留高质量 reads； d. 如果样品来源于宿主（比如人或动物的粪便），而且该宿主本身的基因组已被发表，则去除可能存在的宿主 DNA污染的序列。 软件平台： Seqprep（ https:/ Sickle（ h

10、ttps:/ BWA（ http:/bio-；去宿主污染） 。 数据拼接 方法 参考文献 ： Qin J, Li R, Raes J, et al. A human gut microbial gene catalogue established by metagenomic sequencingJ. Nature, 2010, 464(7285): 59-65. Qin J, Li Y, Cai Z, et al. A metagenome-wide association study of gut microbiota in type 2 diabetesJ. Nature, 2012,

11、490(7418): 55-60. 6 Metagenome 3.2 基因预测 及 丰度分析 3.2.1 基因预测 使用 MetaGene和 GeneMark对 拼接结果中的 contig进行 ORF预测。选择长度大于等于 100bp的基因，并将其翻译为氨基酸序列。 表 1 ORF的预测 结果 统计 Sample ORFs Total length （ bp） Average length （ bp） Max （ bp） Min （ bp） XXX XXX XXX XXX XXX XXX 注 ： Sample：样品名称； ORFs： ORFs 的序列条数； Total length： ORFs

12、 的总序列长度； Average length： ORFs 的平均序列长度； Max：最长 ORF 的序列长度； Min：最短 ORF 的序列长度。 7 Metagenome 3.2.2 非冗余基因集的构建 将所有样品预测出来的基因序列，用 CD-HIT软件进行聚类，每个类取最长的基因作为代表序列，构建非冗余基因集。 将样品所有的高质量 reads 与非冗余基因集进行比对，统计基因在对应样品中的丰度信息，从而构建 geneprofile。比对软件： SOAPaligner。 表 2 基因 在 各个样品中 的 丰度 统计 Sample1 Sample 2 Sample 3 Sample 4 Ge

13、ne1 XXX XXX XXX XXX Gene2 XXX XXX XXX XXX 注：第一列是基因名称，之后各列的数值为基因在各个样品中的丰度值。 此表格 只显示部分结果，详见附件。 8 Metagenome 3.2.3 Core-pan基因分析 从基因在各样品中的丰度表出发，可以获得各样品的基因数目信息，通过随机抽取不同数目的样品，可以获得不同数目样品组 合间的基因数目，由此构建和绘制Core和 Pan基因的稀释曲线，图片展示如下： 图 2 Core-pan基因稀释曲线 注： 上图： Core基因稀释曲线； 下图： Pan基因稀释曲线。横坐标表示抽取的样品个数；纵坐标表示抽取的样品组合的基

14、因数目。 9 Metagenome 3.2.4 基因数目差异分析 为了考察组与组间的基因数目差异情况，绘制了组间基因数目差异箱图，展示结果如下： 图 3 组间基因数目差异箱图 注 ：横坐标为各个分组信息；纵坐标为基因数目。 10 Metagenome 为了考察指定样品（组）间的基因数目分布情况，分析不同样品（组）之间的基因共有、特有信息，绘制了韦恩图 (Venn Graph)或花瓣图，展示结果如下： 图 4 基因数目韦恩图（花瓣图）分析 注 ： 当样本（组）数小于 5时，展示韦恩图，当样本（组）数超过 5个时，展示花瓣图 。 图中，每个圈代表一个样品；圈和圈重叠部分的数字代表样品之间共有的基因

15、个数；没有重叠部分的数字代表样品的特有基因个数 。 11 Metagenome 3.3 物种注释 分析 3.3.1 物种 注释 对于已知 物种 的 注释 采用 Metaphlan方法 ， 示意 结果 图 如下 ， 参考文献 ： Segata N, Waldron L, Ballarini A, et al. Metagenomic microbial community profiling using unique clade-specific marker genesJ. Nature methods, 2012, 9(8): 811-814. 图 5 西方人群 中 肠道菌群的 物种 进化树

16、 注 ： HMP和 MetaHIT组 共 224个样品 中 所有 物种 的 进化枝。 圆 的 尺寸与平均丰度对数成比例 ； 颜色 深浅 代表 HMP（ n = 139）和 MetaHIT（ n = 85，仅健康人群）群体中最丰富 分类群的相对富集 程度 。 对于未知 物种 的 注释 采用 MLG方法 ， 示意 结果 图 如下， 参考文献 ： Qin J, Li R, Raes J, et al. A human gut microbial gene catalogue established by metagenomic sequencingJ. Nature, 2010, 464(7285)

17、: 59-65. Nielsen H B, Almeida M, Juncker A S, et al. Identification and assembly of genomes 12 Metagenome and genetic elements in complex metagenomic samples without using reference genomesJ. Nature biotechnology, 2014, 32(8): 822-828. 图 6 基因簇聚类 和 MGS组装概览 注 ： 提取 人类肠道微生物群 中独立 样品的 DNA，进行 鸟枪法测序。然后将单个样品

18、中组装和鉴定的基因整合形成 非冗余基因 集 。通过 对各样品中比对 上的 序列读数 计数， 评估每个基因的丰度分布。如果含有700个或更多个基因，则 该基因簇 被命名为宏基因组物种（ MGS）。 基因数较少 的组仍被称为 CAG。然后提取 比对 到 MGS的 基因及其重叠群的单个样品的序列读数， 用于组装 MGS的基因组 。 13 Metagenome 3.3.2 物种 相对丰度 在 Domain（域）、 Kingdom（界）、 Phylum（门）、 Class（纲）、 Order（目）、Family（科）、 Genus（属）、 Species（种）各个分类 学水平上统计物种在各个样品中的丰度

19、，从而构建相应分类学水平上的 taxonomy profile。 表 3 物种 在各个 样品 中的丰度 Sample 1 Sample 2 Sample 3 Sample 4 Sample 5 Taxon1 XXX XXX XXX XXX XXX Taxon1 XXX XXX XXX XXX XXX Taxon1 XXX XXX XXX XXX XXX 注：第一列是物种名称，之后各列的数值为物种在各个样品中的丰度值。 图 7 单 样品 物种 /所有 样本物种平均值 分布 饼图 14 Metagenome 图 8 所有样品 物种分布 条形图 注：为使视图效果最佳，作图时可将丰度低于 1%的部分合

20、并为 others 在图中显示。 15 Metagenome 3.3.3 优势 物种 分析 丰度 比较大的 优势 物种 可以 反映样本本身的特性或是 样品 与周围环境之间的关系。根据 物种 丰度， 挑选 在门水平 top5、 属水平 top15和 在种水平 top20的 优势物种进行分析。 图 9种 水平的优势 物种 注 ： 横坐标 代表物种名称，纵坐标代表物种丰度。 16 Metagenome 3.4 功能注释 分析 3.4.1 eggNOG功能 注释 eggNOG（ evolutionary genealogy of genes: Non-supervised Orthologous Gr

21、oups，http:/eggnog.embl.de/） 数据库，源于无指导同源群的完整基因组，基于多样化的注释，通过识别基因的共性获得详细的功能描述信息，包括来自原始 COG/KOG 的功能分类，并且提供了基于分类学水平的同源群功能注释。 eggNOG 数据库跟踪了基因组计划的完成情况，同时优先选择高质量的基因组数据来最小化由不完整蛋白质组数据集导致的错误。目前该数据库（ v4.0）包含 170 万个直系同源类群，覆盖了 3686 个物种，给定了 107 个不同的分类级别的同源群。 通过与 eggNOG 数据库进行 blastp 比对，可以获得基因所对应的 COG 注释，即得出 ORF 对应的

22、 NOG number，并根据 NOG number 进行功能归类。 通过将基因序列与 eggNOG数据库比对，统计 COG功能在各个样品中的丰度，从而构建 COG function profile。 图 10 COG功能 分类 17 Metagenome 3.4.2 KEGG注释 KEGG（ Kyoto Encyclopedia of Genes and Genomes，京都基因和基因组百科全书， http:/www.genome.jp/kegg/）是基因组研究方面的公共数据库。 该数据库是系统分析基因功能，联系基因组信息和功能信息的大型知识库。KEGG 将从 NCBI 等数据库中获得的，包

23、括完整和部分测序的基因组序列及其基因序列存储于 KEGG GENES 数据库中；将各种生物学通路信息存储在 PATHWAY 数据库中，生物学通路包括各种代谢通路、合成通路、膜转运、信号传递、细胞周期以及疾病相关通路等。此外 KEGG LIGAND 数据库中收集了各种化学分子、酶及以及酶促反应等相关信息。 在生物体内，代谢产物并不是孤立存在并各自发挥作用的，不同代谢产物之间通过有序的相互协调来一起行使具体的生物学功能。因此， KEGG数据库中丰富的通路信息将有助于我们从系统水平去了解基因的生物学功能，例如代谢途径、遗传信息传递以及细胞过程等一些复杂的生物过程。 运用 序列同源性 比对 将基因集与

24、 KEGG 的基因数据库（ GENES）进行比对，生成各 样本 ko, modul, pathway基因 丰度谱， 然后进行后续分析。 表 4 KEGG丰度 统计 Sample1 Sample 2 Sample 3 Sample 4 KO1 XXX XXX XXX XXX KO2 XXX XXX XXX XXX 注：第一列是 KO号 ，之后各列的数值为 KO在各个样品中的丰度值。 此表格 只显示部分结果，详见附件。 18 Metagenome 图 11 KEGG代谢通路 注释 注：矩形框表示催化反应进行的酶（红色边框：此次测序预测出的酶），框内编号即 EC 编号；圆圈表示代谢产物（酶促反应的反

25、应物或产物）；实线箭头的方向表示酶促反应进行的方向；虚线箭头表示此产物可再通过中间产物与其它代谢途径发生关系；圆角矩形框代表其它代谢途径。 19 Metagenome 3.4.3 基于 CAZy数据库 功能 注释 （高级分析） 碳水化合物活性酶（ Carbohydrate-active enzymes， CAZyme）对地球上所有碳水化合物的合成、降解与修饰起重要作用，因此深入研究 CAZyme，对于了解微生物碳水化合物的代谢机制非常重要。 碳水化合物活性酶数据库（ CAZy， http:/www.cazy.org/）是关于能够合成或者分解复杂碳水化合物和糖复合物的酶类的专业数据库，根据蛋白质

26、结构域中氨基酸序列的相似性，可以将不同物种来源的碳水化合物活性酶分成糖苷水解酶 （ Glycoside Hydrolases ， GHs） ， 糖基转移酶 （ GlycosylTransferases， GTs） ， 多糖裂合酶 （ Polysaccharide Lyases， PLs） ， 碳水化合物酯酶 （ Carbohydrate Esterases， CEs） ， 碳水化合物结合模块 （ Carbohydrate-Binding Modules， CBMs），碳水化合物酯酶（ Auxiliary Activities， AAs）等六大类蛋白质家族。 通过将基因序列与 CAZy 数据库进

27、行比对，可以获得碳水化合物酶类的物种来源、酶功能 EC 分类、基因序列、蛋白质序列及其结构等 信息。 表 5 CAZy 功能 在 各个样品中的丰度 #Category Description Sample1 Sample 2 Sample 3 Sample 4 Sample 5 1 Glycoside Hydrolases （ GH ） XXX XXX XXX XXX XXX 2 Glycosyl Transferases （ GT） XXX XXX XXX XXX XXX 3 Polysaccharide Lyases （ PL） XXX XXX XXX XXX XXX 4 Carbohyd

28、rate Esterases （ CE） XXX XXX XXX XXX XXX 5 Auxiliary Activities （ AA ） XXX XXX XXX XXX XXX 6 CarbohydrateBinding Modules （ CBM） XXX XXX XXX XXX XXX 注：第一列是分类编号，第二列是 CAZy 功能分类的描述，之后各列的数值为 CAZy 功能分类在各个样品中的丰度值。 20 Metagenome 图 12 CAZy功能分类 统计 图 注：横坐标为碳水化合物活性酶分类，纵坐标为相应酶分类的基因数目。 21 Metagenome 3.4.4 基于 ARDB

29、数据库 功能 注释 （高级分析） 抗生素在有害微生物的防治中发挥着重要作用，然而，随着抗生素的大量滥用，微生物通过基因的水平转移或自身的基因变异会对抗生素产生耐受性或抗药性。目前，抗生素抗性（ antibiotic resistance）已经成为全世界公共卫生的重要威胁。 抗生素抗性基因数据库（ ARDB， http:/ardb.cbcb.umd.edu/）收集了不同环境来源的（如肠道、生活废水、河流等）细菌抗药性基因及其抗性谱、作用机制、本体论、 COG 和 CDD 等注释信息，为研究药物作用、环境治理提供研究依据。当前，ARDB 包括了 13,293 个基因、 377 种类型、 257 种

30、抗生素、 632 个基因组、 933 个物种和 124 个属的抗性信息。 通过将环境微生物的序列与 ARDB 数据库进行比对分析，不仅可以获得抗生素抗性基因的种类和数量，也可以获知其功能信息。 ARDB 功能在各个样品中的丰度表如下： 表 6 AR功能 在各样品 中的丰度 #Category Description Sample1 Sample 2 Sample 3 Sample 4 Sample 5 acrb Resistance-nodulation-cell division transporter system. Multidrug resistance efflux pump. XX

31、X XXX XXX XXX XXX baca Undecaprenyl pyrophosphate phosphatase, which consists in the sequestration of Undecaprenyl pyrophosphate. XXX XXX XXX XXX XXX 注：第一列是抗性基因编号， 第二列是抗性基因的功能描述， 之后各列的数值为抗性基因在各个样品中的丰度值。 22 Metagenome 图 13 抗性 基因 功能分类 统计 图 注 ：横坐标为抗生素抗性基因，纵坐标为 相应 的基因数目。 23 Metagenome 3.5组间 差异显著性 分析 3.5

32、.1 PCA分析 PCA 分析（ Principal Component Analysis），即主成分分析，是一种对数据进行简化分析的技术 ， 这种方法可以有效的找出数据中最 “主要 ”的元素和结构 ， 去除噪音和冗余，将原有的复杂数据降维，揭示隐藏在复杂数据背后的简单结构。 通过分析不同样品物种 、 基因、 COG 和 KEGG 代谢功能 等组分，可以反映样品间的差异和距离， PCA 运用方差分解，将多组数据的差异反映在二维坐标图上，坐标轴取能够最大反映方差值的两个特征值。如样品组成越相似，反映在 PCA 图中的距离越近。不同环境间的样品可能表现出分散和聚集的分布情况， PCA 结果中对样品

33、差异性解释度最高的两个或三个成分 可以用于对假设因素进行验证。 PCA分析 ，可以在物种、基因、 COG和 KEGG代谢 功能等 层次 上 反映 不同样品或不同组间的 差异性 。 软件： R 语言 PCA 统计分析和作图。 图 14 PCA分析 注：不同颜色或形状的点代表不同环境或条件下的样本组，横、纵坐标轴的刻度是相对距离，无实际意义。 PC1、 PC2 分别代表对于两组样本微生物组成发生偏移的疑似影响因素， 需要结合样本特征信息归纳总结。 24 Metagenome 3.5.2 Anosim分析 相似性分析（ Anosim）是一种非参数检验，用来检验组间（两组或多组）的差异是否显著大于组内

34、差异，从而判断分组是否有意义。 通过 Anosim分析，可以判断样品在物种 、 基因、 COG或 KEGG代谢 通路上的差异。 软件： R 语言 vegan 包或 QIIME 软件 。 图 15 组间 与组内距离箱线图 注： 横坐标表示所有样品 (Between)以及每个分组 (A、 B)，纵坐标表示 unifrac距离的秩。 Between组比其它每个分组的秩较高时，则表明组间差异大于组内差异。 R介于 (-1， 1)之间， R大于 0，说明组间差异显著； R小于 0，说明组内差异大于组间差异，统计分析的可信度用 P表示， P 0.05表示统计具有显著性。 25 Metagenome 3.5

35、.3 NMDS分析 （高级分析） 非度量多维尺度分析（ NMDS）是一种将多维空间的研究对象（样本或变量）简化到低维空间进行定位、分析和归类，同时又保留对象间原始关系的数据分析方法。适用于无法获得研究对象间精确的相似性或相异性数据，仅能得到他们之间等级关系数据的情形。其基本特征是将对象间的相似性或相异性数据看成点间距离的单调函数，在保持原始数据次序关系的基础上，用新的相同次序的数据 列替换原始数据进行度量型多维尺度分析。 NMDS可以 对不同 样品在物种、基因、 COG和 KEGG代谢 功能等 层次 上 进行分析。 软件： R 语言 vegan 软件包进行 NMDS 分析及作图。 图 16 N

36、MDS分析 注：不同颜色或形状的点代表不同环境或条件下的样本组。 26 Metagenome 3.5.4 PCoA分析 PCoA（ principal co-ordinates analysis）是一种研究数据相似性或差异性的可视化方法，通过一系列的特征值和特征向量进行排序后，选择主要排在前几位的特征值， PCoA 可以找到距离矩阵中最主要的坐标，结果是数据矩阵的一个旋转，它没有改变样品点之间的相互位置关系，只是改变了坐标系统。通过 PCoA可以观察个体或群体间的差异。 利用 PCoA， 可以 对不同 样品在物种、基因、 COG和 KEGG代谢 功能等 层次 上进行 分析。 软件： R 语言

37、PCoA 分析和作 PCoA 图。 图 17 PCoA分析 注： PC1 和 PC2 是两个主坐标成分， PC1 表示尽可能最大解释数据变化的主坐标成分， PC2 为解释余下的变化度中占比例最大的主坐标成分， PC3 等依次类推。 27 Metagenome 3.5.5 Test 分析多个样品时，如果这些样品分属于两个组，则可以进行 Metastats 分析或T-Test分析 。 Metastats分析通过对比两组条件下的多个样品，找出两组中具有显著差异的微生物类型，默认阈值是 0.05，若表格中无相关信息则表示没有显著差异的物种。 T-Test分析，亦称 student t检验（ Stude

38、nts t test），主要用于样本 含量较小（例如 n30），总体标准差未知的正态分布。用 t 分布理论来推论差异发生的概率，从而比较两个平均数的差异是否显著。 T-Test 要求样本 正态分布 ，应用广泛； Metastat 可以不是 正态分布 ，并且会根据样本情况自动调整统计方法。 Metastat 方法 q 值筛选更 严格 ，所以找出的结果 相对 少，建议物种多的情况用 Metastat， 但是如果 Metastat 的 结果并不符合研究，可以参照 T-Test的结果。 表 7 Test分析 结果 注： ID： 菌种 ； mean：均值； p value (an individual measure of the false positive rate) 假阳性概率值，是统计学中常用的判定值，一般来说 p value0.05时差异 显著 ； q value (an individual measurement of the false discovery rate) 假发现率评估值，指本次计算可信度 ； fdr（ false discovery rate） 错误 发现率 。