扫码联系客服
烧脑的RNA-seq数据分析
RNA-seq数据分析的最开始两步是：
1.Reads与基因组或转录组进行序列比对；
2.对各个转录本进行定量分析。
目前应用最广泛的的TopHat2算法是2013年发表的，用于将RNA-seq得到的Reads比对到基因组或者转录组上。TopHat2比对过程可分为以下3部分：
将Reads比对到已知转录本（可选）；
能比对到单个外显子的Reads比对到基因组的外显子；
跨越内含子的Reads切割成没有overlap的片段比对回基因组。
由此可见，TopHat2是尽量将每一个Read的碱基序列作为一个整体进行序列比对的，只在Read跨越内含子时才进行切割，这样就保留了完整Read中的有用信息，但同时也会花费大量的计算时间。例如，用TopHat2比对20个样品，每个含有30M Reads，启用20个核同时工作，需要花费28小时。由于TopHat2只能进行序列比对，接下来的定量分析还需要用配套的Cufflinks软件再运行14小时才算完结。序列比对的耗时耗力以及对计算机硬件的高要求无疑是一件让人头疼的事，尤其是在如今测序数据越来越多的情况下。而且，RNA-seq的测序结果并不是分析一次就结束了，很多研究者会对已发表的测序结果进行汇总和进一步分析，以期得到新的发现，这就要求对序列比对和定量分析的算法进行优化。
脑洞大开的算法优化
一种有效降低计算工作量的做法是把Read分割成k-mers，其实就是将一个相对较长的Read化整为零，分割为有k个连续碱基的k-mer。简单的说，一个长度为n的Read可以分成n-k+1个k-mers，前一个k-mer和后一个k-mer有k-1个碱基的overlap。例如，一个长度为10的Read序列为AGATCGAGTG（当然实际比这个长），将其分为3-mers，就得到了10-3+1，共8个3-mers。
Read: AGATCGAGTG
3-mers: AGA GAT ATC TCG CGA GAG AGT GTG
分割成k-mer后，将每一个k-mer分配到hash表中一个唯一的位置，再进行序列比对。通过这种转换，可以大大提高序列比对的效率。然而这种方法也面临着一个问题，那就是有些k-mers可能会比对到基因组的不同位置上，这就降低了定量分析的准确度。新算法kallisto则有效地解决了这个问题，它应用了一种思想：我们并不需要知道Reads来源于转录本的哪个具体位置，只要知道是哪个转录本就可以精确定量了。
kallisto怎么看？
kallisto也是首先将Reads简化为k-mers，然后使用hash表构造de Bruijn图。虽然hash表的存储需要额外的内存，但根据hash表的特点，各个k-mer只需存储一次，因此消耗内存很少。将de Bruijn图结构与计算机存储方式结合到一起，实现了相互之间的转化。kallisto还应用了一些其它的优化来进一步提高效率，最后用最大期望算法进行定量分析。
下图是kallisto的流程图，咦——似乎并不怎么看得懂。。。我只能告诉你a图中黑色的是Read，彩色的是不同的转录本，b图是de Bruijn图构建过程，剩下的是基于hash表的计算过程，然后就自行脑补吧。
kallisto棒棒的！
算法介绍完了，我们来看看效果如何。不怕不识货，就怕货比货，与其它主流的算法相比，kallisto在准确度上水平相当，还比其中两种算法（Cufflinks和Sailfish）提高了不少（如下图a所示）。但在所用的时间上，kallisto大大优于其它的算法。（b图最右边的那条细线就是。算法测试比对了20个样品，每个含有30M Reads，启用20个核同时工作）。
当然，以上测试是作者自己进行的，难免会有王婆卖瓜的嫌疑，因此，小编还找到另外一个课题组发表在Genome Biology上的一篇文章，专门用7种方法来测试包括kallisto在内的7种算法在序列处理准确性上的优劣（不考虑耗时），结论是包括kallisto在内的5种算法不相上下，而另外2种算法较差。所以我们可以得出结论：kallisto在不损害准确度的情况下，计算速率大大提升！
结尾有彩蛋：
今年2月份发表在Nature Biotechnology上的另外一篇文章介绍了一种算法（Sequence Bloom Tree，SBT），用于在NIH的海量RNA-seq数据库中寻找你感兴趣的短片段序列，该算法比当前的算法快162倍！相信这些算法的应用能够大大节约你的宝贵时间，这样一来，貌似不能按时毕业就只能怪自己了。。。
[1] Nicolas L Bray, et al. "Near-optimal probabilistic RNA-seq quantification."Nature Biotechnology (2016).
[2] Mingxiang Teng, et al. "A benchmark for RNA-seq quantification pipelines." Genome Biology (2016).
[3] Solomon, Brad, and C. Kingsford. "Fast search of thousands of short-read sequencing experiments." Nature Biotechnology (2016).
