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等方法来反复迭代的考虑转录本组装与比到达量估计。

药企,独角兽,苏州。团队长期招人,感兴趣的都可以发邮件聊聊:tiehan@sina.cn
个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn