转录组ref流程工作手册.docx
《转录组ref流程工作手册.docx》由会员分享，可在线阅读，更多相关《转录组ref流程工作手册.docx（32页珍藏版）》请在冰豆网上搜索。
转录组ref流程工作手册
转录组ref流程工作手册
转录组ref流程工作手册
一、Reference流程生物学原理
1.1实验流程
图一：
转录组实验流程
当我们得到样品时，必须对其测序，才能得到分析所需的数据。
测序基本过程：
提取样品总RNA后，用带有Oligo（dT）的磁珠富集真核生物mRNA（若为原核生物，则用试剂盒去除rRNA后进入下一步）。
加入fragmentationbuffer将mRNA打断成短片段，以mRNA为模板，用六碱基随机引物（randomhexamers）合成第一条cDNA链，然后加入缓冲液、dNTPs、RNaseH和DNApolymeraseI合成第二条cDNA链，在经过QiaQuickPCR试剂盒纯化并加EB缓冲液洗脱之后做末端修复并连接测序接头，然后用琼脂糖凝胶电泳进行片段大小选择，最后进行PCR扩增，使用建好的测序文库进行测序。
得到RNA的序列后，又可以找到它的参考序列（物种本身的基因、基因组）时，可以用reference流程对数据进行详细的分析。
Reference后面所有的流程都是基于参考序列进行的，所以选择正确的参考序列十分重要。
1.2信息分析流程
得到测序序列后，即可利用比对软件，将所测序列比对到参考基因或基因组上，并进行后续分析，信息分析流程图如下：
图二：
转录组信息流程
1.2.1原始fq序列简介
测序得到的原始图像数据经basecalling转化为序列数据，我们称之为rawdata或rawreads，结果以fastq文件格式存储，fastq文件为用户得到的最原始文件，里面存储reads的序列以及reads的测序质量。
在fastq格式文件中每个read由四行描述：
@readID
TGGCGGAGGGATTTGAACCC
+
bbbbbbbbabbbbbbbbbbb
每个序列共有4行，第1行和第3行是序列名称（有的fq文件为了节省存储空间会省略第三行“＋”后面的序列名称），由测序仪产生；第2行是序列；第4行是序列的测序质量，每个字符对应第2行每个碱基，第四行每个字符对应的ASCII值减去64，即为该碱基的测序质量值，比如h对应的ASCII值为104，那么其对应的碱基质量值是40。
碱基质量值范围为0到40。
表1为Solexa测序错误率与测序质量值简明对应关系，具体计算公式如下：
Qphred=-10log10（e）
表1Solexa测序错误率与测序质量值简明对应关系
测序错误率
测序质量值
对应字符
5%
13
M
1%
20
T
0.1%
30
^
0.01%
40
h
1.2.2原始fq序列处理
某些原始序列带有adaptor序列，或含有少量低质量序列。
我们首先经过一系列数据处理以去除杂质数据，得到Cleanreads。
按如下步骤进行处理：
1.去除含adaptor的reads
2.去除N的比例大于10%的reads
3.去除低质量reads（质量值Q<=5的碱基数占整个read的50％以上）
4.获得Cleanreads
原始序列数据经过去除杂质后得到的数据称为Cleanreads，后续分析都基于Cleanreads
1.2.3比对
使用短reads比对软件SOAP2/SOAPaligner{Li,2009#155}将cleanreads分别比对到参考基因组和参考基因序列（允许两个碱基错配）。
通过这一步骤，我们可以将测序得到的reads对应到基因及基因组上，后续分析都是基于上述比对结果。
1.2.4基本生物信息分析结果
基本信息分析结果包含以下内容：
1测序数据产量及与Reference比对结果概述
统计数据量的大小，得到测序数据产量；对soap结果进行处理得到测序数据与Reference序列比对的概况。
2评价测序随机性
在转录组实验过程中，首先要通过物理或化学方法将转录本打断成短片段，然后上机测序。
如果打断随机性差，reads偏向于来自基因特定区域，将会直接影响转录组的各项分析结果。
利用reads在基因上的分布来评价打断随机性。
由于不同参考基因有不同长度，我们把reads在基因上的位置标准化到相对位置（reads在基因上的位置与基因长度的比值），然后统计基因的不同位置比对上的reads数。
如果打断随机性好，reads在基因各部位应分布得比较均匀。
3基因覆盖度、测序深度的分布
基因测序覆盖度指每个基因被reads覆盖的百分比，其值等于基因中uniquemappingreads覆盖的碱基数跟基因编码区所有碱基数的比值。
测序深度指基因被reads覆盖的次数，其值等于reads覆盖到基因的碱基数与基因编码区所有碱基数的比值。
4Reads在参考基因组上的分布
该分析主要是以图形方式概括给出Reads在基因组各个位置的分布情况，以及该位置基因的分布情况。
1.2.5高级生物信息分析结果
高级生物信息分析包含以下结果：
1对基因结构进行优化
通过比较测序结果和现有基因注释结果，对基因的5'端或3'端进行延长。
如图三所示，首先，将reads比对到基因组，提取基因组中被uniquemappingreads覆盖的次数大于或等于某阈值（默认为2）且位置连续的区域作为转录活性区（TranscriptionActiveRegion,TAR，图中蓝色方块区域）；然后通过paired-endreads（图中紫色线条）将不同的TAR连接形成潜在的genemodel；最后，通过比较潜在genemodel与现有基因注释的差别，对基因的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）Exonskipping;B）Intronretention;C）Alternative5’splicesite;D）Alternative3’splicesite;E）Alternativefirstexon;F）Alternativelastexon;G）Mutuallyexclusiveexon.下图是我们利用高通量测序数据鉴别出来的7种可变剪切。
图中每个位置的ExP.Level等于log2（Reads数）。
图四：
可变剪切示意图
A）ExonSkipping.基因AK070385发生可变剪切形成两种不同的转录本，第1种转录本比第2种转录组本多一个外显子（exon）,我们将这种外显子称为inclusiveexon,inclusiveexon两侧的两个外显子称为constitutiveexon。
B）Intronretention.基因AK072590发生可变剪切形成两种不同的转录本，第2种转录本由retainedIntron与两侧的外显子一起形成新的外显子。
C）Alternative5’splicesite.基因AK067602发生可变剪切形成两种不同的转录本，它们的3’端剪切位点一致但5’端剪切位点不同。
D）Alternative3’splicesite.基因AK067602发生可变剪切形成两种不同的转录本，它们的5’端剪切位点一致但3’端剪切位点不同。
E）AlternativeFirstExon.基因AK068497发生可变剪切形成两种不同的转录本，它们的不同之处在于第一个外显子不同。
F）AlternativeLastExon.基因AK064908发生可变剪切形成两种不同的转录本，它们的不同之处在于最后一个外显子不同。
G）MutuallyExclusiveExon.基因AK101575发生可变剪切形成两种不同的转录本，两转录本之间相同的外显子称为constitutiveexon，不同的外显子称为inclusiveexon，两个inclusiveexon不能同时存在与同一转录本中，只能分别存在于不同转录本中。
下面，概述检测可变剪切的算法。
首先，我们使用软件“tophat”{Trapnell,2009#1}鉴定转录本的剪切位点（junctionsite）（使用软件默认参数），剪切位点给出了转录本不同外显子的边界及组合关系，如图五，我们检测到三个剪切位点，分别表明Exon1和Exon2连接在一起，Exon2和Exon3连接在一起，Exon1和Exon3连接在一起。
图五剪切位点示意图
然后，通过分析同一基因的所有剪切位点，找出各种可变剪切事件。
分析算法如下：
A）ExonSkipping.
图六ExonSkipping算法示意图
转录本1和转录本2分别同时检测到如图六所示三个剪切位点，可认为转录本1的Exon1、Exon2和Exon3存在ExonSkipping剪切方式；转录本2的Exon1、Exon3和Exon4也存在ExonSkipping剪切方式。
B）IntronRetention
图七IntronRetention算法示意图
如图七所示，1）检测到Junction1的存在，表明在某个成熟mRNA中Exon1和Exon2之间的Intron被剪切下来；2）Exon1和Exon2之间的Intron有90％以上的区域均有uniquemappingreads覆盖，说明在某个成熟mRNA中该intron被保留下来了（考虑到转录的exon通常也不是100％被reads覆盖到，所以在这里以90％为阈值）。
若同时满足以上两个条件，则认为该基因Exon1和Exon2之间存在IntronRetention的可变剪切方式。
C）Alternative5’SpliceSite
图八Alternative5’SpliceSite算法示意图
如图八，一个转录本的Junction1位点被检测到，并且Junction2和Junction3中有一个被检测到（它们共同点是3’剪切位点和Junction1相同，但5’剪切位点和Junction1不同），那么就认为Exon1和Exon2存在Alternative5’SpliceSite的剪切方式。
D）Alternative3’SpliceSite
图九Alternative3’SpliceSite算法示意图
如图九，一个转录本的Junction1位点被检测到，并且Junction2和Junction3中有一个被检测到（它们共同点是5’剪切位点和junction1相同，但3’剪切位点和junction1不同），那么就认为Exon1和Exon2存在Alternative3’SpliceSite的剪切方式。
E）AlternativeFirstExon
图十AlternativeFirstExon算法示意图
如图十，首先，要求检测到如图所示两个junction位点；其次，不能检测到支持Exon1和Exon2与5’端的Exons有连接的junction位点。
要求以上两个条件同时满足，且这种情况出现在转录本的最5’端，但不要求Exon1为这个转录本的第一个外显子，也不要求被junction连接的外显子都是相邻的，如转录本2中的Exon2和Exon4。
所以，图中转录本1的Exon1、Exon2和Exon3存在AlternativeFirstExon的可变剪切方式，转录本2中Exon1、Exon2和Exon4也存在AlternativeFirstExon的可变剪切方式。
F）AlternativeLastExon
图十一AlternativeLastExon算法示意图
如图十一，转录本1为例，首先，要求检测到如图所示两个junction位点（Junction1和Junction2）；其次，不能检测到支持Exon1和Exon2与3’端的Exons有连接的junction位点。
要求以上两个条件同时满足，且这种情况出现在转录本的最3’端，但不要求Exon3为这个转录本的最后一个外显子，也不要求被junction连接的外显子都是相邻的，如转录本2中的Exon1和Exon4。
所以，图中转录本1的Exon1、Exon2和Exon3存在AlternativeLastExon的可变剪切方式，转录本2中Exon1、Exon3和Exon4也存在AlternativeLastExon的可变剪切方式。
G）MutuallyExclusiveExon
图十二MutuallyExclusiveExon算法示意图
检测到如图十二所示四个junction位点，且不能检测到支持Exon2和Exon3有连接位点的junction位点，则认为该转录本的Exon1、Exon2、Exon3和Exon4之间存在MutuallyExclusiveExon的可变剪切方式。
3发现新转录本
现有数据库中对转录本的注释可能还不全面，通过高通量测序我们能检测到新的转录本{Mortazavi,2008#103}。
我们首先从潜在genemodel中挑选出长度大于150bp且平均覆盖度大于2的genemodel，再从中找出位于基因间区域（一个基因3’端下游200bp到下一个基因5’端上游200bp之间的区域）的潜在genemodel作为候选的新转录本。
4基因结构以及Reads在基因组上分布的精确图形
该分析主要是以图形方式概括给出Reads在基因组各个位置的分布情况，以及该位置基因的分布情况。
我们画出Reads在最长的25条染色体上的分布图，该图为SVG矢量图，如果你的浏览器不支持SVG，请安装SVGView插件。
5基因差异表达分析
5.1基因表达量
基因表达量的计算使用RPKM法（ReadsPerKbperMillionreads）{Mortazavi,2008#103}，其计算公式为：
设RPKM（A）为基因A的表达量，则C为唯一比对到基因A的reads数，N为唯一比对到基因组的总reads数，L为基因A编码区的碱基数。
RPKM法能消除基因长度和测序量差异对计算基因表达的影响，计算得到的基因表达量可直接用于比较不同样品间的基因表达差异。
如果一个基因存在多个转录本，则用该基因的最长转录本计算其测序覆盖度和表达量。
5.2差异分析
差异表达分析找出在不同样本间存在差异表达的基因，并对差异表达基因做GO功能分析和KEGGPathway分析。
参照AudicS.等人发表在GenomeResearch上的基于测序的差异基因检测方法{Audic,1997#8}（该文献已被引用超过五百次），我们开发了严格的算法筛选两样本间的差异表达基因。
假设观测到基因A对应的reads数为x，已知在一个大文库中，每个基因的表达量只占所有基因表达量的一小部分，在这种情况下，p（x）的分布服从泊松分布：
已知，样本一中唯一比对到基因组的总reads数为N1，样本二中唯一比对到基因组的总reads数为N2，样本一中唯一比对到基因A的总reads数为x，样本二中唯一比对到基因A的总reads数为y，则基因A在两样本中表达量相等的概率可由以下公式计算：
然后，我们对差异检验的pvalue作多重假设检验校正，通过控制FDR（FalseDiscoveryRate）来决定pvalue的域值。
假设挑选了R个差异表达基因，其中S个是真正有差异表达的基因，另外V个是其实没有差异表达的基因，为假阳性结果。
希望错误比例Q＝V/R平均而言不能超过某个可以容忍的值，比如1％，则在统计时预先设定FDR不能超过0.01（Benjamini,Yekutieli.2001）。
在得到差异检验的FDR值同时，我们根据基因的表达量（RPKM值）计算该基因在不同样本间的差异表达倍数。
FDR值越小，差异倍数越大，则表明表达差异越显著。
在我们的分析中，差异表达基因定义为FDR≤0.001且倍数差异在2倍以上的基因。
得到差异表达基因之后，我们对差异表达基因做GO功能分析和KEGGPathway分析。
GO功能分析一方面给出差异表达基因的GO功能分类注释；另一方面给出差异表达基因的GO功能显著性富集分析。
GO功能分类注释给出具有某个GO功能的基因列表及基因数目统计。
GO功能显著性富集分析给出与基因组背景相比，在差异表达基因中显著富集的GO功能条目，从而给出差异表达基因与哪些生物学功能显著相关。
该分析首先把所有差异表达基因向GeneOntology数据库（http:
//www.geneontology.org/）的各个term映射，计算每个term的基因数目，然后应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著富集的GO条目，其计算公式为
其中，N为所有基因中具有GO注释的基因数目；n为N中差异表达基因的数目；M为所有基因中注释为某特定GOterm的基因数目；m为注释为某特定GOterm的差异表达基因数目。
计算得到的pvalue通过Bonferroni校正之后，以correctedpvalue≤0.05为阈值，满足此条件的GOterm定义为在差异表达基因中显著富集的GOterm。
通过GO功能显著性富集分析能确定差异表达基因行使的主要生物学功能。
我们的GO功能分析同时整合了表达模式聚类分析，研究人员能方便地看到具有某一功能的所有差异基因的表达模式。
例，immuneresponse为在差异表达基因中最显著富集的一个GOterm（表2）。
图十三显示了参与immuneresponse的差异基因的表达模式。
表2 在差异表达基因中显著富集的GO-term
log2Ratio
图十三参与immuneresponse的差异基因表达模式聚类图
KEGGPathway分析
在生物体内，不同基因相互协调行使其生物学功能，基于Pathway的分析有助于更进一步了解基因的生物学功能。
KEGG是有关Pathway的主要公共数据库{Kanehisa,2008#96}，Pathway显著性富集分析以KEGGPathway为单位，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的Pathway。
该分析的计算公式同GO功能显著性富集分析，在这里N为所有基因中具有Pathway注释的基因数目；n为N中差异表达基因的数目；M为所有基因中注释为某特定Pathway的基因数目；m为注释为某特定Pathway的差异表达基因数目。
FDR≤0.05的Pathway定义为在差异表达基因中显著富集的Pathway。
通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。
结果如表3所示。
表3 pathway显著性富集分析列表
各列意义如下：
#
序号
Pathway
通路名
DEGswithpathwayannotation（2085）
注释到该通路的差异表达基因的数目
Allgeneswithpathwayannotation（8986）
注释到该通路的所有基因的数目
Pvalue
超几何检验的P值
Qvalue
Q值（Q≤0.05为在差异表达基因中显著富集的Pathway）
PathwayID
KEGG数据库中的PathwayID
注：
Qvalue≤0.05的pathway在差异表达基因中显著富集，见表中红框所示。
差异表达基因的pathway显著性富集分析不但得到最有意义的pathway列表，点击其中的pathway链接还将得到KEGG数据库中pathway的详细信息，如点击表3第一列第三行的Bcellreceptorsignalingpathway，可以看到如图十四所示的详细信息，上调基因所在位置用红色标记，下调基因所在位置用绿色标记。
图十四KEGG数据库中Bcellreceptorsignalingpathway的详细信息
二、Reference工作流程
工作流程如下：
2.1前期工作
1）创建项目目录：
由于每个子项目都有自己的子项目代码，且名字简洁，建议使用子项目代码为项目创建目录，随着手头做过项目的增加，如果有需要，建议先以时期为依据创建大目录，再在其下创建项目目录；
2）项目记录：
随着项目的增加，所需记得项目各方面信息内容也会增加，如果需要的话，建议使用excel电子表格记录平时的项目信息，以方便查询，包括：
项目名称、子项目代码、项目结果路径、开始时间、阶段性进展、结束时间、截止时间、网址链接等等；
2.2写工作文件
1）文件模板
根据信息任务描述，选好两个文件的模板，放于所创建的项目目录下；
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/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/
cdHSZ09076敲入tab键
查找结果：
dr-xr-xr-x3solexasolexa41Jan2513:
28ARAcqfTARAAPE
dr-xr-xr-x3solexasolexa41Jan2513:
28ARAcqfTBRAAPE
方法3：
（根据子项目代码查找）
cd/share/fqdata10/solexa/
cd*_ARAcqfT_*
查找结果：
dr-xr-xr-x3solexasolexa41Jan2513:
28ARAcqfTARAAPE
dr-xr-xr-x3solexasolexa41Jan2513:
28ARAcqfTBRAA
升级会员 