转录组ref流程工作手册一、Reference 流程生物学原理1.1 实验流程图一：转录组实验流程当我们得到样品时，必须对其测序，才能得到分析所需的数据。
测序基本过程：提取样品总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 IDTGGCGGAGGGATTTGAACCC+bbbbbbbbabbbbbbbbbbb每个序列共有4行，第1行和第3行是序列名称(有的fq文件为了节省存储空间会省略第三行“＋”后面的序列名称)，由测序仪产生；第2行是序列；第4行是序列的测序质量，每个字符对应第2行每个碱基，第四行每个字符对应的ASCII值减去64，即为该碱基的测序质量值，比如h 对应的ASCII值为104，那么其对应的碱基质量值是40。
碱基质量值范围为0到40。
表 1为Solexa测序错误率与测序质量值简明对应关系，具体计算公式如下：Q phred =-10 log10(e)表 1 Solexa测序错误率与测序质量值简明对应关系5% 13 M1% 20 T0.1% 30 ^0.01% 40 h1.2.2原始fq序列处理某些原始序列带有adaptor 序列，或含有少量低质量序列。
我们首先经过一系列数据处理以去除杂质数据，得到Clean reads。
按如下步骤进行处理：1.去除含adaptor的reads2.去除N的比例大于10%的reads3.去除低质量reads（质量值Q <= 5的碱基数占整个read的50％以上）4.获得 Clean reads原始序列数据经过去除杂质后得到的数据称为Clean reads，后续分析都基于Clean reads1.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的剪切方式。