1 引言
RNA-seq 组装的目标就是以测序短序列为基础来重构全长转录本。由于二代测序技术的局限性,在一个测序单元内只能测到比较短的片段。尽管出现了很有前途的三代测序技术,比如太平洋生物技术公司的PacBio,可以通过单分子技术一条片段测得几kb,但是在当前日常的转录组测序中还没有使用(译者注:最近三代测序在全长转录组测序中正火热开展)。所以,在实践中,为了获得全长转录本序列,必须用相互重叠的短片段来组装出来。在理论上,存在两种方案。一种是有参考基因组的转录组组装,前提是有可用的参考基因组信息。RNA-seq的 reads首先被定位到基因组上,组装工作表现为解决定位在基因组上的reads和哪一个转录本匹配。另外一种途径是不依赖其他信息进行从头(de novo)组装。在无参考基因组信息的情况下,组装是基于RNA-seq reads相互间的序列相似性开展的。这两种方案都可以转变成数学计算问题,包括在位图中找到一系列的路径。由于存在天文数字多的可能路径,即使对小的转录组进行组装,也有一定问题。简单而言,通过枚举法去找到最佳值是不可能的,因此要采用各种试探策略和近似策略来开展组装。
转录组组装与基因组组装有很大不同。在基因组组装过程中,read覆盖度通常更加均一化(除了文库制备和测序技术引起的偏好)。而且测序深度的偏差在基因组组装中意味着存在重复序列。与之不同的是,在RNA-seq数据中,基因与基因之间,还有同一基因的不同转录本之间,其表达丰度的差异可能存在高达几个数量级的差异。尽管基因间高度差异化的表达丰度这个特性可以被用于检测并构建不同的转录本,但是也带来了挑战。它要求采用更高的测序深度来检测低表达丰度的基因和稀有事件。为了平衡基因间表达丰度的差异性,需要在实验进程中开展文库的均一化。本文不对这些研究方法展开详细描述,但是有助于牢记组装质量是同计算方法及数据的取得方式相关的。因为测序技术只是把RNA-seq文库内容转变成计算机数字化形式,所以获得高质量数据的一个关键要素就是文库制备。
在开展组装前必须对数据继续质量控制处理。在本章中,我们各选择了两个软件包来进行以定位为基础的有参考组装和从头组装。它们全都是公开发布的非商业软件。同其它计算方法一样,我们要意识到组装的结果是依赖于数据与计算方法整合的。典型地,每一个方法中的参数都是可以调整,因此,依赖这些方法和参数,即便使用同一数据,也会产生很不一样的组装结果。
本章最初简介组装中的问题和相应的解决方法。,我们将逐一介绍被选出的四种软件包,并采用同一套数据集来分别实践演示。这一套数据集来自ENCODE项目,包含一个个体的双端测序reads。为了减少数据量,只有定位到人类第18号染色体上的reads被用于测试。双端测序的reads是从 (
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeCaltechRna-Seq/)这个测序文件中提取的。提取出的这套数据是只含有344000对reads的小数据集,所以运行组装软件所用的时间不会太长。
2 方法
RNA-seq组装起源于二十世纪九十年代初的表达序列标签(EST)测序[1]。ESTs技术就需要聚类与组装[2]。聚类就是把相似的EST序列通过计算两两间的重叠度来分组成簇。簇的成员范围是可以自定义的,比如,可以指定两条序列之间存在大于40 bp的重叠区,且在这个重叠区序列一致性大于95%的为一簇。把reads聚类后,再在每一个聚类簇内独立地开展组装。尽管详细的过程方法已经改变,转录组装依旧包含了主要的两步:(1)找出属于同一基因组位点的reads;(2)将每一个位点上的reads构建出一个位图来展现转录本。ESTs同目前的RNA-seq数据有一个很大的不同之处,就在于典型的EST只是转录本的部分片段信息,而现今的高通量数据能展现全长转录本。尽管现在的二代测序平台的单个reads不能获得全长转录本,但是大量高覆盖度的数据使得构建全长转录本成为可能。
2.1 转录组装是不同于基因组组装
在早期EST时期,基因组组装程序也可以用于转录组。尽管技术上是可行的,但实际上转录组不再采用基因组组装程序。现在的基因组组装与转录组组装在基础层面就有很大不同。除了测序深度均一性不同,主要区别在于基因组组装中理想的输出结果是每一条线性序列代表基因组的一段区域,而在转录组中在同一基因组位点可以有多个异构体,即同一个基因的一个外显子可能在不同的转录本中同其它外显子产生多种组合。因此,在转录组组装中,基因自然表述为一个位图,并用节点表示外显子,弧线表示选择性剪切事件。连接节点的分枝对应于选择性剪切。一个独立的转录本是一个单分子,它应该被一条线性序列所表述,它在位图中由一条包含节点的路径来表示。一个外显子节点在一个异构体中只能存在一次,但是这个外显子可以在不同的异构体中存在多次。在一个外显子的位图上所有可能的路径包括了所有可能的异构体。这些可能的路径数量巨大,但是只有少数是在转录组中真实存在的。转录组装中一个挑战性的工作就是将所有可能的候选真实的异构体中找出来。重申一下,这个问题的根源在于测序产生的是短片段序列。如果我们可以把一个转录本的序列一次测完整,这个问题将不再存在。我们采用短序列片段来构建长序列时这种组合问题才会暴露出来。
2.2 重构转录本的复杂性
为了展现重构转录本的复杂性,我们用下面一个例子来说明。我们假定存在这样一个基因,它有三个外显子,那么由一个外显子,两个外显子和三个外显子分别构成的全部可能异构体数量分别是3,3和1个,总数量为7个。针对含有N个外显子的基因,其全部可能出现的转录本总数量的通用计算公式为:
这就是说,在一个异构体中对每个外显子而言都存在两种可能,存在与不存在。
包括了没有外显子即空异构体的情况,所以需要在总数量中减去1(对应于公式中的 k=0 )。这只是全部可能的异构体数量。在转录组中,不论多少种异构体存在都是可能的。尽管选择性剪切提供了多种组合的可能,但是在真实的转录组中并不是所有的都存在的。异构体的存在与否所构成的不同转录本数据集,沿用上面的公式,可以用下面的公式计算出来数据集的可能的总套数:
可以看出,单外显子数量N增长的时候,总套数增长的非常迅猛。当N=1,2,3,4或者5的时候,所有可能的转录本总套数分别为1,7,127,32767和2147483647。这也表明当外显子数量超过4个时,通过枚举法来把所有的可能都检测变得不具备实践性。
2.3 组装进程
目前有两种策略来进行转录本的重构,即图位法和从头组装。它们都需要利用RNA-seq数据来对每一个位点进行位图的构建。位图是充当一个起点来解析异构体。两种策略都面临着同一个问题,就是如何把数据分离从而使得一个位图标识一个位点。
本书第四章已经介绍了定位法。任何容许分离reads的方法都可以被用于将RNA-seq的reads匹配到基因组上。如果基因结构信息已知,就是说知道哪个外显子属于哪个基因。如果基因结构信息未知,定位上的reads必须要被分割以便标识基因位点。一个外显子位图,也叫做一个剪切位图,首先被用于构建每个位点,之后的任务就是在每个位图中查找一系列路径,每个路径就代表了一条异构体。
通过限制外显子位图中的链接,可以减少异构体的数目。每一个链接代表一个外显子结合点。在一个完整的链接位图中个,由于所有的节点都有弧线相连,代表了所有的异构体都是可能的存在的。计算任务就是选择最同数据相吻合的一个拓扑学位图。那些在RNA-seq测序序列中得不到支持的选择性剪切事件都会被丢掉,只有位图中的必要的链接会被保留。保留一条弧线所需要的支持证据包括分割的read序列和双端信息。对于一个分割的read而言,如果该read的前端序列被定位到一个外显子而它的末端序列被定位到另外一个外显子,这一现象就支持我们认为这两个外显子是在同一个转录本序列中的邻近位置。对于双端测序而言,一对reads的两端,一端定位到一个外显子而另一端定位到另外一个外显子,也支持上述结论。对一个外显子剪切点而言,存在分割的read比一对双端reads更具有支持力。因为对双端定位的一对reads而言,为了避免出现这两个外显子虽然在同一个转录本内,但是它们两个外显子之间还可能存在其它外显子等情况,必须要利用测序文库插入片段大小信息来确认这两个外显子是真的结合在一起的。插入片段大小分布是直接依赖于测序文库构建过程的。通常,如果文库内的插入片段大小的方差比较大,对每一对reads就使用平均插入片段大小,因为对每一对都无法准确的预测其插入片段大小。
从头组装有两个基本的途径:(1)计算reads间成对的重叠区从而给出组装位图的拓扑结构,(2)构建de Bruijn位图,即用一套k-mers及它们间的连接性来描绘整个测序数据。作为一种数学模型,de Bruijn位图在测序技术出现以前就被引入到计算中[3],而在基因组组装中由Pevzner等人第一次使用[4]。从头组装的目标是从代表了基因组或转录组的一部分原始信息的组装位图中提取到尽可能长的连续片段(contigs,即重叠群)。在1990年代的人类基因组计划时,测序长度相对较长(他们采用的是Sanger测序),当然测序的总量比现在的二代技术要少。当时基因组组装是以测序片段间的重叠为基础延伸的,这种策略也叫做重叠-排列-一致(OLC),也代表了组装过程的三个阶段。尽管计算所有reads间的两两重叠区是很费时,但是从方法上讲这是组装问题中最容易的部分。主要的困难在于组合学:如何确定多个组合在一起的序列在位图的上排列(译者注:即重复序列的定位问题)。虽然可以构建一种最佳算法用来解决组装问题,但是对任何一套实际数据而言它的运行时间都将非常长,所以必须采用多种启发式和近似值算法(5)。当测序产生的数据量增加,同时测序片段变短,这种情况下采用de Bruijn位图策略就变得更流行。在转录组组装中,现在大部分的方法都应用了de Bruijn位图发。然而,还是有一些例外,比如MIRA EST组装器[6]就是基于OLC范式的。
2.4de Bruijn图
de Bruijn位图de Bruijn位图上的每个节点是同(k-1)-mer序列相关联。如果存在一个k-mer,它的前端序列是节点A的(k-1)-mer,末端序列是节点B的(k-1)-mer,那么节点A和B就连接在一起。用这种方法,k-mer在de Bruijn位图中制造边线[7]。在位图中,路径代表序列,即使是一条测序读长也被分成多个相连的节点,第一个节点就是从序列的起始位置开始的包含(k-1)-mer长,第二个节点就是从序列的第二个碱基开始的包含(k-1)-mer长,等等类似。每一个k-mer在位图中只被使用一次,用来做两个节点间的边线。如果两条序列有一个相同的k-mer,则它们两条序列可以共享边线。这就不需要明确地计算reads间两两比对结果,就可以获得序列间的重叠关系信息。相比在所有reads间计算重叠关系,构建de Bruijn位图是更加直白和快速的。它包括从reads中简单提取所有k-mer并将代表(k-1)-mer的节点连接起来。随之而来的挑战就是如何找到位图中的合适路径来代表真实的转录本。测序错误造成的尖端(tips)在位图中就是死胡同,泡状结构(bubbles)将使位图结构复杂化。泡状结构的形成是由于位图中的分枝同该位图中的另外一部分反向重合。一些泡状结构是由于测序错误产生的,一些事由于选择性剪切造成的,例如,一个外显子在基因中部,但是它可能存在于某一异构体中却在另外一个异构体中发生了外显子跳跃。这就使得在一个位图中有两条路径,它们有共同的前端与末端,但是在中部产生了分枝。在位图中寻找路径时要利用到双端测序信息和k-mer信息。通过k-mer丰度我们可以给位图中的边线计算权重,从而减少错误路径的出现。K-mer长度会影响位图的复杂度。直白地说,k-mer必须比read长度短,但是如果太短,由于节点特异性降低会造成生成的位图变得非常密。相反,如果k-mer较大,必须要有足够的测序数据才能构建好位图。为了解决这个问题,选出一个合适的k值,可以尝试用不同的k值来进行多次组装,从而选出最优的一次组装结果或者把不同k-mer下组装出来的contigs(重叠群)组合在一起成一套结果[8-10]。
2.5 丰度信息的应用
如果一套候选异构体具有合理的大小,那就可以用RNA-seq的丰度信息来解析这些异构体。原因就是在同一个转录本中所有的外显子都应该具有相同的表达丰度。一个转录本是一个分子,因此如果在文库制备与测序(及定位)过程中没有失误,测序reads应该均匀地分布和覆盖整条转录本。如果存在差异,比如一些外显子具有更深的测序深度,这意味着这些外显子也存在于其它异构体中。
针对一套修正过的异构体,预测它们的相对丰度是可能的。最佳的方法是获得最有说明力的丰度值。这个具体步骤第一步是给丰度设定初始值,例如,将丰度在所有异构体间平均地分割,然后运用最大期望算法(EM)迭代计算。EM算法是1970年代被开发的[11],在转录组数据中用这个算法详细描述在文献[12]中。这个最佳的迭代过程有两步:计算期望(E)和最大化(M)。在技术期望这一步骤中,根据异构体的丰度,所有的reads被成比例地分配到每个异构体上。在最大化这一步骤中,异构体的相对丰度被重新计算。这两步被重复进行直到估算的丰度值不再改变,这也就是说,算法完成了聚集。这个方法是针对给定的一套异构体,因此若有新的异构体被加入这套测试集,所有的值都可能改变。通常而言,由于EM算法找出一个局部最优解,在存在多个局部最优解的情况下,最终方案需要依赖起始值。然而,只能有一个最大化,所有局部最优解也就是全局最优解[13]。基础的EM算法可以在多方面被修订。例如,在iReckon软件[14]中,EM算法被正则化,以便降低假性转录本重构的数量。
广告:
转录组分析视频教程,欢迎登陆网易云课堂观看:
扫码获取
更多精彩内容,欢迎关注生信人
TCGA | 小工具 | 数据库 |组装| 注释 | 基因家族 | Pvalue
基因预测 |bestorf | sci | NAR | 在线工具 | 生存分析 | 热图
生信不死 | 初学者 | circRNA | 一箭画心| 十二生肖 | circos
舞台|基因组 | 黄金测序 | 套路 | 杂谈组装 | 进化 | 测序简史