[0071] 在以下方程式(3)-(5)中概述最大值函数B[j,k]的变量参数,其中MISMATCH_PENALTY、MATCH_BONUS、INSERTION_PENALTY、DELETION_PENALTY和OPENING_PENALTY都是常量,并且除MATCH_BONUS以外都是负数。通过以下方程式(3)给出匹配变量参数p[j,k]:
[0072] 如果S[j]≠A[k],则p[j,k]=max(p[[j-1,k-1],i[j-1,k-1],d[j-1,k-1])+MISMATCH_PENALTY(3)
[0073] 如果S[j]=A[k],则p[j,k]=max(p[[j-1,k-],i[j-1,k-1],d[j-1,k-1])+MATCH_BONUS,
[0074] 通过以下方程式(4)给出插入变量参数i[j,k]:
[0075] i[j,k]=max(p[j-l,k]+OPENING_PENALTY,i[j-l,k],d[j-l,k]+(4)
[0076] OPENING_PENALTY)+INSERTION_PENALTY
[0077] 并且缺失变量参数d[j,k]通过以下方程式(5)给出:
[0078] d[j,k]=max(p[j,k-1]+OPENING_PENALTY,i[j,k-1]+(5)
[0079] OPENING_PENALTY,d[j,k-1])+DELETION_PENALTY
[0080] 对于所有三个变量参数,将[0,0]元素设置为零以确保回溯完成,即,p[0,0]=i[0,0]=d[0,0]=0。
[0081] 评分参数是稍微任意的,并可经调整以实现计算的行为。DNA的评分参数设置的一个示例(Huang,《第3章:生物序列比较和比对(Bio-Sequence Comparison and Alignment)》,当前顶级比较分子生物学(Curr Top Comp Mol Biol.)丛书,马萨诸塞州剑桥市(Cambridge,Mass.):麻省理工学院出版社(The MIT Press),2002)将是:
[0082] MATCH_BONUS:10
[0083] MISMATCH_PENALTY:-20
[0084] INSERTION_PENALTY:-40
[0085] OPENING_PENALTY:-10
[0086] DELETION_PENALTY:-5
[0087] 以上空位罚分(INSERTION_PENALTY、OPENING_PENALTY)之间的关系有助于限制空位开放的数目,即,支持通过设置高于空位开放成本的空位插入罚分来将空位集合在一起。当然,MISMATCH_PENALTY、MATCH_BONUS、INSERTION_PENALTY、OPENING_PENALTY和DELETION_PENALTY之间可能存在可替代的关系。
[0088] 在一些实施例中,本发明的方法及系统并入多维比对算法。本发明的多维算法提供了序列信息的“向后看”类型分析(如在史密斯-沃特曼中),其中通过包含多个通路和多个节点的多维空间进行向后看。多维算法可以被用于将RNA测序读数与DAG类型参考进行比对。该比对算法通过关于包含在DAG(例如,参考序列构造)上的位置处的每个序列识别最大评分为Ci,j识别最大值。事实上,通过在先前位置处“向后”看,有可能横越多个可能的路径识别最佳比对。
[0089] 在上面所论述的读数(亦称“字符串”)和有向非循环图表(DAG)上执行本发明的算法。出于定义该算法的目的,假设S是要比对的字符串,并且假设D是S将进行比对的有向非循环图表。以从1开始的索引对字符串S的元素加括号。因此,如果S是字符串ATCGAA,那么S[1]=A、S[4]=G等。
[0090] 在某些实施例中,对于DAG,节点的序列的每个字母将被表示为独立元素d。d的前趋被定义为:
[0091] (i)如果d不是其节点的序列的首字母,则其节点中在d之前的字母是其(唯一)前趋;
[0092] (ii)如果d是其节点的序列的首字母,则为d的节点的
父节点的任何节点(例如,基因组中上游的所有外显子)的序列的最后一个字母是d的前趋。
[0093] 所有前趋的组继而表示为P[d]。
[0094] 为了找出“最优”比对,算法寻求M[j,d](S的前j个元素与在d之前(且包含d)的DAG的一部分的最佳比对的评分)的值。该步骤类似于在以上的方程式1中发现Hi,j。具体来说,确定M[j,d]包括找出a、i、e和0的最大值,如下面所定义:
[0095] M[j,d]=max{a,i,e,0}(6)
[0096] 其中
[0097] 对于在P[d]中的p*,e=max{M[j,p*]+DELETE_PENALTY}
[0098] i=M[j-1,d]+INSERT_PENALTY
[0099] 对于在P[d]中的p*,如果S[j]=d,则a=max{M[j-1,p*]+MATCH_SCORE};
[0100] 对于在P[d]中的p*,如果S[j]≠d,则max{M[j-1,p*]+MISMATCH_PENALTY}[0101] 如上面所描述,e是S的前j个字符与DAG直到(但不包含)d的部分的比对加上另外的DELETE_PENALTY的最高值。因此,如果d不是节点的序列的首字母,则仅存在一个前趋p,并且S的前j个字符与DAG(直到且包含p)的比对评分等效于M[j,p]+DELETE_PENALTY。在其中d是其节点的序列的首字母的实例中,可以存在多个可能的前趋,并且因为DELETE_PENALTY是常量,所以求[M[j,p*]+DELETE_PENALTY]的最大值和选择具有与S的前j个字符比对的最高比对评分的前趋相同。
[0102] 在方程式(6)中,i是字符串S的前j-1个字符与DAG的直到并且包含d的部分的比对加上INSERT_PENALTY,其类似于SW中的插入变量参数的定义(参看方程式1)。
[0103] 另外,a是S的前j个字符与DAG的直到但不包含d的部分比对的最高值,加上或者MATCH_SCORE(如果S的第j个字符与字符d相同)或者MISMATCH_PENALTY(如果S的第j个字符与字符d不同)。如同e一样,这意味着如果d不是其节点的序列的首字母,则仅存在一个前趋,即p。这意味着a是S的前j-1个字符与DAG(直到并且包含p)的比对评分,即,M[j-1,p](再加上或者MISMATCH_PENALTY或者MATCH_SCORE)取决于d与S的第j个字符是否匹配。在其中d是其节点的序列的首字母的实例中,可以存在多个可能的前趋。在该情况下,求{M[j,p*]+MISMATCH_PENALTY或MATCH_SCORE}的最大值和以下相同:选择具有与S的前j-1个字符的最高比对评分(即,候选M[j-1,p*]变量参数的最高值)的前趋并且取决于d与S的第j个字符是否匹配而加上或者MISMATCH_PENALTY或者MATCH_SCORE。
[0104] 再次,与在SW算法中相同,罚分例如DELETE_PENALTY、INSERT_PENALTY、MATCH_SCORE和MISMATCH_PENALTY可以被调整以促进与较少空位等的比对。
[0105] 如以上方程式中所描述,该算法通过不仅计算该元素的插入、缺失和匹配评分,而且向后看(逆着DAG的方向)到DAG上的任何之前节点以找出最大评分,来找出每个读数的最佳(即,最大)值。因此,该算法能够详细研究穿过DAG的含有已知突变的不同路径。因为图表是有向的,所以逆着图表方向移动的回溯跟随朝向图表起点的优选异构体,并且最佳比对评分识别高度确定性的最可能比对。
[0106] 图1示出与DAG比较的读数的图形表示。图1的顶部区域呈现序列读数“ATCGAA”以及参考序列TTGGATATGGG(序列ID号:1)和已知的插入事件TTGGATCGAATTATGGG(序列ID号:2),其中插入是带下划线的。图1的中间示出读数、序列ID号:l和序列ID号:2之间的关系,保持每个序列线性以辅助观察关系。图1的底部区域示出使用DAG构造的比对。在所描绘的DAG中,即使沿不同路径,但可以通过从DAG的5'端到DAG的3'端读取来读取序列ID号1和序列ID号2两者。如所描绘的,序列读数被示出为与上部路径比对。
[0107] 图2示出对应于比较的实际矩阵。相似于史密斯-沃特曼技术,本发明的所例示算法识别最高评分,并且实行回溯以识别读数的正确位置。图1和图2还突出本发明产生字符串与该构造的实际匹配。在其中序列读数包含未被包含在DAG中的变异体的情况下,将通过空位、插入等报告经比对的序列。
[0108] 将RNA测序读数与经注释的参考的DAG进行比对允许发现转录组中的异构体。另外,已经发生的不同形式的可变剪接的可能性(例如,一个密码子已在相关外显子中的一个中被跳过的可能性)可以被“免费”预期。
[0109] 外显子可自动地分成两个节点,一个对应于第一密码子,并且另一个对应于外显子的其余部分。这允许在DAG中表示两种类型的剪接,即使还未观测到它们。鉴于一个密码子的与比对的差值的价值将往往不足以惩罚候选比对以诱导我们寻找替代方式,预测这些情况可以是有用的。(最大的此罚分将是三个插入或三个缺失的罚分,并且由于机率,所以常常罚分将是较小的。)
[0110] 将转录组表示为DAG准许将异构体自然地解释为穿过DAG的路径。与往常一样,当系统和方法表示数学结构中的物理项时,益处是大量的。例如,本发明的系统和方法允许序列以直观的方式可视化。例如,可以通过示出穿过DAG的两个路径来比较序列。
[0111] 另外,本发明的方法有益于利用图表定理和算法。穿过DAG的路径的分析是在现代的数学和计算机科学中充满生机的研究区域。因此,利用该研究来改进转录组研究的能力是有价值的。例如,如果边缘被给定与由读数跨越的接合点的概率对应的权重,则路径的总权重可以携载关于实现不同异构体的先验可能性的信息。穿过DAG的路径最大化(或最小化)权重是已知的,并且最
短路径算法可以被用于快速地找出那些最短路径。
[0112] 图3提供了本发明的方法的图式。如上面所论述的,一旦DAG内置有以规范的顺序连接的外显子,则获得双端读数。然后,双端读数与DAG比对,并且在读数和包括表示RNA序列的节点和连接多对节点的边缘的DAG之间找出比对(具有满足预定指标的比对评分)。预定指标可以是例如最高评分比对。比对可以通过上述改良的史密斯-沃特曼算法进行。对于一对双端读数,比对可以包含基于穿过包含比对的有向非循环数据结构的路径,识别转录组中的转录本。通过使用插入大小作为比对的约束,可以有效地实行比对。双端读数可以被用于构建DAG(例如,添加边缘)。
[0113] 该方法适合于分析由RNA测序获得的读数。在某些实施例中,有向非循环数据结构基本上表示至少一个
染色体的所有已知的特征。可以在找出比对之前创建和存储有向非循环数据结构。在一些实施例中,通过将每个序列读数与穿过有向非循环数据结构的可能路径中的至少大部分进行比较来找出比对。方法可以包含基于找出的比对将多个序列读数装配成重叠群。
[0114] 3.使用插入大小作为比对的约束
[0115] 其中插入大小信息可以与DAG组合的一种方式是使用插入大小作为比对的约束。通过使用该方法,速度和精确性都有可能增加。因为该方法防止丢弃与比对相关的重要信息,所以与单独比对双端配偶相比,该方法可以更精确。应当注意的是,独立比对两个双端配偶近似地相当于在基因组中L个可能的位置(如果基因组具有长度L)之间选择两个位置。
因此,比对包括在大约L^2种可能性之间进行选择。然而,第一比对被用于强有力地将另一比对约束到在第一比对周围的常量大小k的一些区域,然后仅有大约kL种可能性——在L而不是L^2的量级上。这将更加快速。
[0116] 图4-图6例示使用插入大小信息以约束比对。
[0117] 图4例示任意结构的假设的DAG和定义插入大小k的一对假设的双端读数。在如所绘制的DAG中,节点是方
块,并且边缘是线。在所描绘的实施例中,从具有均值k bp的近似正态分布获取插入大小k。将该对的第一成员与DAG比对(例如,使用上述史密斯-沃特曼或改良的史密斯-沃特曼)。一般来说,双端包括片段,该片段的每端上具有衔接子,并且该片段具有从片段-衔接子边界朝向彼此朝内延伸的双端读数。在双端读数之间可存在空位。插入长度可以意味着片段的长度——即,从一个片段-衔接子边界到另一个片段-衔接子边界的距离。参见例如Lindgreen,2012,《AdapterRemoval:下一代序列读数的简单清理(AdapterRemoval:easy cleaning of next-generation sequence reads)》,《BMC研究纪要《(BMC Res Notes)》5:337(例如,Lindgreen,图1,画面C),其内容以引用方式并入。
[0118] 在图5中,最左边黑色方块表示外显子,该对的第一成员在外显子处比对。确定了一组距离,它们占期望的插入大小值中的大部分(例如,99.95%)。例如,其中插入大小被正态分布,该组距离可以是k的3个标准差内的那些距离。所确定的该组距离被用于选择下游节点,且因此选择候选路径。所选择的下游节点表示该对双端读数的第二成员的潜在比对位点。在图5中,该组下游节点还被填充在黑色方块中。开始于上游节点且结束于下游节点中的一个的路径是候选路径,将该对双端读数与候选路径比对。
[0119] 图6例示来自图4的DAG的候选路径。应当注意的是,图6还呈现包含在图4的DAG内的DAG。图6中例示了所有候选路径且仅例示了候选路径。可以将定义插入大小k(来自图4)的假设的双端读数对与图6中呈现的该DAG进行比对以节省计算量。在没有认识到对于比对的该约束的情况下,该对的每个读数可以沿着长度L的基因组并且在复杂的DAG(例如,图4中的DAG近似地为横越其长度4重的)中潜在地与任何位置比对,可存在多得多的位置(即,平均来说,如果DAG是横越任何基因位点n重的,则比对一对双端读数相当于从(nL)^2个选项中进行选择)。根据所例示的方法实行与DAG的比对可以使非常大的计划计算上成本小得多,产生具有来自RNA测序计划的高容量双端读数数据的大得多的吞吐量。
[0120] 4.对可变剪接的变种的敏感度。
[0121] 当读数跨越外显子的接合点时,当前技术使得难以或不可能在几种类型的情况之间进行区分。一种此情况是已发生突变的情况。即使在没有突变存在的情况下,当因为已经发生一种形式的可变剪接(其中密码子的一部分已被包含或跳过),所以读数可以显得与接合点的任一侧都不匹配时,出现另一种此情况。此外,测序误差呈现另一种此情况。
[0122] 考虑以下序列作为示例:
[0123] AUGCAUUUCCGUGGAUAGCGAUGCAUACUCACGAUG AUGAAAAUGCAUCAGAAAUAG(序列ID号:3)
[0124] 其中明文区域表示外显子,第一带下划线的部分表示内含子,并且第二带下划线的部分是由于可变剪接所以可能或可能不包含在第二外显子中的区域。因此,成熟RNA可采取以下两种形式中的任一种:
[0125] AUGCAUUUCCGUGGAUAGAUGAAAAUGCAUCAGAAAUAG(序列ID号:4)
[0126] AUGCAUUUCCGUGGAUAGAUGCAUCAGAAAUAG(序列ID号:5)
[0127] 现在假设我们获得以下读数:TGGATAGATGAA(序列ID号:6)。此读数可以具有在重要的方面不同的若干因果关系历史记录。
[0128] 图7呈现包括这些序列的三个可替代的历史记录。在横越图7的顶部第三个所描绘的第一情况中,采取来自序列ID号:4的读数,并且没有突变或测序误差已发生。在第二情况中,采取来自序列ID号:5的读数,并且测序误差已发生。在第三情况中,采取来自转录序列ID号:5的读数,并且读数精确地反映突变。
[0129] 如果线性参考序列仅由明文区域(图7中短划的下划线)组成,则很有可能比对评分将非常高——通过很多匹配升高且仅遭受一个不匹配罚分——使用中的算法将不使比对标记为可能错误。并非检查用于剪接信息的一些其它文件或数据库,而是将简单地假定突变或单个测序误差。与DAG参考进行比对的确较好。
[0130] 图8描绘了序列ID号:3与DAG参考(由于没有了内含子,所以现在出现在图8中作为序列ID号:4)的比对。这里,因为精确的比对将具有与该DAG任何可能的比对的最高评分,所以精确的比对立即被恢复。由于DAG的本质,图8还描绘了:
[0131] AUGCAUUUCCGUGGAUAGAUGCAUCAGAAAUAG(序列ID号:5)
[0132] 鉴于外显子的小部分遭受由可变剪接引起的包含或除去的频繁程度,这是很一般的现象并且将常常发生。当这些情况出现时,包括线性参考的系统将通常无法检测错误的特性描述,而防止这些错误的唯一方式是与如此多的线性参考进行比对,以至使过程过分地昂贵(在时间和金钱两方面)。
[0133] 5.正确地使用双端数据:改良DAG。
[0134] 本发明提供另外地利用插入大小和准确比对之间的关系。当处理插入大小作为对给定比对或一对比对的约束时,可以检视当前范例。这是对该情况的不完整的描述。将插入大小检视作为对于比对和转录本的组合的约束可以是更精确的。
[0135] 本发明的方法包括不仅使用插入大小作为对比对的品质控制,还使用插入大小作为对于转录本的品质控制。如果两个双端配偶中的每个在转录DAG中具有高评分比对,并且如果该对比对对应于非常不合情理的插入大小,则与关于比对中的每个的正确性调整我们的置信相比,关于转录调整置信(例如,通过将边缘添加到图表)常常将是更合理的。
[0136] 这可以通过——无论何时两个双端配偶与DAC比对——计算插入大小或由该比对暗示的大小来实施。(因为读数可能与若干特征一致,所以可能包括不止一个大小。)如果最小的此插入大小和插入大小的期望值的商大于一些常量(可能为2),则检查DAG以找出一对未连接的节点,如果该对未连接的节点连接,则将必然伴有异构体的存在,对于该异构体,该异构体的所暗示的插入大小和插入大小的预期值的商在0.5和1.5之间。如果此对节点存在,则用边缘连接它们。如果若干个此对存在,则选择产生最接近于1的比率的那个对并且连接它。
[0137] 在某些方面中,本发明提供一种分析转录组的方法,该方法包括从转录组获得至少一对双端读数、找出在该对的第一读数和在有向非循环数据结构(该数据结构具有表示RNA序列诸如外显子或转录本的节点和连接每对节点的边缘)中的节点之间的具有最佳评分的比对,识别包含节点的候选路径,该节点通过具有与该对双端读数的插入长度大体上类似的长度的路径连接到下游节点,和将双端读数与候选路径比对以确定最佳评分比对,以及改良DAG以在必要时包含新的边缘或者优选地创建由一对双端读数
指定的候选路径。本领域中的技术人员将认识到,可以由任何合适的指标定义具有与该对双端读数的插入长度大体上类似的长度的路径。例如,大体上类似可以意指相同数量的核苷酸,或者可以优选地意指+10或+100核苷酸。在某些实施例中,大体上类似可以由用户定义以意指在长度上彼此差距在10%内,或优选地彼此差距在25%内,或彼此差距在50%内,或彼此差距在100%内。在优选的实施例中,长度大体上类似意指一个的长度不超过另一个的长度的两倍。该方法适合于许多双端读数,诸如来自cDNA片段的那些读数(例如,表示单个转录组的cDNA片段)。节点和下游节点可以表示一对外显子,该对双端读数从该对外显子获得。方法可以包含基于所确定的最佳评分比对来识别异构体。在一些实施例中,所识别的异构体是新颖的异构体。可以更新有向非循环数据结构以表示新颖的异构体(例如,通过添加新的边缘)。方法可以包含从包括该对双端读数的任何比对计算中排除不是候选路径的任何路径。本发明的其它方面通过例如可用于推断异构体频率的系统和方法提供了表达分析。
[0138] 6.使用双端数据:推断异构体频率。
[0139] 强有力的技术不是起因于使用双端结构以正确地比对,而是起因于使用比对以用从比对推断插入大小的中间步骤推断异构体频率。如果多对外显子出现在样本中的频率是已知的(参见例如Harrow,2012,《基因组研究(Genome Research)》22:1760),则那些频率可以被用于构造异构体出现在样本中的频率。此外,不同的多对外显子通常将被不同的距离分隔开。因此,在双端配偶之间观测到的距离可以(稍微间接地)被用于
跟踪外显子频率。
[0140] 例示性示例如下。对于每对双端配偶,该对可以首先可选地与DAG参考进行比对。注意且跟踪每个对成员上的参考点(例如,最左边位置)之间的距离。这可以从与DAG的比对获得,虽然应当注意推断异构体频率的方法可以不需要将该对与DAG比对。
[0141] 用于注意和跟踪在每个被比对的对的参考点之间的距离的一个合适的方法是通过使用散列表。在距离散列中,关键字可以都是在1和1000之间的距离,并且每个关键字的值被初始化为0。然后,可以通过查找对应于距离的关键字的值并递增该值来完成递增。应注意,可以极其快速地完成散列表的这类运算——实际上,在常数时间内。另外,应该注意的是,通过现有的生物信息学
软件包诸如BioPerl很好地支持散列表。然后,每对外显子与对应于最接近外显子之间真实距离的距离的值的总和相关联。(也就是说:将为接近那些外显子之间距离的距离所观测到的所有“命中次数”加起来。在样本中,相对于的所有对的外显子,对的该数量和对的总数量的商将非常接近该对外显子发生的频率)。
[0142] 可从该数据用代数方法确定异构体频率。方程式的优选的形式可以取决于相关异构体的结构,但一般技术获得如下:
[0143] 假设有三个异构体,其中一个包含外显子A、B和C,其中另一个包含外显子A、B和D,并且其中第三个包含外显子B、C和D。外显子对的频率——其中[AB]是在所有的所比对的读数之间的外显子对AB的频率——从所比对的读数获得(即,在比对中,每当外显子i被连接到外显子j时,对应的频率递增1)。
[0144] 然后,如果三个异构体的频率是(按顺序)[A]=f,[B]=g和[C]=h,则外显子-对的频率应该如下:
[0145] [AB]=f+g
[0146] [W]=f+h
[0147] [X]=g+h
[0148] 当左侧值是从比对已知的时,可以解算异构体频率。一般来说,本发明的方法包含编写将外显子-对-频率与异构体-频率相关的方程式,并且使用前者以获得后者。本领域中的技术人员将认识到用于解算所描绘的方程式的代数方法。为了给出例示性示例,应指出上面的代数运算产生:
[0149] f=[AB]-g,
[0150] 以及:
[0151] h=[BD]-g
[0152] 因此:
[0153] [BC]=f+[BD]-g
[0154] 以及:
[0155] [BC]=[AB]-g+[BD]-g
[0156] 重排提供:
[0157] [BC]/2([AB]+[BD])=g
[0158] 由于[BC]、[AB]和[BD]都是从比对已知的,因此计算g。然后解算f和h是数学上最简单的。
[0159] 在一些情况下,可以出现“
自由度”问题,其中要解算的变量的数量不小于输入变量的数量。用于解决该显而易见的问题的任何合适的技术可以被采用,并且可以取决于特定的情况。可以找出的合适技术包含使变量标准化;解算超过一组异构体频率且事后选择一个;在外部获得至少一个输入值(例如,从生物样本获得),或者提供输入信息(例如,从文献或假设提供)。
[0160] 重要的是要注意关键步骤的图形解释:假设产生在x轴上具有距离和在y轴上具有频率的图表。然后,外显子-对将出现作为该图表上的峰值,并且经常那些峰值中的一些还将容易地可解释为异构体-频率。
[0161] 因此,本发明提供通过从转录组获得多个双端读数对且可选地将那些读数与参考数据结构进行比对,推断转录组中异构体频率的方法。数据结构包含来自基因组的每个被表示为节点的多个外显子。多对节点或外显子在它们的近端处通过边缘连接。确定双端读数的所比对的对之间的距离和所比对的对之间的距离的频率。确定一组异构体路径和异构体频率,以使得表示以异构体频率穿过结构的异构体路径引起多对外显子以所比对的对之间的距离的频率被包含在该异构体路径中。所确定的该组异构体路径中的每个路径可以表示存在于生物体中的转录异构体。该比对的对之间的距离可以是从每对的第一比对成员的上游端到该对的第二比对成员的下游端的路径长度。在某些实施例中,确定该组异构体路径包含发现新颖的异构体路径。可以更新参考数据结构以包含新颖的异构体路径。
[0162] 可以用代数方法例如通过解算一组线性方程式实行确定该组异构体路径。在一些实施例中,在解算该组线性方程式(例如,所以可以确定线性方程式中的其余的变量)之前,从生物学数据获得异构体频率中的一个。另外地或可替代地,在确定一组异构体路径和异构体频率之前,可以使在所比对的对之间的距离的频率标准化。
[0163] 7.合并先验外显子计数分布和异构体计数分布。
[0164] 在一些实施例中,本发明利用其中有用于异构体计数的相对频率的先验数据的情况。转录本的异构体分布可以被推断与蛋白质异构体的分布相关。即使不完全地,也已经彻底地研究了关于给定基因的外显子计数分布和异构体计数分布的问题,以及这些分布对关于序列的其它事实(例如,其是编码序列还是非编码序列)敏感的方式。
[0165] 可以从文献或者通过实行分析获得异构体频率。例如,可以从如在上面第6节“使用双端数据:推断异构体频率”中所描述的双端读数RNA测序读数推断异构体频率。在一些实施例中,来自双端读数数据的插入大小用于推断异构体频率,如下面第8节“使用插入大小更有效地推断异构体概率”中所描述。在某些实施例中,从源诸如数据库或文献获得异构体频率。
[0166] 图9表示由SwissProt
氨基酸序列的研究获得的蛋白质异构体分布,并且是来自Nakao等人,2005,《人类可替代的蛋白质异构体的大规模分析:模式分类和与亚细胞
定位信号的相关性(Large-scale analysis of human alternative protein isoforms:pattern classification and correlation with subcellular localization signals)》,《核酸研究(Nucl Ac Res)33(8):2355-2363的第三个图的再现。另外,Chang等人已经使用编程方法推导各种可变剪接形式的定量分布。Chang等人,2005,《来自EST数据库的可变剪接的定量分析中可变剪接图表的应用(The application of alternative splicing graphs in quantitative analysis of alternative splicing form from EST database)》,《计算机技术应用国际期刊(Int J.Comp.Appl.Tech)》22(1):14。
[0167] 由图9呈现的信息表示可以被用于辅助将读数映射到DAG的异构体计数分布。类似地,可以获得外显子计数分布信息。
[0168] 例如,有每转录本的外显子的数量的显示的分布。参见例如Harrow等人,2012,《GENCODE:用于编码计划的参考人类基因组注释(GENCODE:The reference human genome annotation for The ENCODE Project)》《,基因组研究(Genome Res)》22:1760-1774。蛋白质-编码基因位点处的转录本在每转录本4个外显子处显示峰值,而IncRNA在两个外显子处示出不同的峰值(参见Harrow等人,2012的图2)。可以用给出相对频率而不是绝对频率的y轴概率性地解释来自Harrow等人的外显子计数分布或来自Nakao等人的异构体分布此信息可以被用于告知关于候选比对对应于实际核苷酸序列的似然性的判定。
现有技术主要依赖于比对的评分。本发明提供包含异构体数量的先验知识的方式。
[0169] 例如,对于大的读数组,75%可以明确地映射到一组8个异构体{I1、I2、...、I8}中的区域,而其另外25%可以被解释为或者映射I5或者映射到对应于另一异构体I9的外显子的一些其它组合。在该情况下,各种类型的事实都与如何比对读数的最后25%的决策相关。此类相关事实包含(i)与I5比对和与I9比对的评分;(ii)我们可能具有的关于调控机制的任何进一步的知识;以及(III)有8个或9个异构体是否更可能为先验的。
[0170] 当将读数映射到DAG参考时,本发明合并异构体分布信息。可以以若干种不同的方式使用异构体频率分布信息。在映射读数中包含先验异构体数量的第一示例性方法包含:(1)为基因组的相关区域构造DAG(该DAG的边缘可以是加权的或未加权的)以及(2)将DAG的每个边缘与指示边缘的使用的值相关联。如果以跨越该边缘的方式所比对的读数的数量小于一些变量U,则设置使用“异常”。该变量直观地反映边缘是否已经被“使用”。估计该变量的一种方式可能是通过将其设置为等于(样本中读数的数量)*(外显子长度的合理估计和平均读数长度的商)*(.005)。
[0171] 进一步地,该方法包含(3)将基因组的相关区域与对应于使用
阈值T的变量相关联,可能设置T等于用于图表的边缘的数量的期望值加上用于图表的边缘的数量的标准差。这些数量可以通过以下步骤计算:估计DAG中边缘的数量,该DAG具有与为所讨论的基因而期望的经验外显子频率、异构体频率或这两者对应的节点和最大路径的数量。经验频率可以例如从Harrow等人,2012(用于外显子计数分布)或者从Nakao等人,2005(用于异构体分布)获得。
[0172] 该方法包含(4)改良用于史密斯-沃特曼或其它算法的评分方案;如果“使用的”边缘(即,布尔型使用值是“正常的”或“真的”边缘)的总数N超过T,则通过一些罚分惩罚未使用的边缘的遍历。可以使该罚分随着N增加而增加。用于罚分的合理的公式可以是max(0,(N-T)*(插入罚分)*2)。0项出现以避免负罚分——即,以避免使用新边缘的加分——但此些加分在一些情形中可能是合适的,在该情形中,不需要采取0的最大值和后者表达的值。
[0173] 为了避免呈比对对读数被处理的顺序敏感的形式的不被接受的路径依赖性,可以重复比对若干次且比较结果。这不影响算法的渐进运行时间。
[0174] 在映射读数中包含先验异构体数量的可替代的示例性方法包含:
[0175] (i)比对每个读数与DAG;
[0176] (ii)将所观测的异构体数量与期望的异构体数量(例如,给定关于基因的背景信息,从异构体数量的先验概率分布获得的)比较;以及
[0177] (iii)如果所观测的数量和所预期的数量之间的差值超过某一阈值,则重新比对被分配给最罕见的候选异构体的读数。在重新比对中,用根据差值的量值的权重将与罕见的候选异构体的任何比对罚分。罚分可以是max(0,(N-T)*(插入罚分)*2),其中N是“使用的”边缘(即,布尔型使用值是“正常的”或“真的”边缘)的总数,并且T是使用阈值。在一些实施例中,T是用于图表的边缘的数量的期望值加上用于图表的边缘的数量的标准差。
[0178] 通过应用使用异构体数量的先验概率分布的方法,将推动读数被映射到DAG以产生接近先验分布的分布。
[0179] 8.更有效地使用插入大小以推断异构体概率。
[0180] 假设在不同的异构体中出现两个外显子E1和E2,并且在那两个异构体中外显子E1和E2通过距离L1和L2被分隔开。将偶尔发生获得一对双端读数,其一个配偶映射到E1,并且另一个映射到E2。然后,将该对归属于哪个异构体的选择可以被缩减为鉴于配偶的从外显子边缘的偏移是读数之间的真实距离,是否推断L1或L2或那些距离中的一个移位适当的量的问题。
[0181] 因为我们通常具有与双端测序运行中插入大小的分布有关的良好的背景信息,所以将异构体-分配问题转换为关于插入大小的问题是有价值的。鉴于通常将没有其它理由相信更有可能从一个异构体或另一个异构体采集读数,我们可以实行简单的贝叶斯计算,以根据分布中插入大小L1和L2的概率的比率将读数概率性地分配给两个异构体。
[0182] 这是对当前技术的显著改进,其该情况下总体上检视读数为携载用于两个异构体的等量证据权重。因此,以所描述的方式合并插入大小是不忽略相关证明的方式。
[0183] 因此,在这些不明确的情况下,我们可以计数异构体,如下:
[0184] (1)计算插入大小的概率分布。(由P(s)指代具有插入大小s的读数的概率。)应注意这仅需要计算一次,并且这可以在预处理步骤中完成。
[0185] (2)如果两个双端配偶中的每个映射在外显子内部:
[0186] (2a)制作包含那些外显子的所有异构体的列表。
[0187] (2b)将每个异构体与将由已经从该异构体导出的读数暗示的插入大小相关联。
[0188] (2c)为(2b)中所计算的插入大小中的每个计算P(s)。
[0189] (2d)将每个异构体与P(s)除以(2c)中所计算的所有P(s)值的总和的相关值相关联。
[0190] 应该注意的是,本文所描述的技术将以不影响使用“下游”工具(也就是说,将读数计数作为输入且给出其它有用的输出的工具)的能力的方式改进异构体的读数计数。
[0191] 9.构造“最大DAG”。
[0192] 在一些实施例中,构造DAG可以是有益的,其中,在此DAG中每对节点用从表示“较早”或“左侧”外显子到“较迟”或“右侧”外显子的节点的边缘连接。
[0193] 使用Gencode或类似的经注释的转录组,创建节点以表示外显子和其它特征。Gencode提供足够信息以将这些节点与此信息相关联作为节点的序列、其在参考序列中占据的位置(并且因此其长度)等。因此,使用包括耦接到非暂时性存储器的处理器的计算机系统,来自基因组的多个特征中的每个被表示为节点。对于该多个特征中的每对,创建边缘,边缘通过该对的两个成员的近端连接该两个成员。此移动的生物学解释是猜测表示遵守外显子出现在基因组中的顺序的那些外显子的每个组合的异构体的存在,或者至少促进发现该异构体。对遵守外显子出现在基因组中的顺序敏感既有感应性的原因又有生物化学的原因。
[0194] 例如,报告表明在多外显子前mRNA中的外显子总是按顺序保持。参见例如Black,2005《,剪接难题的简单答案(A simple answer for a splicing conundrum)》,《美国国家科学院院刊(PNAS)》102:4927-8。没有被任何特定的机制束缚,从理论上可以讲,剪接机械设备辨识内含子的开始和结束,以及分
支点腺嘌呤或接近剪接接合点的其它保守序列。剪接机械设备或剪接体沿着前mRNA从5'到3'移动,移除内含子且剪接外显子。因此看来存在用于以5'到3'顺序连接所有外显子的现象学以及文献
支撑,在5'到3'顺序中,所有的外显子看起来沿着基因组DNA链。还参见Shao等人,2006,《人类中外显子重复接收、外显子加扰和反式剪接的生物信息学分析(Bioinformatic analysis of exon repetition,exon scrambling and trans-splicing in humans)》,《生物信息学(Bioinformatics)》22:692-
698。对于某些假设,本发明的方法遵循人类前mRNA被剪接成保留由基因组序列定义的外显子顺序的线性分子的假设。
[0195] 另外,在期望在有向非循环数据结构中包含外显子加扰的实例(参见例如Carrington等人,1985,《在伴刀豆球蛋白A的转译后修饰期间发生多肽接合(Polypeptide ligation occurs during post-translational modification of concanavalin A)》,《自然(Nature)》313:64-67)的情况下,可以通过包含双重复外显子作为“虚构”大体上保持外显子以它们的原始顺序剪接的假设(或者,在可替代的方案中,可以废除该假设)。
[0196] 因此,在例如基因组包含处于以下顺序:ABCDEFGHI的以下外显子的情况下[0197] 使用短划线作为边缘且字母作为节点,完整的该组边缘可以被表示为:
[0198] A-C;A-D;A-E;A-F;A-G;A-H;A-I;B-C;B-D;B-E;B-F;B-G;B-H;B-I;C-D;C-E;C-F;C-G;C-H;C-I;D-E;D-F;D-G;D-H;D-I;E-F;E-G;E-H;E-I;F-G;F-H;F-I;G-H;G-I。
[0199] 果这足以证明发现所有的异构体而不是一种,并且该一种异构体——给出示例:ADEHG
[0200] 则,借助仿真的或虚构的G——称它为G',该非循环数据结构可以是“钉齿状的”。然后,可以由A-D-E-H-G'表示该异构体,并且结构不需要包含从H到G的任何边缘。然而,可以发现可以忽略外显子加扰和由所有包含的边缘保持的基因组外显子顺序。
[0201] 可以通过使用来自经注释的转录组的输入数据,首先制备本发明的DAG。同样地,从头开始,DAG将包含所有已知的异构体。如果制成最大的DAG,则DAG将包含保持基因组外显子顺序的所有可能的边缘,并且然后,将因此包含既表示已知的异构体又表示新颖的异构体的节点和边缘。在存储器中节点和边缘被存储为有向非循环数据结构,并且方法包含将双端序列读数与包括由边缘的子集连接的特征的子集的路径进行比对,从而识别对应于路径的异构体。
[0202] 当将新的数据与DAG(例如,双端RNA测序读数)比对时,那些数据可以包含在经注释的参考中未被识别的外显子。因为在该示例中外显子是新颖的,所以该外显子是异构体。使用DAG作为参考理想地适合于该情况,因为将新颖的异构体与现有的DAG进行比对还可以引导在DAG中创建新的节点的以表示新颖的外显子。因此,在所识别的异构体是新颖的情况下,比对转录序列可以包括在有向非循环数据结构中创建新节点。还可以为剪接变异体创建新节点,例如,其中一个外显子通过单个密码子不同于另一个外显子。
[0203] 优选地,在分析双端读数之前且独立于分析双端读数,初始地全部或部分地构建最大DAG。因此,方法可以包含在比对转录本之前,为多个特征中的每对创建边缘。出于与有向非循环数据结构的本质相关的原因,数据结构中表示的参考可以任意大的,例如,数百或数千外显子、内含子或两者。
[0204] 10.本发明的系统
[0205] 图10例示用于实施本文所描述的方法的计算机系统401。本发明的系统可以包含图10中所示的任一个部件或任何数量的部件。一般来说,系统401可以包含能够通过网络415彼此通信的计算机433和
服务器计算机409。另外,可以可选地从数据库405(例如,本地或远程)获得数据。在一些实施例中,系统包含用于获得测序数据的仪器455,仪器455可以被耦接到测序仪计算机451用于序列读数的初始处理。
[0206] 在一些实施例中,通过并行处理实行方法,并且服务器409包含具有并行架构的多个处理器,即,能够将多个序列(例如,RNA测序读数)与参考序列构造(例如,DAG)进行比较的处理器和存储器的分布式网络。系统包括多个处理器,这些处理器同时执行多个读数与参考序列构造之间的多个比较。虽然可能有其它混合配置,但是并行计算机中的主存储器通常或者在单个地址空间中的所有处理元件之间共享,或者呈分布式,即,每个处理元件具有其自身的本地地址空间。(分布式存储器是指这样的事实:存储器以逻辑方式分布,但通常暗示它也以物理方式分布。)分布式共享存储器和存储器虚拟化组合这两种方法,其中处理元件具有其自身的本地存储器以及对非本地处理器上的存储器的存取权。对本地存储器的存取通常比对非本地存储器的存取更快。
[0207] 其中可以用相等时延和带宽存取主存储器的每个元件的计算
机架构被称为均匀存储器存取(UMA)系统。通常,UMA系统只能通过共享的存储器系统来实现,其中存储器并非以物理方式分布。不具有该特性的系统被称为不均匀存储器存取(NUMA)架构。分布式存储器系统具有不均匀存储器存取。
[0208] 可以以若干种方式在
硬件中实施处理器-处理器和处理器-存储器通信,这些方式包含经由共享的(或者多端口的或者多路复用的)存储器、纵横
开关、共享的总线或无数拓扑的互连网络(包含星形、环形、树形、超立方体、胖超立方体(在一个节点处具有不止一个处理器的超立方体))或n维网格。
[0209] 基于互连网络的并行计算机必须合并路由以实现在并非直接连接的节点之间的消息传递。用于处理器之间的通信的媒体很可能在大型多处理器机器中分层。此类资源在市面上可购买用于专用用途,或可以经由例如亚马逊云计算的“云”存取这些资源。
[0210] 计算机总体上包含通过总线耦接到存储器和输入-输出(I/O)机构的处理器,存储器可以包含RAM或ROM,并且优选地包含存储指令的至少一个有形的非暂时性媒体,该指令可执行以致使系统实行本文所描述的功能。本领域中的技术人员将认识到,根据需要或者最适合于实行本发明的方法,本发明的系统包含一个或多个处理器(例如,中央处理单元(CPU)、
图形处理单元(GPU)等)、计算机可读存储设备(例如,主存储器、静态存储器等),或其组合,它们通过总线彼此通信。
[0211] 处理器可以是本领域中已知的任何合适的处理器,诸如由英特尔(加利福尼亚州圣克拉拉(Santa Clara,CA))以商标XEON E7出售的处理器,或由AMD(加利福尼亚州桑尼维尔(Sunnyvale,CA))以商标OPTERON 6200出售的处理器。
[0212] 根据本发明的输入/输出设备可以包含视频显示单元(例如,
液晶显示器(LCD)或
阴极射线管(CRT)监视器)、字母数字输入设备(例如,
键盘)、
光标控制设备(例如,
鼠标或触控板)、磁盘
驱动器单元、信号生成设备(例如,扬声器)、
触摸屏、
加速计、麦克
风、蜂窝无线电频率天线,以及网络
接口设备,网络接口设备可以是例如网络接口卡(NIC)、Wi-Fi卡或蜂窝
调制解调器。
[0213] 以引用的方式并入
[0214] 贯穿本公开已经参考并且引用了其它文档,诸如专利、专利申请、专利公开案、期刊、书籍、论文、网络内容。所有此类文档在此以全文引用的方式并入本文中用于所有目的。
[0215] 等效物
[0216] 除本文示出且描述的那些之外,本领域的技术人员将从本文档的完整内容(包含对在本文中引用的科学和专利文献的参考)中明白本发明的各种改良及本发明的许多进一步的实施例。本文中的主题含有重要信息、范例和指南,它们可适于本发明在本发明的各种实施例及其等效物中的实践。