诺禾致源宏转录组报告诺禾致源宏转录组报告 北京诺禾致源生物信息科技有限公司 I 目录目录 1 概述概述1 2 项目流程项目流程2 2.1 实验上机流程 2 2.1.1Total RNA 样品检测.3 2.1.2文库构建及库检. 3 2.1.3上机测序. 4 2.2 信息分析流程. 4 3分析及结果分析及结果. 6 3.1数据预处理 6 3.2De novo 组装7 3.3 物种注释 8 3.4功能注释 9 3.4.1eggNOG/COG 注释 10 3.4.2KEGG 注释. 11 3.4.3CAZy 注释 12 3.5 基因表达水平分析. 13 3.5.1 参考序列比对 13 3.5.2 基因表达水平统计表 14 3.5.3基因表达差异分析. 15 3.5.4 差异基因 GO 富集分析15 3.5.5 差异基因 KEGG 富集分析.17 3.6多样品之间的比较分析 19 3.6.1多样品间 eggNOG/KEGG/CaZy 功能比较.20 3.6.2多样品间功能聚类. 21 II 3.6.3多样品间功能的 PCoA 分析25 4 参考资料参考资料. 27 1 1 概述概述 在地球生物圈中，微生物扮演着极为重要的角色，它们的活动影响着自然环境的营养 循环， 土壤肥力， 有机质的分解， 以及物种与能量之间的交换。 人类对微生物的研究从Antoni van Leeuwenhoek 发明显微镜开始的数百年中，主要基于纯培养的研究方式，而在数以万 亿计的微生物种类中，仅 0.1%1%的物种可培养，极大地限制了对微生物多样性资源的研 究和开发。 宏转录组学(Metatranscriptomics)兴起于宏基因组之后， 从整体水平上研究某一特定环 境，特定时期群体生命全部基因组转录情况以及转录调控规律，它以生态环境中的全部 RNA 为研究对象，避开了微生物分离培养困难的问题，能有效的扩展微生物资源的利用空 间。2006 年，Leiniger 等首次使用 454 测序技术对一个复杂微生物群落的宏转录组进行研 究。与宏基因组学相比较，宏转录组学能从转录水平研究复杂微生物群落变化，能更好的 挖掘潜在的新基因。 近年来，随着测序技术和信息技术的快速发展，利用新一代测序技术(Next Generation Sequencing)研究宏转录组，能快速准确的得到大量生物数据和丰富的微生物研究信息，从 而成为研究微生物多样性和群落特征的重要手段。如致力于研究微生物与人类疾病健康关 系的人体微生物组计划(HMP, Human Microbiome Project, http:/www.hmpdacc.org/ )，研究 全球微生物组成和分布的全球微生物组计划(EMP, Earth Microbiome Project, http:/www.earthmicrobiome.org/ )都主要利用高通量测序技术进行研究。 2 2 项目流程项目流程 2.1 实验上机流程 从环境（如土壤、海洋、淡水、肠道等）中采集实验样品，进行必要的处理后，将样 品以适当形式（如原始采样样品、已提取的 RNA 样品等）送往我公司。我公司将对接收 到的样品进行必要的预实验处理，执行严格的样品质控。 检测合格的样品，经过片段筛选、建库并做相应的实验检测。检测合格的文库将采用 Illumina Hiseq 高通量测序平台进行测序， 测序得到的下机数据(Raw Data)将用于后期信息 分析。 图图 1 宏转录组样品的实验上机流程宏转录组样品的实验上机流程 3 2.1.1Total RNA 样品检测样品检测 诺禾致源对 RNA 样品的检测主要包括 4 种方法： (1) 琼脂糖凝胶电泳分析 RNA 降解程度以及是否有污染 (2) Nanodrop 检测 RNA 的纯度(OD260/280 比值) (3) Qubit 对 RNA 浓度进行精确定量 (4) Agilent 2100 精确检测 RNA 的完整性 2.1.2文库构建及库检文库构建及库检 样品检测合格后，通过试剂盒去除 rRNA。随后加入 fragmentation buffer 将 mRNA 打 断成短片段，以 mRNA 为模板，用六碱基随机引物(random hexamers)合成一链 cDNA，然 后加入缓冲液、dNTPs 和 DNA polymerase I 和 RNase H 合成二链 cDNA， 再用 AMPure XP beads 纯化双链 cDNA。纯化的双链 cDNA 先进行末端修复、加 A 尾并连接测序接头，再 用 AMPure XP beads 进行片段大小选择。最后进行 PCR 扩增， 并用 AMPure XP beads 纯化 PCR 产物，得到最终的文库。文库构建完成后，先使用 Qubit2.0 进行初步定量，稀释文库 至 2ng/ul，随后使用 Agilent 2100 对文库的 insert size 进行检测，insert size 符合预期后，使 用 Q-PCR 方法对文库的有效浓度进行准确定量(文库有效浓度 2nM)， 以保证文库质量。 文库构建原理图如下： 4 图图 2文库构建流程文库构建流程 2.1.3上机测序上机测序 库检合格后，把不同文库按照有效浓度及目标下机数据量的需求 pooling 后进行 Illumina HiSeq 测序。 2.2 信息分析流程信息分析流程 测序得到的原始数据(Raw Data)，会存在一定比例的低质量数据，为了保证后续信息 分析结果的准确可靠，首先要对原始数据进行预处理，得到有效数据(Clean Data)。 5 然后基于有效数据进行物种分类分析和复杂度分析以及基因的表达丰度分析；再进行 拼接与组装，进行代谢通路(KEGG)，同源基因簇(eggNOG)，碳水化合物酶(CAZy)等功 能注释，全面了解样品中的微生物组成结构和功能注释信息。 最后，基于以上分析结果，可以进行多样品比较分析，如聚类分析，PCoA 分析等， 挖掘出样品之间的物种和功能差异。 图图 3宏转录组分析流程宏转录组分析流程 6 3分析及结果分析及结果 3.1数据预处理数据预处理 测序产生的原始数据（Raw Data）存在一定比例低质量数据，为了保证后续分析的结 果准确可靠，首先对原始的测序数据进行预处理，获取用于后续分析的有效数据（Clean Data）。预处理方法参见方法。处理步骤如下： 1) 去除质量值5 的碱基数达到一定比例的 reads （默认 reads 长度的 40%， 设置为 40） ； 2) 去除含 N 的碱基数目达到一定比例的 reads（默认 reads 长度的 10%，设置为 10）； 3) 去除 Adapter 污染 （默认 Adapter 序列与 reads 序列有 15 bp 的 overlap， 设置为 15） ； 4) 在有宿主污染可能性的前提下，需与宿主数据库进行比对，过滤掉可能是宿主污 染的 reads（默认设置比对一致性90%的 reads 为宿主污染）。测序数据处理结果见表 1 表表 1 数据预处理统计表数据预处理统计表 SampleRaw ReadsClean ReadsClean BasesError (%)Q20(%)Q30(%)GC(%) RL1_121221813190864512.86G0.0298.6996.0252.03 RL1_221221813190864512.86G0.0495.1987.5052.21 RL2_121580733191265322.87G0.0298.7596.1652.89 RL2_221580733191265322.87G0.0494.6686.3153.11 RL3_123504202231447053.47G0.0298.5096.0352.12 RL3_223504202231447053.47G0.0494.2887.1852.31 Sample：样品名。 Raw reads：统计原始序列数据，以四行为一个单位，统计每个文件的测序序列的个数。 Clean reads：计算方法同 Raw Reads、Raw bases，只是统计的文件为过滤后的测序数据。后续的生物信 息分析都是基于 Clean reads。 Clean bases：测序序列的个数乘以测序序列的长度，并转化为以 G 为单位。 Error rate：碱基错误率。 Q20、Q30：Phred 数值大于 20、30 的碱基占总体碱基的百分比。 GC content：碱基 G 和 C 的数量总和占总的碱基数量的百分比。 7 3.2De novo 组装组装 针对每个样品经预处理得到的 Clean Reads， 先使用 NCBI 的 rRNA、 tRNA 以及 SILVA 数据库进行比对分离出来宏基因组中 rRNA 序列，剩下的 mRNA 序列则使用拼接软件 Trinity(version: r20140413p1)分别进行从头组装，然后对所有样品的序列整合并使用 CD-HIT-EST 去冗余（设定序列一致性阈值为 0.95），得到 unigene 集合。 表表 2 组装结果统计组装结果统计 AssemblyStatistics Number of Unigene493186 Large Unigene(=1000bp)49938 Max Unigene length(bp)46320 Mean Unigene length(bp)596 N50 length(bp)650 图图 4组装组装 unigene 长度分布统计长度分布统计 8 3.3 物种注释物种注释 通过与 Nr （ NCBI non-redundant protein sequences）库进行 BLAST 比对（evalue 1e-5），由于每一条序列可能会有多个比对结果，得到多个不同的分类级别，为了保证其 生物意义，采取 LCA 算法(应用 MEGAN 软件的系统分类(Huson, Daniel H., et al,2011)， 将出现第一个分支前的分类级别，作为该序列的物种注释信息。从门水平上的相对丰度表 出发，选取出在各样品中的最大相对丰度排名前 10 的门类，并将其余的物种设置为 Others，绘制出各样品对应的物种注释结果在门水平的统计图。 图图 5物种注释结果在门水平的统计图物种注释结果在门水平的统计图 纵坐标为注释到某类型的物种的相对比例； 横坐标为样品名称； 各颜色区块对应的物种类别见右侧 图例。 根据所有样品在属水平的物种注释及丰度信息，选取丰度排名前 35 的属及它们在每 个样品中的丰度信息绘制热图，并从分类信息和样品间差异两个层面进行聚类，方便而找 出研究样品中聚集较多的物种或样品，结果展示见图 6。 9 图图 6物种丰度聚类图物种丰度聚类图 横坐标为样品信息，纵坐标为物种注释信息，图中左侧为物种聚类树；上方为样品聚类树；中间热 图对应的值为每一行物种相对丰度经过标准化处理后 得到的 Z 值，即一个样品在某个分类上的 Z 值 为样品在该分类上的相对丰度和所有样品在该分类的平均相对丰度的差除以所有样品在该分类上的标 准差所得到的值。 3.4功能注释功能注释 组装得到的转录组与不同功能的注释数据库进行蛋白序列BLAST比对进行功能注释。 由于每一条序列比对结果可能不止一条，为保证后续研究的生物意义，再从每条序列的比 对结果进行筛选，并且计算 Reference 与 Query 中每个基因的覆盖比率 BCR（The BLAST 10 Coverage Ratio），保证每条比对记录中的 BCR (Ref.)和 BCR (Que.)大于 40%，然后根据各 数据库的特点统计汇总，最终得到对应的功能注释信息。 Reference 与 Query 基因的 BCR 值计算公式如下，其中 Match 为二者比对有效长度， Length (R)为 Reference 基因长度，Length (Q)为 Query 基因长度。 BCR (Ref.