转录组ref流程工作手册.docx

1、转录组ref流程工作手册转录组ref流程工作手册一、Reference 流程生物学原理1.1 实验流程图一：转录组实验流程当我们得到样品时，必须对其测序，才能得到分析所需的数据。测序基本过程：提取样品总RNA后，用带有Oligo(dT)的磁珠富集真核生物mRNA（若为原核生物，则用试剂盒去除rRNA后进入下一步）。加入fragmentation buffer将mRNA打断成短片段，以mRNA为模板，用六碱基随机引物（random hexamers）合成第一条cDNA链，然后加入缓冲液、dNTPs、RNase H 和DNA polymerase I合成第二条cDNA链，在经过QiaQuick P

2、CR试剂盒纯化并加EB缓冲液洗脱之后做末端修复并连接测序接头，然后用琼脂糖凝胶电泳进行片段大小选择，最后进行PCR扩增，使用建好的测序文库进行测序。得到RNA的序列后，又可以找到它的参考序列（物种本身的基因、基因组）时，可以用reference流程对数据进行详细的分析。Reference后面所有的流程都是基于参考序列进行的，所以选择正确的参考序列十分重要。1.2信息分析流程 得到测序序列后，即可利用比对软件，将所测序列比对到参考基因或基因组上，并进行后续分析，信息分析流程图如下：图二：转录组信息流程1.2.1原始fq序列简介测序得到的原始图像数据经base calling转化为序列数据，我们称

3、之为raw data或raw reads，结果以fastq文件格式存储，fastq文件为用户得到的最原始文件，里面存储reads的序列以及reads的测序质量。在fastq格式文件中每个read由四行描述：read IDTGGCGGAGGGATTTGAACCC+bbbbbbbbabbbbbbbbbbb每个序列共有4行，第1行和第3行是序列名称(有的fq文件为了节省存储空间会省略第三行“”后面的序列名称)，由测序仪产生；第2行是序列；第4行是序列的测序质量，每个字符对应第2行每个碱基，第四行每个字符对应的ASCII值减去64，即为该碱基的测序质量值，比如h 对应的ASCII值为104，那么其对应

4、的碱基质量值是40。碱基质量值范围为0到40。表1为Solexa测序错误率与测序质量值简明对应关系，具体计算公式如下：Qphred =-10 log10(e)表 1 Solexa测序错误率与测序质量值简明对应关系测序错误率测序质量值对应字符5%13M1%20T0.1%300.01%40h1.2.2原始fq序列处理 某些原始序列带有adaptor 序列，或含有少量低质量序列。我们首先经过一系列数据处理以去除杂质数据，得到Clean reads。按如下步骤进行处理：1.去除含adaptor的reads2.去除N的比例大于10%的reads3.去除低质量reads（质量值Q = 5的碱基数占整个re

5、ad的50以上）4.获得 Clean reads原始序列数据经过去除杂质后得到的数据称为Clean reads，后续分析都基于Clean reads1.2.3比对使用短reads比对软件SOAP2/SOAPalignerLi, 2009 #155将clean reads分别比对到参考基因组和参考基因序列（允许两个碱基错配）。通过这一步骤，我们可以将测序得到的reads对应到基因及基因组上，后续分析都是基于上述比对结果。1.2.4基本生物信息分析结果基本信息分析结果包含以下内容：1 测序数据产量及与 Reference 比对结果概述 统计数据量的大小，得到测序数据产量；对soap结果进行处理得到

6、测序数据与Reference序列比对的概况。2 评价测序随机性 在转录组实验过程中，首先要通过物理或化学方法将转录本打断成短片段，然后上机测序。如果打断随机性差，reads偏向于来自基因特定区域，将会直接影响转录组的各项分析结果。利用reads在基因上的分布来评价打断随机性。由于不同参考基因有不同长度，我们把reads在基因上的位置标准化到相对位置（reads在基因上的位置与基因长度的比值），然后统计基因的不同位置比对上的reads数。如果打断随机性好，reads在基因各部位应分布得比较均匀。3 基因覆盖度、测序深度的分布基因测序覆盖度指每个基因被reads覆盖的百分比，其值等于基因中uniq

7、ue mapping reads覆盖的碱基数跟基因编码区所有碱基数的比值。测序深度指基因被reads覆盖的次数，其值等于reads覆盖到基因的碱基数与基因编码区所有碱基数的比值。4 Reads 在参考基因组上的分布该分析主要是以图形方式概括给出Reads在基因组各个位置的分布情况，以及该位置基因的分布情况。1.2.5高级生物信息分析结果高级生物信息分析包含以下结果：1 对基因结构进行优化通过比较测序结果和现有基因注释结果，对基因的 5端或3端进行延长。如图三所示，首先，将reads比对到基因组，提取基因组中被unique mapping reads覆盖的次数大于或等于某阈值（默认为2）且位置连

8、续的区域作为转录活性区(Transcription Active Region, TAR，图中蓝色方块区域)；然后通过paired-end reads（图中紫色线条）将不同的TAR连接形成潜在的gene model；最后，通过比较潜在gene model与现有基因注释的差别，对基因的5端和3端进行延长（图中表现的仅是基因3端发生延长的情况）。图三：基因结构优化2 鉴定基因的可变剪切可变剪切使一个基因产生多个mRNA转录本，不同mRNA可能翻译成不同蛋白。因此，通过可变剪切一个基因可能产生多个蛋白，极大地增加了蛋白多样性Black, 2003 #6Stamm, 2005 #21;Lareau,

9、2004 #22。虽然已知可变剪切在真核生物中普遍存在，但我们可能仍低估了可变剪切的比例，最近，基于高通量测序的可变剪切研究在人Pan, 2008 #3 Wang, 2008 #4 Sultan, 2008 #5、小鼠Tang, 2009 #18;Mortazavi, 2008 #19、拟南芥Filichkin, #156中发现了很多新的可变剪切事件。在生物体内，主要存在7种可变剪切类型：A）Exon skipping; B）Intron retention; C) Alternative 5 splice site; D) Alternative 3 splice site; E) Alte

10、rnative first exon; F) Alternative last exon; G) Mutually exclusive exon. 下图是我们利用高通量测序数据鉴别出来的7种可变剪切。图中每个位置的ExP.Level等于log2(Reads数)。图四：可变剪切示意图 A) Exon Skipping. 基因AK070385发生可变剪切形成两种不同的转录本，第1种转录本比第2种转录组本多一个外显子(exon), 我们将这种外显子称为inclusive exon, inclusive exon两侧的两个外显子称为 constitutive exon。B) Intron retent

11、ion. 基因AK072590发生可变剪切形成两种不同的转录本，第2种转录本由retained Intron 与两侧的外显子一起形成新的外显子。C) Alternative 5 splice site. 基因AK067602发生可变剪切形成两种不同的转录本，它们的3端剪切位点一致但5 端剪切位点不同。D) Alternative 3 splice site. 基因AK067602发生可变剪切形成两种不同的转录本，它们的5端剪切位点一致但3 端剪切位点不同。E) Alternative First Exon. 基因AK068497发生可变剪切形成两种不同的转录本，它们的不同之处在于第一个外显子不

12、同。F) Alternative Last Exon. 基因AK064908发生可变剪切形成两种不同的转录本，它们的不同之处在于最后一个外显子不同。G) Mutually Exclusive Exon. 基因AK101575发生可变剪切形成两种不同的转录本，两转录本之间相同的外显子称为constitutive exon，不同的外显子称为inclusive exon，两个inclusive exon不能同时存在与同一转录本中，只能分别存在于不同转录本中。 下面，概述检测可变剪切的算法。首先，我们使用软件“tophat” Trapnell, 2009 #1鉴定转录本的剪切位点 (junction

13、site)（使用软件默认参数），剪切位点给出了转录本不同外显子的边界及组合关系，如图五，我们检测到三个剪切位点，分别表明Exon1和Exon2连接在一起，Exon2和Exon3连接在一起，Exon1和Exon3连接在一起。图 五 剪切位点示意图然后，通过分析同一基因的所有剪切位点，找出各种可变剪切事件。分析算法如下：A) Exon Skipping.图 六 Exon Skipping算法示意图转录本1和转录本2分别同时检测到如图六所示三个剪切位点，可认为转录本1的Exon1、Exon2和Exon3存在Exon Skipping剪切方式；转录本2的Exon1、Exon3和Exon4也存在Exon

14、 Skipping剪切方式。B) Intron Retention图 七 Intron Retention算法示意图如图七所示，1）检测到Junction 1的存在，表明在某个成熟mRNA中Exon1和Exon2之间的Intron被剪切下来；2）Exon1和Exon2之间的Intron有90以上的区域均有unique mapping reads覆盖，说明在某个成熟mRNA中该intron被保留下来了（考虑到转录的exon通常也不是100被reads覆盖到，所以在这里以90为阈值）。若同时满足以上两个条件，则认为该基因Exon1和Exon2之间存在Intron Retention的可变剪切方式。

15、C) Alternative 5 Splice Site图 八 Alternative 5 Splice Site算法示意图如图八，一个转录本的Junction 1位点被检测到，并且Junction 2 和Junction 3 中有一个被检测到（它们共同点是3剪切位点和Junction 1相同，但5剪切位点和Junction 1不同），那么就认为Exon1和Exon2 存在Alternative 5 Splice Site的剪切方式。D) Alternative 3 Splice Site图 九 Alternative 3 Splice Site算法示意图如图九，一个转录本的Junction

16、1位点被检测到，并且Junction 2 和Junction 3 中有一个被检测到（它们共同点是5剪切位点和junction 1相同，但3剪切位点和junction 1不同），那么就认为Exon1和Exon2 存在Alternative 3 Splice Site的剪切方式。E) Alternative First Exon图 十 Alternative First Exon算法示意图如图十，首先，要求检测到如图所示两个junction位点；其次，不能检测到支持Exon1和Exon2与5端的Exons有连接的junction位点。要求以上两个条件同时满足，且这种情况出现在转录本的最5端，但不要

17、求Exon1为这个转录本的第一个外显子，也不要求被junction连接的外显子都是相邻的，如转录本2中的Exon2和Exon4。所以，图中转录本1的Exon1、Exon2和Exon3存在Alternative First Exon的可变剪切方式，转录本2中Exon1、Exon2和Exon4也存在Alternative First Exon的可变剪切方式。F) Alternative Last Exon图 十一 Alternative Last Exon算法示意图如图十一，转录本1为例，首先，要求检测到如图所示两个junction位点(Junction1和Junction2)；其次，不能检测到支

18、持Exon1和Exon2与3端的Exons有连接的junction位点。要求以上两个条件同时满足，且这种情况出现在转录本的最3端，但不要求Exon3为这个转录本的最后一个外显子，也不要求被junction连接的外显子都是相邻的，如转录本2中的Exon1和Exon4。所以，图中转录本1的Exon1、Exon2和Exon3存在Alternative Last Exon的可变剪切方式，转录本2中Exon1、Exon3和Exon4也存在Alternative Last Exon的可变剪切方式。G) Mutually Exclusive Exon图 十二 Mutually Exclusive Exon算

19、法示意图检测到如图十二所示四个junction位点，且不能检测到支持Exon2和Exon3有连接位点的junction位点，则认为该转录本的Exon1、Exon2、Exon3和Exon4之间存在Mutually Exclusive Exon的可变剪切方式。3 发现新转录本现有数据库中对转录本的注释可能还不全面，通过高通量测序我们能检测到新的转录本Mortazavi, 2008 #103。我们首先从潜在gene model中挑选出长度大于150bp且平均覆盖度大于2的gene model，再从中找出位于基因间区域（一个基因3端下游200bp到下一个基因5端上游200bp之间的区域）的潜在gene

20、 model作为候选的新转录本。4 基因结构以及 Reads 在基因组上分布的精确图形该分析主要是以图形方式概括给出Reads在基因组各个位置的分布情况，以及该位置基因的分布情况。我们画出Reads在最长的25条染色体上的分布图，该图为 SVG 矢量图，如果你的浏览器不支持 SVG，请安装 SVGView 插件。5 基因差异表达分析5.1基因表达量基因表达量的计算使用RPKM法（Reads Per Kb per Million reads）Mortazavi, 2008 #103，其计算公式为：设RPKM(A)为基因A的表达量，则C为唯一比对到基因A的reads数，N为唯一比对到基因组的总re

21、ads数，L为基因A编码区的碱基数。RPKM法能消除基因长度和测序量差异对计算基因表达的影响，计算得到的基因表达量可直接用于比较不同样品间的基因表达差异。如果一个基因存在多个转录本，则用该基因的最长转录本计算其测序覆盖度和表达量。 5.2差异分析 差异表达分析找出在不同样本间存在差异表达的基因，并对差异表达基因做GO功能分析和KEGG Pathway分析。 参照Audic S.等人发表在Genome Research上的基于测序的差异基因检测方法Audic, 1997 #8（该文献已被引用超过五百次），我们开发了严格的算法筛选两样本间的差异表达基因。假设观测到基因A对应的reads数为x，已知

22、在一个大文库中，每个基因的表达量只占所有基因表达量的一小部分，在这种情况下，p(x)的分布服从泊松分布：已知，样本一中唯一比对到基因组的总reads数为N1，样本二中唯一比对到基因组的总reads数为N2，样本一中唯一比对到基因A的总reads数为x，样本二中唯一比对到基因A的总reads数为y，则基因A在两样本中表达量相等的概率可由以下公式计算：然后，我们对差异检验的p value作多重假设检验校正，通过控制FDR（False Discovery Rate）来决定p value的域值。假设挑选了R个差异表达基因，其中S个是真正有差异表达的基因，另外V个是其实没有差异表达的基因，为假阳性结果。

23、希望错误比例QV/R平均而言不能超过某个可以容忍的值，比如1，则在统计时预先设定FDR不能超过0.01（Benjamini, Yekutieli. 2001）。在得到差异检验的FDR值同时，我们根据基因的表达量（RPKM值）计算该基因在不同样本间的差异表达倍数。FDR值越小，差异倍数越大，则表明表达差异越显著。在我们的分析中，差异表达基因定义为FDR0.001且倍数差异在2倍以上的基因。得到差异表达基因之后，我们对差异表达基因做GO功能分析和KEGG Pathway分析。GO功能分析一方面给出差异表达基因的GO功能分类注释；另一方面给出差异表达基因的GO功能显著性富集分析。GO功能分类注释给出

24、具有某个GO功能的基因列表及基因数目统计。GO功能显著性富集分析给出与基因组背景相比，在差异表达基因中显著富集的GO功能条目，从而给出差异表达基因与哪些生物学功能显著相关。该分析首先把所有差异表达基因向Gene Ontology数据库（http:/www.geneontology.org/）的各个term映射，计算每个term的基因数目，然后应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著富集的GO条目，其计算公式为其中，N为所有基因中具有GO注释的基因数目；n为N中差异表达基因的数目；M为所有基因中注释为某特定GO term的基因数目；m为注释为某特定GO term的差异表达基

25、因数目。计算得到的p value通过Bonferroni校正之后，以corrected p value0.05为阈值，满足此条件的GO term定义为在差异表达基因中显著富集的GO term。通过GO功能显著性富集分析能确定差异表达基因行使的主要生物学功能。我们的GO功能分析同时整合了表达模式聚类分析，研究人员能方便地看到具有某一功能的所有差异基因的表达模式。例，immune response为在差异表达基因中最显著富集的一个GO term（表2）。图十三显示了参与immune response的差异基因的表达模式。表2 在差异表达基因中显著富集的GO-termlog2Ratio图 十三 参与

26、immune response的差异基因表达模式聚类图KEGG Pathway分析在生物体内，不同基因相互协调行使其生物学功能，基于Pathway的分析有助于更进一步了解基因的生物学功能。KEGG是有关Pathway的主要公共数据库Kanehisa, 2008 #96，Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的Pathway。该分析的计算公式同GO功能显著性富集分析，在这里N为所有基因中具有Pathway注释的基因数目；n为N中差异表达基因的数目；M为所有基因中注释为某特定Pathway的基因数目；m为注

27、释为某特定Pathway的差异表达基因数目。FDR0.05的Pathway定义为在差异表达基因中显著富集的Pathway。通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。结果如表3所示。表3pathway显著性富集分析列表各列意义如下：#序号Pathway通路名DEGs with pathway annotation (2085)注释到该通路的差异表达基因的数目All genes with pathway annotation (8986)注释到该通路的所有基因的数目Pvalue超几何检验的 P 值QvalueQ 值（Q 0.05 为在差异表达基因中显著富

28、集的 Pathway）Pathway IDKEGG 数据库中的 Pathway ID注：Qvalue 0.05 的 pathway 在差异表达基因中显著富集，见表中红框所示。差异表达基因的 pathway 显著性富集分析不但得到最有意义的 pathway 列表，点击其中的 pathway 链接还将得到KEGG 数据库中 pathway 的详细信息，如点击表3第一列第三行的 B cell receptor signaling pathway，可以看到如图十四所示的详细信息，上调基因所在位置用红色标记，下调基因所在位置用绿色标记。图 十四 KEGG 数据库中 B cell receptor sig

29、naling pathway 的详细信息二、Reference 工作流程工作流程如下：2.1前期工作1）创建项目目录：由于每个子项目都有自己的子项目代码，且名字简洁，建议使用子项目代码为项目创建目录，随着手头做过项目的增加，如果有需要，建议先以时期为依据创建大目录，再在其下创建项目目录；2) 项目记录：随着项目的增加，所需记得项目各方面信息内容也会增加，如果需要的话，建议使用excel电子表格记录平时的项目信息，以方便查询，包括：项目名称、子项目代码、项目结果路径、开始时间、阶段性进展、结束时间、截止时间、网址链接等等；2.2写工作文件1）文件模板根据信息任务描述，选好两个文件的模板，放于所创

30、建的项目目录下；2）找fq文件 方法1：(根据文库名查找)find /share/fqdata10/solexa/ -name *ARAcqfTARAAPE*fq查找结果：/share/fqdata10/solexa/HSZ09076_ARAcqfT_transcriptome_Transcriptome/ARAcqfTARAAPE/100114_I649_0002_FC42T26AAXX/100114_I649_FC42T26AAXX_L7_ARAcqfTARAAPE/100114_I649_FC42T26AAXX_L7_ARAcqfTARAAPE_1.fq/share/fqdata10/s

31、olexa/HSZ09076_ARAcqfT_transcriptome_Transcriptome/ARAcqfTARAAPE/100114_I649_0002_FC42T26AAXX/100114_I649_FC42T26AAXX_L7_ARAcqfTARAAPE/100114_I649_FC42T26AAXX_L7_ARAcqfTARAAPE_2.fq方法2：(根据项目编号查找)cd /share/fqdata10/solexa/cd HSZ09076 敲入tab键查找结果：dr-xr-xr-x 3 solexa solexa 41 Jan 25 13:28 ARAcqfTARAAPEdr-xr-xr-x 3 solexa solexa 41 Jan 25 13:28 ARAcqfTBRAAPE方法3：(根据子项目代码查找)cd /share/fqdata10/solexa/cd *_ARAcqfT_*查找结果：dr-xr-xr-x 3 solexa solexa 41 Jan 25 13:28 ARAcqfTARAAPEdr-xr-xr-x 3 solexa solexa 41 Jan 25 13:28 ARAcqfTBRAAPE数据存放路径：一般在以下几