【T】每日一生信--RNA测序数据回帖和组装
(2013-11-27 20:59:50)
[<cite style="font-style: normal; line-height: 23px; padding: 0px 20px 0px 32px; display: inline-block !important; height: 23px; background-repeat: no-repeat; background-position: right -28px; min-width: 1px; overflow-x: visible; white-space: nowrap; text-align: left; font-size: 12px !important; background-image: url("../../../image/8/8_50/images/sg_newsp.png"); width: 34px;">
转载▼</cite>](javascript:;)
| 标签：
readsmapping
joinexon
splitreads
cufflinks
| 分类： RNAseq分析 |
已整理到新地址
<wbr> <wbr> Reads mapping往往被作为测序深度分析的第一步，其质量的好坏与速度的快慢，都会直接对后续的分析工作产生影响。同样是基于深度测序技术，在reads长度，数量，质量等方面，RNA-Seq产生的reads与之前基因组重测产生的DNA reads具有相似的性质；长度短，数量大，质量参差不齐，错误率较高。然而，RNA-Seq测序数据来源于RNA转录本，在DNA转录成mRNA的过程中，内含子被切掉，外显子会在剪切位点连接到一起。对于这些跨过剪切位点的reads，也就是所谓的junction reads,如果不把他们从中间断开，就无法准确贴到基因组上，这些junction reads是确定剪切位点的直接证据，对于正确重构转录本结构至关重要。
<wbr> <wbr> 如上图，途中横跨exion <wbr>1和exon 3的junction reads,据可以直接支持存在外显子1和3直接相连，中间不包括exon 2的转录本的存在。我们的算法在mapping时，需要意识到junction site和intron的存在，从而正确处理这些junction reads.具体来说，目前解决这个问题主要有两类策略。
<wbr> <wbr> 其一是join exon,这个策略的第一步是基于已知转录本中所有的exon来构建所有可能的junctions.这个library中的junction未必是已知的，而是所有可能的组合。例如对于4个exons,就会有6种组合。而后，在mapping时，首先直接采取之前DNA reads类似的unspliced方式基因组比对，而非junction reads map到基因组。而对于直接无法map的junction reads，则再与第一步构建的junction library进行比对。join exon策略使之上相当于对之前DNA reads mapping算法的一个“补丁”（patch）.通过构建所有肯能的Junction库，join exon策略可能发现新的splicing isoform。然而，对于以前未知的exon,这个策略就无能为力了。
<wbr> <wbr> 为了克服这个问题，可以使用split reads策略。与之前类似，slit reads策略在mapping时，也是首先采取之前DNA reads类似的nspliced方式将非junction reads map到基因组。对于无法直接map的junction reads,则参照之前BLAST的方法，切成若干长度为K的种子（seed）,再利用这些seed重试mapping,也就是在更小的粒度上尝试寻找junction site,最后将临近的mapped seeds重新组合起来，得到最终的全read的alignment.与之前的join exon策略相比,split reads策略因为要在更短的seeds mapping而影响了速度。但这个策略不依赖于先验（prior）的exon注释，可以发现新的exon乃至新的基因。
<wbr> <wbr> <wbr>事实上，目前常用的RNA-Seq工具通常会组合两个策略，以平衡检测的灵敏度与速度。例如TopHat 2工具中，首先利用join exon策略检测已知的junction site,而后再利用split reads策略来发现新的junctions,值得注意的是，TopHat 2针对不同的阶段采用了不同的索引（index）,可以进一步提高mapping的速度。
<wbr> <wbr> 完成mapping只是RNA-Seq数据分析的第一步，之后，我们需要将这些reads组装成转录本，并对每个转录本估计相应的表达量。在包括junction reads在内的所有reads均正确map之后，我们就可以将转录本组织问题描述为一个对有向图（directed graph）的遍历问题。通过对图中的边给予不同的权值（weight）作为约束（constrain）,我们利用图论（graph theory）中的寻路算法（pateh finding algorithm）找到一条或多条符合约束的最优路径以及其对应的转录本序列。
<wbr> <wbr>接下来以常用的Cufflinks为例，简要的说一下相关的想法。
<wbr> <wbr>Cufflinks是一个准对RNA-Seq数据进行转录本组装和表达分析的工具软件。
<wbr> <wbr>首先，我们假设，我们不知道这里面有这样的三个转录本结构，我们只看到了这些reads,首先Cufflinks会去找那些不可能出现在同一转录本中的片段（fragment）,比如说这里黄色和蓝色的片段，他们就不可能出现在同一个转录本里，如果他们来自同一个转录本，黄色就应该在蓝色相应的位置断开，而不是直接跨过去。同理，红色，黄色，蓝色的片段都是亮亮不相容的，而同一颜色的片段则是彼此相容的。通过将每个相容片段作为节点，并与它最近且相容的片段相连，就可以得到重叠图（overlap graph）.基于精简原则（parsimony principle）,Cufflinks在途中选择能覆盖所有reads额路径中互补相连且最少的一组路径，作为最优路径。来据此得到最终的三个转录本的集合。原则上，在转录本组织正确完成，且每个exon水的表达量均利用我们在上一单元中提到的归一（normalization）正确处理后，转录本的表达量可以直接从各个exon的表达量推出。
<wbr> <wbr> 例如，我们假定从基因组上的三个exon可以转录出两个转录本t1和t2,同时，每个exon表达水平在归一化后可以确定为e1=20,e2=40,e3=60,那么根据转录本的结构，我们就可以直接得到转录本表达水平与exon表达水平的关系。例如，exon1只在转录本1中出现，那么它的所有的表达就完全由转录本1来贡献，类似的，exon3在转录本1和2中同时出现，因此这两个转录本都对其表达有贡献。所以，我们认为，外显子1的表达量就对于转录本1的表达量。而外显子3的表达量应该对应为转录本1和转录本2表达量的和。以此类推，我们就可以推出转录本1和2各自的表达量：20，40
当然，考虑到对reads的分配本质上由转录本拼接算法决定的，在实际中这个问题会更加复杂。比如说，在Cufflinks中，对于reads的分配还会同时考虑长度分布等因素。事实上，为了更为准确的估计表达量，通常会采用EM等方法来反复迭代的考虑转录本组装与比到达量估计。
