转录组ref流程工作手册 - 图文

转录组ref流程工作手册

一、Reference 流程生物学原理 1.1 实验流程

Total RNAeukaryoteEnriched mRNA by OligoTprokaryoteRemove rRNARNA fragment(200~700 bp)Random hexamer primed cDNA synthesisSize selection,then PCR amplificationSolexa Sequencing 图一：转录组实验流程

当我们得到样品时，必须对其测序，才能得到分析所需的数据。测序基本过程：提取样品总RNA后，用带有Oligo(dT)的磁珠富集真核生物mRNA（若为原核生物，则用试剂盒去除rRNA后进入下一步）。加入fragmentation buffer将mRNA打断成短片段，以mRNA为模板，用六碱基随机引物（random hexamers）合成第一条cDNA链，然后加入缓冲液、dNTPs、RNase H 和DNA polymerase I合成第二条cDNA链，在经过QiaQuick PCR试剂盒纯化并加EB缓冲液洗脱之后做末端修复并连接测序接头，然后用琼脂糖凝胶电泳进行片段大小选择，最后进行PCR扩增，使用建好的测序文库进行测序。

得到RNA的序列后，又可以找到它的参考序列（物种本身的基因、基因组）

时，可以用reference流程对数据进行详细的分析。Reference后面所有的流程都是基于参考序列进行的，所以选择正确的参考序列十分重要。

1.2信息分析流程

得到测序序列后，即可利用比对软件，将所测序列比对到参考基因或基因组上，并进行后续分析，信息分析流程图如下：

图二：转录组信息流程

1.2.1原始fq序列简介

测序得到的原始图像数据经base calling转化为序列数据，我们称之为raw data或raw reads，结果以fastq文件格式存储，fastq文件为用户得到的最原始文件，里面存储reads的序列以及reads的测序质量。在fastq格式文件中每个read由四行描述： @read ID

TGGCGGAGGGATTTGAACCC

+

bbbbbbbbabbbbbbbbbbb

每个序列共有4行，第1行和第3行是序列名称(有的fq文件为了节省存储空间会省略第三行“＋”后面的序列名称)，由测序仪产生；第2行是序列；第4行是序列的测序质量，每个字符对应第2行每个碱基，第四行每个字符对应的ASCII值减去64，即为该碱基的测序质量值，比如h 对应的ASCII值为104，那么其对应的碱基质量值是40。碱基质量值范围为0到40。表 1为Solexa测序错误率与测序质量值简明对应关系，具体计算公式如下：

Qphred =-10 log10(e)

表 1 Solexa测序错误率与测序质量值简明对应关系

测序错误率 5% 1% 0.1% 0.01% 测序质量值 13 20 30 40 对应字符 M T ^ h

1.2.2原始fq序列处理

某些原始序列带有adaptor 序列，或含有少量低质量序列。我们首先经过一系列数据处理以去除杂质数据，得到Clean reads。

按如下步骤进行处理：

1. 去除含adaptor的reads

2. 去除N的比例大于10%的reads

3. 去除低质量reads（质量值Q <= 5的碱基数占整个read的50％以上） 4. 获得 Clean reads

原始序列数据经过去除杂质后得到的数据称为Clean reads，后续分析都基于Clean reads

1.2.3比对

使用短reads比对软件SOAP2/SOAPaligner{Li, 2009 #155}将clean reads分别比对到参考基因组和参考基因序列（允许两个碱基错配）。

通过这一步骤，我们可以将测序得到的reads对应到基因及基因组上，后续分析都是基于上述比对结果。

1.2.4基本生物信息分析结果

基本信息分析结果包含以下内容：

1 测序数据产量及与 Reference 比对结果概述

统计数据量的大小，得到测序数据产量；对soap结果进行处理得到测序数据与Reference序列比对的概况。 2 评价测序随机性

在转录组实验过程中，首先要通过物理或化学方法将转录本打断成短片段，然后

上机测序。如果打断随机性差，reads偏向于来自基因特定区域，将会直接影响转录组的各项分析结果。

利用reads在基因上的分布来评价打断随机性。由于不同参考基因有不同长度，我们把reads在基因上的位置标准化到相对位置（reads在基因上的位置与基因长度的比值），然后统计基因的不同位置比对上的reads数。如果打断随机性好，reads在基因各部位应分布得比较均匀。

3 基因覆盖度、测序深度的分布

基因测序覆盖度指每个基因被reads覆盖的百分比，其值等于基因中unique mapping reads覆盖的碱基数跟基因编码区所有碱基数的比值。测序深度指基因被reads覆盖的次数，其值等于reads覆盖到基因的碱基数与基因编码区所有碱基数的比值。 4 Reads 在参考基因组上的分布

该分析主要是以图形方式概括给出Reads在基因组各个位置的分布情况，以及该位置基因的分布情况。

1.2.5高级生物信息分析结果

高级生物信息分析包含以下结果：

1 对基因结构进行优化

通过比较测序结果和现有基因注释结果，对基因的 5'端或3'端进行延长。如图三所示，首先，将reads比对到基因组，提取基因组中被unique mapping reads覆盖的次数大于或等于某阈值（默认为2）且位置连续的区域作为转录活性区(Transcription Active Region, TAR，图中蓝色方块区域)；然后通过paired-end reads（图中紫色线条）将不同的TAR连接形成潜在的gene model；最后，通过比较潜在gene model与现有基因注释的差别，对基因的5'端和3'端进行延长（图中表现的仅是基因3’端发生延长的情况）。

图三：基因结构优化

2 鉴定基因的可变剪切

可变剪切使一个基因产生多个mRNA转录本，不同mRNA可能翻译成不同蛋白。因此，通过可变剪切一个基因可能产生多个蛋白，极大地增加了蛋白多样性{Black, 2003 #6}{Stamm, 2005 #21;Lareau, 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) Alternative 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 retention. 基因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发生可变剪切形成两种不同的转录本，它们的不同之处在于第一个外显子不同。F) Alternative Last Exon. 基因AK064908发生可变剪切形成两种不同的转录本，它们的不同

之处在于最后一个外显子不同。G) Mutually Exclusive Exon. 基因AK101575发生可变剪切形成两种不同的转录本，两转录本之间相同的外显子称为constitutive exon，不同的外显子称为inclusive exon，两个inclusive exon不能同时存在与同一转录本中，只能分别存在于不同转录本中。

下面，概述检测可变剪切的算法。首先，我们使用软件“tophat” {Trapnell, 2009 #1}鉴定转录本的剪切位点 (junction 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 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的可变剪切方式。

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 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’端，但不要求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)；其次，不能检测到支持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算法示意图

检测到如图十二所示四个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 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为唯一比对到基因组的总reads数，L为基因A编码区的碱基数。RPKM法能消除基因长度和测序量差异对计算基因表达的影响，计算得到的基因表达量可直接用于比较不同样品间的基因表达差异。

如果一个基因存在多个转录本，则用该基因的最长转录本计算其测序覆盖度和表达量。

5.2差异分析

差异表达分析找出在不同样本间存在差异表达的基因，并对差异表达基因做GO功能

分析和KEGG Pathway分析。

参照Audic S.等人发表在Genome Research上的基于测序的差异基因检测方法{Audic, 1997 #8}（该文献已被引用超过五百次），我们开发了严格的算法筛选两样本间的差异表达基因。

假设观测到基因A对应的reads数为x，已知在一个大文库中，每个基因的表达量只占所有基因表达量的一小部分，在这种情况下，p(x)的分布服从泊松分布：

已知，样本一中唯一比对到基因组的总reads数为N1，样本二中唯一比对到基因组的总reads数为N2，样本一中唯一比对到基因A的总reads数为x，样本二中唯一比对到基因A的总reads数为y，则基因A在两样本中表达量相等的概率可由以下公式计算：

然后，我们对差异检验的p value作多重假设检验校正，通过控制FDR（False Discovery Rate）来决定p value的域值。假设挑选了R个差异表达基因，其中S个是真正有差异表达的基因，另外V个是其实没有差异表达的基因，为假阳性结果。希望错误比例Q＝V/R平均而言不能超过某个可以容忍的值，比如1％，则在统计时预先设定FDR不能超过0.01（Benjamini, Yekutieli. 2001）。在得到差异检验的FDR值同时，我们根据基因的表达量（RPKM值）计算该基因在不同样本间的差异表达倍数。FDR值越小，差异倍数越大，则表明表达差异越显著。在我们的分析中，差异表达基因定义为FDR≤0.001且倍数差异在2倍以上的基因。 得到差异表达基因之后，我们对差异表达基因做GO功能分析和KEGG Pathway分析。 GO功能分析一方面给出差异表达基因的GO功能分类注释；另一方面给出差异表达基因的GO功能显著性富集分析。

GO功能分类注释给出具有某个GO功能的基因列表及基因数目统计。

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

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

表2 在差异表达基因中显著富集的GO-term

log2Ratio

图 十三 参与immune response的差异基因表达模式聚类图

KEGG Pathway分析 在生物体内，不同基因相互协调行使其生物学功能，基于Pathway的分析有助于更进一步了解基因的生物学功能。KEGG是有关Pathway的主要公共数据库{Kanehisa, 2008 #96}，Pathway显著性富集分析以KEGG Pathway为单位，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的Pathway。该分析的计算公式同GO功能显著性富集分析，在这里N为所有基因中具有Pathway注释的基因数目；n为N中差异表达基因的数目；M为所有基因中注释为某特定Pathway的基因数目；m为注释为某特定Pathway的差异表达基因数目。FDR≤0.05的Pathway定义为在差异表达基因中显著富集的Pathway。通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。结果如表3所示。

表3 pathway显著性富集分析列表

各列意义如下： # Pathway DEGs with pathway annotation (2085) 序号 通路名 注释到该通路的差异表达基因的数目 All genes with pathway annotation 注释到该通路的所有基因的数目 (8986) Pvalue Qvalue Pathway ID 超几何检验的 P 值 Q 值（Q ≤ 0.05 为在差异表达基因中显著富集的 Pathway） KEGG 数据库中的 Pathway ID 注：Qvalue ≤ 0.05 的 pathway 在差异表达基因中显著富集，见表中红框所示。

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

图 十四 KEGG 数据库中 B cell receptor signaling pathway 的详细信息

二、Reference 工作流程

工作流程如下：

2.1前期工作

1） 创建项目目录：

由于每个子项目都有自己的子项目代码，且名字简洁，建议使用子项目代码为项目创建目录，随着手头做过项目的增加，如果有需要，建议先以时期为依据创建大目录，再在其下创建项目目录； 2) 项目记录：

随着项目的增加，所需记得项目各方面信息内容也会增加，如果需要的话，建议使用excel电子表格记录平时的项目信息，以方便查询，包括：项目名称、子项目代码、项目结果路径、开始时间、阶段性进展、结束时间、截止时间、网址链接等等；

2.2写工作文件

1）文件模板

根据信息任务描述，选好两个文件的模板，放于所创建的项目目录下； 2）找fq文件

方法1：(根据文库名查找)

find /share/fqdata10/solexa/ -name \查找结果：

/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/solexa/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 ARAcqfTARAAPE dr-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 ARAcqfTARAAPE dr-xr-xr-x 3 solexa solexa 41 Jan 25 13:28 ARAcqfTBRAAPE

数据存放路径：

一般在以下几个库中：

/share/fqdata12/solexa/ （2010年2-3月的数据） /share/fqdata10/solexa/ （2010年1-2月的数据）

/share/fastdata1/solexa （09年11月份下机的数据）

/share/solid2/solexa-work/Project_solexa_fq （09年10-11月份下机的数据） /share/solid1/solexa-work/Project_solexa_fq （09年9-10月份下机的数据）

以下是09年9月之前可以查找的：

/share/raid007/ solexa-work/Project_solexa_fq /share/raid009/ solexa-work/Project_solexa_fq /share/raid7/ solexa-work/Project_solexa_fq

3）找参考序列（包括参考基因组、参考基因、psl文件）

如合作伙伴提供参考序列，则使用合作伙伴提供的参考序列。如合作伙伴未

提供，找到相关数据后，将链接发送给合作伙伴确认可行后方能使用。

本文来源：https://www.bwwdw.com/article/rgta.html