RNA测序数据回帖和组装

Reads mapping往往被作为测序深度分析的第一步，其质量的好坏与速度的快慢，都会直接对后续的分析工作产生影响。同样是基于深度测序技术，在reads长度，数量，质量等方面，RNA-Seq产生的reads与之前基因组重测产生的DNA reads具有相似的性质；长度短，数量大，质量参差不齐，错误率较高。然而，RNA-Seq测序数据来源于RNA转录本，在DNA转录成mRNA的过程中，内含子被切掉，外显子会在剪切位点连接到一起。对于这些跨过剪切位点的reads，也就是所谓的junction reads,如果不把他们从中间断开，就无法准确贴到基因组上，这些junction reads是确定剪切位点的直接证据，对于正确重构转录本结构至关重要。

如上图，途中横跨exion 1和exon 3的junction reads,据可以直接支持存在外显子1和3直接相连，中间不包括exon 2的转录本的存在。我们的算法在mapping时，需要意识到junction site和intron的存在，从而正确处理这些junction reads.具体来说，目前解决这个问题主要有两类策略。

其一是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,这个策略就无能为力了。

为了克服这个问题，可以使用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乃至新的基因。

事实上，目前常用的RNA-Seq工具通常会组合两个策略，以平衡检测的灵敏度与速度。例如TopHat 2工具中，首先利用join exon策略检测已知的junction site,而后再利用split reads策略来发现新的junctions,值得注意的是，TopHat 2针对不同的阶段采用了不同的索引（index）,可以进一步提高mapping的速度。

完成mapping只是RNA-Seq数据分析的第一步，之后，我们需要将这些reads组装成转录本，并对每个转录本估计相应的表达量。在包括junction reads在内的所有reads均正确map之后，我们就可以将转录本组织问题描述为一个对有向图（directed graph）的遍历问题。通过对图中的边给予不同的权值（weight）作为约束（constrain）,我们利用图论（graph theory）中的寻路算法（pateh finding algorithm）找到一条或多条符合约束的最优路径以及其对应的转录本序列。

接下来以常用的Cufflinks为例，简要的说一下相关的想法。

Cufflinks是一个准对RNA-Seq数据进行转录本组装和表达分析的工具软件。

首先，我们假设，我们不知道这里面有这样的三个转录本结构，我们只看到了这些reads,首先Cufflinks会去找那些不可能出现在同一转录本中的片段（fragment）,比如说这里黄色和蓝色的片段，他们就不可能出现在同一个转录本里，如果他们来自同一个转录本，黄色就应该在蓝色相应的位置断开，而不是直接跨过去。同理，红色，黄色，蓝色的片段都是亮亮不相容的，而同一颜色的片段则是彼此相容的。通过将每个相容片段作为节点，并与它最近且相容的片段相连，就可以得到重叠图（overlap graph）.基于精简原则（parsimony principle）,Cufflinks在途中选择能覆盖所有reads额路径中互补相连且最少的一组路径，作为最优路径。来据此得到最终的三个转录本的集合。原则上，在转录本组织正确完成，且每个exon水的表达量均利用我们在上一单元中提到的归一（normalization）正确处理后，转录本的表达量可以直接从各个exon的表达量推出。

例如，我们假定从基因组上的三个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 每日一生信–RNA测序数据回帖和组装

当然，考虑到对reads的分配本质上由转录本拼接算法决定的，在实际中这个问题会更加复杂。比如说，在Cufflinks中，对于reads的分配还会同时考虑长度分布等因素。事实上，为了更为准确的估计表达量，通常会采用EM等方法来反复迭代的考虑转录本组装与比到达量估计。