CN105830078A

CN105830078A - 用于在有向非循环结构中使用双端数据的系统和方法

Info

Publication number: CN105830078A
Application number: CN201480067827.2A
Authority: CN
Inventors: 丹尼斯·库拉尔; 南森·梅维斯
Original assignee: Seven Bridges Genomics Inc
Current assignee: Seven Bridges Genomics Inc
Priority date: 2013-10-21
Filing date: 2014-10-15
Publication date: 2016-08-03
Anticipated expiration: 2034-10-15
Also published as: EP3061022A1; WO2015061099A1; KR20160073406A; AU2014340461B2; EP3061022B1; US20150112658A1; WO2015061103A1; US20150310167A1; JP2016539443A; SG11201603044SA; US9063914B2; CA2927723C; US9092402B2; CN105830078B; US20150302145A1; US10055539B2; CA2927723A1; US20150112602A1; US10204207B2

Abstract

本申请涉及一种分析转录组的方法，该方法包括从来自生物体的转录组获得至少一对双端读数，找出在该对的第一读数和有向非循环数据结构(该数据结构具有表示RNA序列诸如外显子或转录本的节点和连接多对节点的边缘)中的节点之间的具有最佳评分的比对，识别包含节点的候选路径，该节点通过具有与该对双端读数的插入长度大体上类似的长度的路径连接到下游节点，并且将双端读数与候选路径比对以确定最佳评分比对。

Description

用于在有向非循环结构中使用双端数据的系统和方法

相关申请的交叉参考

本申请要求2014年1月17日提交的美国专利申请序列号14/157,979的权益和优先权，该美国专利申请要求2013年10月21日提交的美国临时专利申请序列号61/893,467的权益和优先权，该临时专利申请的内容以引用方式并入本文。

序列表

本申请包含已通过网络电子提交系统(EFS-Web)以ASCII格式提交的序列表，并且以其全部内容以引用方式并入本文。2014年10月14日创建的ASCII格式化的序列表被命名为SBG-009-01WO-2-Sequences_ST25，并且大小为1,217字节。

技术领域

本发明大体上涉及生物信息学，并且尤其涉及使用有向非循环数据结构进行分析。

背景技术

对于特定的发育阶段或生理状态，转录组是细胞中完整的一组转录本和它们的数量。转录组对生物体的功能、发育和疾病极为重要。生物体的本质和活力来源于所有的转录本(包含mRNA、非编码RNA和小RNA)、它们的表达水平、基因剪接模式和转录后修饰。事实上，人类基因组中与期望相比大得多的部分现在已知为被转录的。参见Bertone等人，2004，《具有基因组嵌合阵列的人类转录序列的全局识别(Globalidentificationofhumantranscribedsequenceswithgenometilingarrays)》，《科学(Science)》306：2242-2246。

已经报道基于下一代测序(NGS)技术的用于转录组分析的方法。参见，例如，Wang等人，2009，《RNA测序：用于转录组学的革命性的工具(RNA-Seq:arevolutionarytoolfortranscriptomics)》，《遗传学自然评论(NatRevGenet)》10(l)：57-632009。RNA测序(将RNA测序)方法承诺快速生成高容量的转录组数据。

然而，RNA测序面临信息挑战，诸如存储和分析大量数据，这必须被解决以更好地使用读数。比对和分析RNA测序读数不仅在涉及的信息的本质上呈现出问题，还在容量上呈现出问题。具有大数据集的读映射是计算上的挑战，并且用于差异表达的分析方法才开始出现。参见Garber等人，2011，《使用RNA测序的用于转录组注释和量化的计算方法(ComputationalmethodsfortranscriptomeannotationandquantificationusingRNA-Seq)》，《自然方法(NatMeth)》8(6)：469-477。映射来自RNA测序的大量读数的计算需求被将来可获得的大量参考数据大大倍增。

例如，GENCODE联合体力图识别已经被报道的人类基因组中的所有基因特征。Harrow等人，2012，《GENCODE：用于编码计划的参考人类基因组注释(GENCODE:ThereferencehumangenomeannotationforTheENCODEProject)》，《基因组研究(GenomeRes)》22：1760-1774。由GENCODE计划包含的材料的量是庞大的。不仅仅要继续添加新的蛋白质编码基因位点，而且经注释的可变剪接转录本的数量稳定地增加。GENCODE7版本包含超过20,000个蛋白质-编码和几乎10,000个长的非编码RNA基因位点(IncRNA)，以及未在其它源中展现的超过33,000个编码转录本。GENCODE还包含其它特征诸如非转录区域(UTR)、长基因间非编码RNA(lincRNA)基因、短非编码RNA和可变剪接模式。即使RNA测序研究开始仅具有有限量的新数据，来自源诸如GENCODE的潜在参考数据的容量也是如此大，使得理解新数据在计算上是具有挑战性的。

发明内容

本发明大体上提供了用于使用组织为有向非循环图(DAG)或类似的数据结构的参考对双端读数进行分析的系统和方法。参考中的特征诸如外显子和内含子在DAG中提供节点，并且那些特征通过边缘以它们的规范的基因顺序彼此链接。DAG可以按比例绘制为任何大小，并且事实上，在第一实例中可以被填入来自外来经注释的参考的导入序列。在新数据包含双端读数的情况下，一对中的每个读数可以潜在地沿长度L的基因组且在复杂的DAG中与任何位置比对，可存在多得多的位置(即，平均来说，如果DAG是横越任何基因位点n重的，则比对一对双端读数相当于从(nL)^2个选项中进行选择)。然而，插入长度信息提供了用于将潜在的比对约束在DAG内的手段——不仅是在路径内远离的距离中，还通过排除路径、将选项的数量限制为相对于L是线性的。双端读数还可以提供其它信息。例如，双端读数可以提供用于修饰DAG的基础，可以被用于推断异构体频率、或者其它。事实上，具有插入长度信息的双端读数集携载信息的品质，该品质非常适合于聚焦基于DAG的比对和更新DAG以表示可能没有被发现或注意到的实际上出现的固有异构体。因此，双端读数特别适合于RNA测序计划，允许相对于很大的参考集分析很大的RNA测序数据集。使用双端RNA测序读数发现转录组的所有的异构体并不是计算上难处理的，并且转录组学可以以先前不可能的步伐和吞吐量移动。可以从一个或任何数量的生物体研究一个或任何数量的转录组以增加对生物体功能、发育和疾病的理解。

本发明的某方面涉及创建包含参考的特征(例如，外显子、内含子等)的DAG，参考的特征以基因组顺序通过它们的近端边缘连接。该DAG数据结构具有转录组学中的应用，包含例如识别异构体。此DAG参考被用于涉及包含配对读数的RNA测序数据的转录组学中。可以使用包含作为节点的特征(例如从经注释的源基因组或转录组创建的)的DAG参考来分析RNA测序数据。本发明的系统和方法提供了使用双端读数的改进的比对、使用双端读数信息构建DAG、推断异构体概率或频率，或其组合。

在某些方面中，本发明提供了分析转录组的方法，该方法包括从来自生物体的转录组获得至少一对双端读数，找出在该对的第一读数和有向非循环数据结构中的节点之间的具有最佳评分的比对(该数据结构具有表示RNA序列诸如外显子或转录本的节点和连接多对节点的边缘)，识别包含节点的候选路径，该节点通过具有与该对双端读数的插入长度大体上类似的长度的路径连接到下游节点，以及将双端读数与候选路径比对以确定最佳评分比对。方法可以包含从包括该对双端读数的任何比对计算中排除不是候选路径的任何路径。

该方法适合于任何数量的双端读数。节点和下游节点可以表示一对外显子，该对双端读数从该对外显子获得。该方法可以包含基于所确定的最佳评分比对来识别异构体。在一些实施例中，所识别的异构体是新颖的异构体。有向非循环数据结构可以被更新以表示新颖的异构体(例如，通过添加至少一个新节点或边缘)。

本发明的相关方面提供用于转录组学的计算机系统。该系统具有耦接到非暂时性存储器的处理器，并且被用于储存有向非循环数据结构，该有向非循环数据结构具有表示RNA序列的节点和连接多对节点的边缘。系统可操作以从转录组获得一对双端读数，在该对的第一读数和数据结构中的节点之间找出具有最佳评分的比对，识别包含节点的候选路径，该节点通过具有与该对双端读数的插入长度大体上类似的长度的路径连接到下游节点，以及将该双端读数与该候选路径比对以确定最佳评分比对。优选地，系统从随后的比对尝试中排除所有非候选路径。在某些实施例中，该系统进一步可操作以从cDNA片段中获得多个双端读数。优选地，系统可操作以为多个双端读数对中的每个确定最佳评分比对。cDNA片段可以表示单个转录组。系统可以基于所确定的最佳评分比对来识别异构体，诸如新颖的异构体。系统可操作以更新有向非循环数据结构以表示新颖的异构体。在一些实施例中，更新有向非循环数据结构以表示新颖的异构体包括添加至少一个新节点。节点和下游节点可以表示一对外显子，双端读数对从该外显子对获得。系统可操作以从包括该对双端读数的任何比对计算中排除不是候选路径的任何路径。

本发明的其它方面通过例如可用于推断异构体频率的系统和方法提供了表达分析。

在某些方面中，本发明提供推断转录组中异构体频率的方法。该方法包含从生物体的转录组获得多个双端读数对，并且可选地将那些读数与参考数据结构进行比对。数据结构包含来自基因组的每个被表示为节点的多个特征。多对节点或特征在它们的近端处通过边缘连接(即，在对内，两个对成员在基因组上在它们的近端处通过边缘彼此连接；由于任何对成员还可以是其它对的成员，所以特征的任一端或两端可以通过一个或任何数量的边缘被连接到其它特征)。确定双端读数的所比对的对之间的距离和该所比对的对之间的距离的频率。确定一组异构体路径和异构体频率，以使得表示以异构体频率穿过结构的异构体路径引起多对特征以所比对的对之间的距离的频率被包含在该异构体路径中。所确定的该组异构体路径中的每个路径可以表示存在于生物体中的转录异构体。所比对的对之间的距离可以是从每对的第一比对成员的上游端到该对的第二比对成员的下游端的路径长度。在某些实施例中，确定该组异构体路径包含发现新颖的异构体路径。可以更新参考数据结构以包含新颖的异构体路径。

可以用代数方法例如通过解算一组线性方程式实行确定该组异构体路径。在一些实施例中，在解算该组线性方程式(例如，所以可以确定线性方程式中的其余的变量)之前，从生物学数据获得异构体频率中的一个。另外地或可替代地，在确定一组异构体路径和异构体频率之前，可以使所比对的对之间的距离的频率标准化。

可以从现有的基因特征数据库构建数据结构。例如，每个节点可以表示来自此数据库的特征。

在相关的方面中，本发明提供推断转录组中异构体频率的计算机系统。系统包含耦接到非暂时性存储器的处理器，并且存储参考数据结构，该参考数据结构包括来自基因组的每个被表示为节点的多个特征，其中该多个特征的对在它们的近端处通过边缘连接。系统可以被用于获得或确定来自转录组的多个双端读数对之间的距离和该对之间的距离的频率。然后可以确定一组异构体路径和异构体频率，以使得表示以异构体频率穿过结构的异构体路径引起多对特征以所比对的对之间的距离的频率被包含在该异构体路径中。特征可以包含外显子，并且系统进一步可操作以将多个双端读数对与参考数据结构进行比对。所确定的该组异构体路径中的每个路径可以表示存在于生物体中的转录异构体。确定该组异构体路径可以包含发现新颖的异构体路径。系统可以更新参考数据结构以包含新颖的异构体路径。所比对的对之间的距离可以包括从每对的第一比对成员的上游端到该对的第二比对成员的下游端的路径长度。在一些实施例中，确定该组异构体路径包括解算一组线性方程式。在某些实施例中，确定该组异构体路径包括在解算该组线性方程式之前，从生物学数据获得异构体频率中的一个。系统可操作以在确定该组异构体路径和异构体频率之前，使在所比对的对之间的距离的频率标准化。数据结构可以包含表示远程基因特征数据库中的每个经注释的特征的节点。

附图说明

图1示出与DAG比较的读数的图形表示。

图2示出对应于比较的实际矩阵。

图3提供了本发明的方法的图式。

图4图示了假设的DAG和双端读数对。

图5示出穿过假设的DAG的候选路径。

图6示出来自假设的DAG的候选路径。

图7呈现包括序列的三个可替代的历史记录。

图8根据本发明描绘了比对。

图9示出异构体频率分布。

图10例示了用于实施本发明的方法的系统。

具体实施方式

应当理解，一般来说，基因不是用于代码单个、特定蛋白质。代替地，在不同的时间以不同的比例转录不同子组的给定的基因的外显子。由这些不同的转录本创建的不同的蛋白质被称为异构体，并且相关地，不同的转录本本身也被称为异构体。如上面所论述的，现代生物信息学已经提供了潜在地用作参考的至少部分地经注释的转录组。例如，Gencode计划已经不仅编录基因，还编录人类基因组中的外显子(连同其它区域)。

本文所描述的系统和方法被用于有效地利用那些经注释的转录组中的信息。实际上，本发明包含将双端读数与经注释的转录组进行比较理想地适合于分析方法的洞见，在该分析方法中序列数据被存储为基于DAG的数据结构。

本发明的一些实施例提供用经注释的转录组中的数据创建有向非循环数据结构的方法。在此数据结构中，来自经注释的转录组的特征(例如，外显子和内含子)被表示为节点，该节点通过边缘连接。因此，由穿过基于DAG的数据结构的对应路径表示来自现实世界转录组的异构体。

本发明的某方面涉及创建包含来自一个或多个已知参考的特征诸如内含子和外显子的DAG。在本领域中DAG被理解为是指可以被呈现为图表的数据以及呈现数据的图表。本发明提供用于将DAG存储为可以由计算机系统读取以用于生物信息学处理或用于呈现为图表的数据的方法。可以以任何合适的格式保存DAG，该合适格式包含例如节点和边缘的列表、矩阵或表示矩阵的表、一批阵列或表示矩阵的类似的变量结构，以内置有图表语法的语言保存，以用于图表目的的通用标记语言保存，或其它合适格式。

在一些实施例中，DAG被存储为节点和边缘的列表。此方式的其中一种是创建包含所有节点和所有边缘的文本文件，该所有节点具有分配给每个节点的ID，该所有边缘中的每个具有开始和结束节点的节点ID。因此，例如，如果为两个语句“SeeJanerun”和“Run,Janerun”创建DAG，则可以创建不分大小写的文本文件。可以使用任何合适的格式。例如，文本文件可包含逗号分隔值。将该DAG命名为“Jane”用于将来的参考，以该格式，DAG“Jane”可以读作如下：1see、2run、3jane、4run，1-3、2-3、3-4。本领域的技术人员将理解，该结构可适用于图7和图8中所表示的内含子和外显子，以及下面的随附论述。

在某些实施例中，DAG被存储为表示矩阵(或一批阵列或表示矩阵的类似的变量结构)的表，其中N×N矩阵中的(i,j)项指示节点i和节点j被连接(其中N是以基因组顺序包含节点的矢量)。为了DAG是非循环的，仅需要所有的非零项在对角线上面(假定节点以基因组顺序表示)。在二进制情况下，0项表示从节点i到节点j不存在边缘，并且1项表示从i到j有一个边缘。本领域的技术人员将理解矩阵结构允许除0到1之外的值与边缘相关联。例如，任何项可以是指示权重或所使用的次数，反映世界中所观测的数据的一些固有品质的数值。矩阵可以被写入到文本文件作为表或一系列线性的行(例如，首先是行1，接着是分隔符等)，因此提供简单的串行化结构。

在定义节点之后，使矩阵DAG串行化的一种有用的方式将是为项使用逗号分隔值。使用该格式，DAG“Jane＂”将包含与上述相同的节点定义，接着是矩阵项。该格式可以读作：

1see、2run、3jane、4run

,,1,\,,1,\,,,1\,,,

其中简单地省略了零(0)的项，并且“/”是换行符。

本发明的实施例包含以内置有图表语法的语言存储DAG。例如，具有被称为图表可视化(Graphviz)的图表可视化软件包的DOT语言提供了可以被用于存储具有辅助信息且可以使用从Graphviz网站商购获得的多个工具转化成图表文件格式的数据结构。Graphviz是开源图表可视化软件。图表可视化是一种将结构信息表示为抽象的图表和网络的图式的方式。其具有在联网、生物信息学、软件工程、数据库和网页设计、机器学习中的应用，以及在用于其它技术领域的可视界面中的应用。Graphviz布局规划采取以简单文本语言为形式的图表描述，并且以有用的格式(诸如用于网页的图像和SVG；用于包含在其它文档中的PDF或附言；或者在交互式图表浏览器中的显示)制作图式。

在相关的实施例中，DAG以用于图表目的的通用标记语言或其它格式被存储。根据上面线性文本文件或逗号分隔矩阵的描述，本领域的技术人员将认识到，语言诸如XML可以被用于(扩展到)创建定义节点和它们的报头或ID、边缘、权重等的标签(标记)。然而，DAG被结构化且被存储，本发明的实施例包括使用节点表示特征诸如外显子和内含子。这为分析RNA测序读数与发现、识别和表示异构体提供了有用的工具。

如果节点表示特征或特征的片段，则可以由穿过这些片段的路径表示异构体。由将一些先前的外显子连接到一些之后的外显子的边缘表示外显子是“被跳过的”。本文所呈现的是用于构造DAG以表示基因的可变剪接或异构体的技术。在Lee和Wang，2005，《可变剪接的生物信息学分析(Bioinformaticsanalysisofalternativesplicing)》，《简明生物信息学(BriefBioinf)》6(l)：23-33；Heber等人，2002，《剪接图表和EST装配问题(SplicinggraphsandESTassemblyproblems)》，《生物信息学(Bioinformatics)》18Suppl：sl81-188；Leipzig等人，2004，《可变剪接巷道(ASG)：桥接基因组和转录组之间的空位(Thealternativesplicinggallery(ASG):Bridgingthegapbetweengenomeandtranscriptome)》，《核酸研究(NuclAcRes)》23(13)：3977-2983；以及LeGault和Dewey，2013，《来自具有概率性的剪接图表的RNA测序数据的可变剪接的推断(InferenceofalternativesplicingfromRNA-Seqdatawithprobabilisticsplicegraphs)》，《生物信息学(Bioinformatics)》29(18)：2300-2310中论述了可变剪接，这些文献中的每个文献的内容以引用方式并入。可以在Florea等人，2005，《具有AIR的基因和可变剪接注释(GeneandalternativesplicingannotationwithAIR)》，《基因组研究(GenomeResearch)》15：54-66；Kim等人，2005，《ECgene：用于可变剪接的基于基因组的EST群集和基因建模(ECgene:Genome-basedESTclusteringandgenemodelingforalternativesplicing)》，《基因组研究(GenomeResearch)》15：566-576；以及Xing等人，2006，《用于来自剪接图表的全长异构体的概率性重构的期望-最大化算法(Anexpectation-maximizationalgorithmforprobabilisticreconstructionsoffull-lengthisoformsfromsplicegraphs)》，《核酸研究(NucleicAcidsResearch)》，34，3150-3160中找出另外的论述，这些文献中的每个文献的内容以引用方式并入。

1.使用双端读数

本发明涉及包括RNA测序数据的转录组学，该RNA测序数据包含双端读数。本发明的系统和方法提供使用双端读数的改进的比对。在某些实施例中，双端读数被用于约束比对，由此大大降低RNA测序计划的计算成本。

很多测序数据是双端数据。存在其中此类数据的特征应该在比对中体现他们自身的方式：双端配偶在比对中被分配的位置之间的距离应该大致对应于由测序运行给定的插入大小。一般来说，由插入大小给定的信息涉及比对算法试图揭示的信息。

用于产生成对的核酸读数的若干方法在本领域中是已知的。通常，该方法包括选择性地将核酸样本片段化从而得到已知大小分布的序列的一些形式。通常用官能团或通过自环化来保护片段化的样本的端部，于是剩余的核酸样本材料被移除。然后，端部被扩增且测序，导致被认为彼此之间产生一些距离的核酸读数的群体。

配偶对样本制备成套装置是例如从NexteraIllumina(加利福尼亚州圣迭戈(SanDiego，CA))商购的。为了制备成对的核酸用于在此平台上测序，核酸样本被片段化成长度在2kb和5kb之间的段。然后，使用亲和柱分隔，片段端被生物素标记以促进之后的恢复。生物素标记的端部被接合以创建环状DNA段。此时，来自样本的非环化的核酸被消化且随后被移除。然后，剩余的环化的DNA被再片段化以产生适合于群集和测序的片段。可以用亲和柱或链霉亲和素包被的小珠从环化的DNA的剩余片段移除被生物素标记的原始端。然后，恢复的端进行端修复、聚腺苷酸加尾，并且在扩增和测序之前添加测序衔接子。

2.比对

本发明提供用于将序列读数与表示经注释的参考的有向非循环数据结构进行比对的方法和系统。使用本发明的比对算法，即使存在与RNA测序结果相关联的大数量和大体上全面的参考，读数也可以被快速地映射。与在其它情况下相同，DAG是比对的理想目标。与之后在某个比对的工作流程中使用经注释的转录组相比，实行与由经注释的转录组构造的DAG进行比对是有效得多的。

通过使用DAG作为参考获得大量益处。例如，和与一个参考进行比对，并且然后试图根据经注释的转录组调整某个比对的结果相比，与DAG进行比对更精确。这主要是因为前一方法在用于初始比对的序列和转录组中所表示的另一序列之间强加了不自然的不对称。和试图与每种物理可能性(此类可能性的数量将在数个接合点中总体上指数增长)的线性序列进行比对相比，与潜在地表示所有的相关物理可能性的目标进行比对是计算上有效得多的。

本发明的实施例包含将一个或多个读数与参考数据结构进行比对。该比对可以包含基于DAG的比对，其找出最优地表示读数的数据结构的边缘和节点。该比对还可以包含以成对的方式将读数与来自DAG的节点的外显子比对的实例。

成对比对总体上包括沿目标的一部分放置一个序列，根据算法引入空位，对两个序列匹配的程度评分，和优选地沿着参考对不同的位置重复进行该成对比对。最优评分匹配被认为是比对，并且表示关于序列数据表示的内容的推断。在一些实施例中，对一对核酸序列的比对进行评分包括为置换和插入缺失的概率设置值。当比对了单个碱基时，匹配或不匹配以置换概率贡献于比对评分，置换概率可以是例如对于匹配为1以及对于不匹配为-0.33。插入缺失以将空位罚分从比对评分中扣除，空位罚分可以是例如-1。空位罚分和置换概率可以基于关于序列如何演化的经验知识或先验假设。它们的值影响得到的比对。特别地，空位罚分和置换概率之间的关系影响置换或插入缺失是否在得到的比对中将是有利的。

正式地说，比对表示两个序列x和y之间的推断关系。例如，在一些实施例中，序列x和y的比对A将x和y分别映射到可以包含空隙的另外两个字符串x'和y'，使得：(i)lx'l＝ly'l；(ii)从x'和y'移除空隙应该分别回到x和y；以及(III)对于任何i，x'[i]和y'[i]无法两个都是空隙。

空位是x'或y'中的任一者中的连续空隙的最大子串。比对A可以包含以下三种区域：(i)匹配的对(例如，x'[i]＝y'[i])；(ii)不匹配的对，(例如，x'[i]≠y'[i]，并且两者均不是空隙)；或(III)空位(例如，或者x'[i..j]或者y'[i..j]是空位)。在某些实施例中，仅匹配的对具有高的正评分a。在一些实施例中，不匹配的对总体上具有负评分b，并且长度r的空位还具有负评分g+rs，其中g，s<0。对于DNA，一个通用评分方案(例如，由BLAST所使用的)使得评分a＝l、评分b＝-3、g＝-5和s＝-2。比对A的评分是所有的匹配的对、不匹配的对和空位的评分的总和。x和y的比对评分可以被定义为在x和y的所有可能的比对之中的最大评分。

在一些实施例中，任何对具有由置换概率的4×4矩阵B定义的评分a。例如，B(i,i)＝1和0<B(i,j)i<>j<1是一个可能的评分系统。例如，在与颠换相比变换被认为是更加生物学上可能的情况下，矩阵B可包含B(C,T)＝7和B(A,T)＝3，或者期望的或由本领域中已知的方法确定的任何其它组值。

根据本发明的一些实施例，比对包含成对比对。一般来说，成对比对包括——具有m个字符的序列Q(查询)和n个字符的参考基因组T(目标)——找出和评估Q和T之间的可能的局部比对。对于任何l<i<n和l<j<m，计算T[h..i]和Q[k..j]的最大可能比对评分(即，在位置i处结束的T的任何子串和在位置j处结束的Q的任何子串的最优比对评分)，其中h<i且k<j。这可以包含检查所有具有cm个字符的子串，其中c是根据相似模型的常量，并且将每个子串单独与Q比对。每个比对被评分，并且具有优选的评分的比对被接受为比对。本领域中的技术人员将理解，存在序列比对的精确算法和近似算法。精确算法将找出最高评分的比对，但是在计算上会昂贵的。两种众所周知的精确算法是尼德曼-翁施算法(Needleman-Wunsch)(《分子生物学杂志(JMolBiol)》，48(3)：443-453，1970)和史密斯-沃特曼算法(Smith-Waterman)(《分子生物学杂志(JMolBiol)》，147(1)：195-197，1981；《数学进展(Adv.inMath.)》20(3)，367-387，1976)。后藤(Gotoh)(《分子生物学杂志(JMolBiol)》，162(3)，705-708，1982)对史密斯-沃特曼算法的进一步改进将计算时间从O(m2n)减少到O(mn)，其中m和n是被比较的序列大小，该改进更可修正为并行处理。在生物信息学领域，正是后藤的改良算法通常被称为史密斯-沃特曼算法。史密斯-沃特曼方法用于将较大序列集与较大参考序列进行比对，因为可更普遍并且更便宜地获得并行计算资源。参见例如亚马逊(Amazon)的云计算资源。本文所参考的所有期刊文章以其全部内容以引用的方式并入。

史密斯-沃特曼(SW)算法通过奖励序列中的碱基之间的重叠并且对序列之间的空位罚分来比对线性序列。史密斯-沃特曼算法还与尼德曼-翁施算法不同，不同之处在于SW不要求较短序列跨越描述较长序列的字母字符串。也就是说，SW不假定一个序列是另一个序列的全部内容的读数。此外，因为SW并不一定找出横跨字符串全长的比对，所以局部比对可以在两个序列内的任何地方开始和结束。

在一些实施例中，根据点矩阵法、动态规划法或单字法，进行成对比对。动态规划法一般实施史密斯-沃特曼(SW)算法或尼德曼-翁施(NW)算法。根据NW算法的比对总体上根据具有线性空位罚分d的相似矩阵S(a,b)(例如，诸如前述矩阵B)对比对的字符进行评分。矩阵S(a,b)总体上提供置换概率。SW算法类似于NW算法，但是任何负评分矩阵单元被设置为0。在美国专利5,701,256和美国公布2009/0119313中更详细地描述了SW算法和NW算法及其实施方式，该专利和公布以其全部内容以引用方式并入本文。

实施史密斯-沃特曼算法的一个版本的比对程序是MUMmer，MUMmer可以从由Geeknet(弗吉尼亚州费尔法克斯(Fairfax，VA))维护的SourceForge网站商购获得。MUMmer是用于快速比对基因组规模序列的系统(Kurtz、S等人，《基因组生物学(GenomeBiology)》，5：R12(2004)；Delcher、A.L.等人，《核酸研究(Nucl.AcidsRes.)》，27：11(1999))。例如，MUMmer3.0可以在2.4GHzLinux台式计算机上使用78MB存储器在13.7秒内找出在一对5-兆碱基基因组之间的所有20-碱基对或更长的精确匹配。MUMmer可以处理来自鸟枪法测序计划的100个或1000个重叠群，并且将使用系统包含的NUCmer程序将100个或1000个重叠群与另一组重叠群或参考进行比对。如果对于DNA序列比对来说物种太相异而不能检测相似性，则PROmer程序可以基于两者输入序列的六框翻译生成比对。

其它示例性比对程序包含：核苷酸数据库的高效大规模比对(EfficientLarge-ScaleAlignmentofNucleotideDatabases(ELAND))或序列和变异的同感评估(ELANDv2componentoftheConsensusAssessmentofSequenceandVariation(CASAVA))软件(加利福尼亚州圣迭戈的Illumina(Illumina,SanDiego,CA))的ELANDv2部件；RealTimeGenomics公司的实施基因组学(RTG)调查器(加利福尼亚州旧金山(SanFrancisco,CA))；Novocraft的Novoalign(马来西亚雪兰莪州(Selangor,Malaysia))；Exonerate，欧洲生物信息研究所(EuropeanBioinformaticsInstitute)(英国辛克斯顿(Hinxton,UK))(Slater，G.和Birney,E.，《BMC生物信息学(BMCBioinformatics)》6：31(2005))，都柏林大学的ClustalOmega，(爱尔兰都柏林(Dublin,Ireland))(SieversF等人，《分子系统生物学(MolSystBiol)》7，article539(2011))；都柏林大学的ClustalW或ClustalX(爱尔兰都柏林(Dublin,Ireland))(LarkinM.A等人，《生物信息学(Bioinformatics)》，23，2947-2948(2007)；和FASTA，欧洲生物信息研究所(EuropeanBioinformaticsInstitute)(英国辛克斯顿(Hinxton,UK))(PearsonW.R等人，《美国国家科学院院刊(PNAS)》85(8)：2444-8(1988)；Lipman,D.J.，《科学(Science)》227(4693)：1435-41(1985))。

如上面所论述的，当将RNA测序读数与有向非循环经注释的参考基因组进行比对时，实施SW比对算法或(下面进一步更详细地论述的)的改良的版本可以是优选的或期望的。

根据以下方程式(1)，对于表示长度n和m的两个字符串的n×m矩阵H，易于表达SW算法：

H_k0＝H_ol＝0(对于0≤k≤n且0≤l≤m)(1)

H_ij＝max{H_i-1,j-1+s(a_i,b_j),H_i-1-W_in,H_i,j-1-W_del,0}

(对于1≤i≤n且1≤j≤m)

在以上方程式中，s(a_i,b_j)表示匹配奖分(当a_i＝b_j时)或不匹配罚分(当a_i≠b_j时)，并且对插入和缺失分别给出罚分W_in和W_del。在大多数情况下，所得矩阵具有为零的许多元素。这种表示使得更容易在矩阵中从高到低、从右到左回溯，因此识别比对。

在已经用评分完全填充矩阵后，SW算法实行回溯以确定比对。开始于矩阵中的最大值，算法将基于三个值(H_i-1,j-1、H_i-1,j或H_i,j-1)被用于计算每个单元的最终最大值来进行回溯。当达到零时回溯停止。最佳评分比对可以包含大于插入和缺失的最小可能数量的可能数量，同时包含远小于置换的最大可能数量的可能数量。

当作为SW或SW-后藤应用时，该技术使用动态规划算法来实行分别具有大小m和n的两个字符串S和A的局部序列比对。此动态规划技术采用表或矩阵来保存匹配评分并避免对连续单元的重新计算。可以根据序列的字母索引字符串的每个元素，也就是说，如果S是字符串ATCGAA，则S[l]＝A。

代替将最佳比对表示为H_i,j(上面)，在以下方程式(2)中将最佳比对表示为B[j,k]：

B[j,k]＝max(p[j,k],i[j,k],d[j,k],0)(对于0<j≤m，0<k≤n)(2)

在以下方程式(3)-(5)中概述最大值函数B[j,k]的变量参数，其中MISMATCH_PENALTY、MATCH_BONUS、INSERTION_PENALTY、DELETION_PENALTY和OPENING_PENALTY都是常量，并且除MATCH_BONUS以外都是负数。通过以下方程式(3)给出匹配变量参数p[j,k]：

如果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)

如果S[j]＝A[k]，则p[j,k]＝max(p[[j-1,k-],i[j-1,k-1],d[j-1,k-1])+MATCH_BONUS,

通过以下方程式(4)给出插入变量参数i[j,k]：

i[j,k]＝max(p[j-l,k]+OPENING_PENALTY,i[j-l,k],d[j-l,k]+(4)

OPENING_PENALTY)+INSERTION_PENALTY

并且缺失变量参数d[j,k]通过以下方程式(5)给出：

d[j,k]＝max(p[j,k-1]+OPENING_PENALTY,i[j,k-1]+(5)

OPENING_PENALTY,d[j,k-1])+DELETION_PENALTY

对于所有三个变量参数，将[0,0]元素设置为零以确保回溯完成，即，p[0,0]＝i[0,0]＝d[0,0]＝0。

评分参数是稍微任意的，并可经调整以实现计算的行为。DNA的评分参数设置的一个示例(Huang，《第3章：生物序列比较和比对(Bio-SequenceComparisonandAlignment)》，当前顶级比较分子生物学(CurrTopCompMolBiol.)丛书，马萨诸塞州剑桥市(Cambridge,Mass.)：麻省理工学院出版社(TheMITPress)，2002)将是：

MATCH_BONUS：10

MISMATCH_PENALTY：-20

INSERTION_PENALTY：-40

OPENING_PENALTY：-10

DELETION_PENALTY：-5

以上空位罚分(INSERTION_PENALTY、OPENING_PENALTY)之间的关系有助于限制空位开放的数目，即，支持通过设置高于空位开放成本的空位插入罚分来将空位集合在一起。当然，MISMATCH_PENALTY、MATCH_BONUS、INSERTION_PENALTY、OPENING_PENALTY和DELETION_PENALTY之间可能存在可替代的关系。

在一些实施例中，本发明的方法及系统并入多维比对算法。本发明的多维算法提供了序列信息的“向后看”类型分析(如在史密斯-沃特曼中)，其中通过包含多个通路和多个节点的多维空间进行向后看。多维算法可以被用于将RNA测序读数与DAG类型参考进行比对。该比对算法通过关于包含在DAG(例如，参考序列构造)上的位置处的每个序列识别最大评分为C_i,j识别最大值。事实上，通过在先前位置处“向后”看，有可能横越多个可能的路径识别最佳比对。

在上面所论述的读数(亦称“字符串”)和有向非循环图表(DAG)上执行本发明的算法。出于定义该算法的目的，假设S是要比对的字符串，并且假设D是S将进行比对的有向非循环图表。以从1开始的索引对字符串S的元素加括号。因此，如果S是字符串ATCGAA，那么S[1]＝A、S[4]＝G等。

在某些实施例中，对于DAG，节点的序列的每个字母将被表示为独立元素d。d的前趋被定义为：

(i)如果d不是其节点的序列的首字母，则其节点中在d之前的字母是其(唯一)前趋；

(ii)如果d是其节点的序列的首字母，则为d的节点的父节点的任何节点(例如，基因组中上游的所有外显子)的序列的最后一个字母是d的前趋。

所有前趋的组继而表示为P[d]。

为了找出“最优”比对，算法寻求M[j,d](S的前j个元素与在d之前(且包含d)的DAG的一部分的最佳比对的评分)的值。该步骤类似于在以上的方程式1中发现H_i,j。具体来说，确定M[j,d]包括找出a、i、e和0的最大值，如下面所定义：

M[j,d]＝max{a,i,e,0}(6)

其中

对于在P[d]中的p*，e＝max{M[j,p*]+DELETE_PENALTY}

i＝M[j-1,d]+INSERT_PENALTY

对于在P[d]中的p*，如果S[j]＝d，则a＝max{M[j-1,p*]+MATCH_SCORE}；

对于在P[d]中的p*，如果S[j]≠d，则max{M[j-1,p*]+MISMATCH_PENALTY}

如上面所描述，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个字符比对的最高比对评分的前趋相同。

在方程式(6)中，i是字符串S的前j-1个字符与DAG的直到并且包含d的部分的比对加上INSERT_PENALTY，其类似于SW中的插入变量参数的定义(参看方程式1)。

另外，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。

再次，与在SW算法中相同，罚分例如DELETE_PENALTY、INSERT_PENALTY、MATCH_SCORE和MISMATCH_PENALTY可以被调整以促进与较少空位等的比对。

如以上方程式中所描述，该算法通过不仅计算该元素的插入、缺失和匹配评分，而且向后看(逆着DAG的方向)到DAG上的任何之前节点以找出最大评分，来找出每个读数的最佳(即，最大)值。因此，该算法能够详细研究穿过DAG的含有已知突变的不同路径。因为图表是有向的，所以逆着图表方向移动的回溯跟随朝向图表起点的优选异构体，并且最佳比对评分识别高度确定性的最可能比对。

图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两者。如所描绘的，序列读数被示出为与上部路径比对。

图2示出对应于比较的实际矩阵。相似于史密斯-沃特曼技术，本发明的所例示算法识别最高评分，并且实行回溯以识别读数的正确位置。图1和图2还突出本发明产生字符串与该构造的实际匹配。在其中序列读数包含未被包含在DAG中的变异体的情况下，将通过空位、插入等报告经比对的序列。

将RNA测序读数与经注释的参考的DAG进行比对允许发现转录组中的异构体。另外，已经发生的不同形式的可变剪接的可能性(例如，一个密码子已在相关外显子中的一个中被跳过的可能性)可以被“免费”预期。

外显子可自动地分成两个节点，一个对应于第一密码子，并且另一个对应于外显子的其余部分。这允许在DAG中表示两种类型的剪接，即使还未观测到它们。鉴于一个密码子的与比对的差值的价值将往往不足以惩罚候选比对以诱导我们寻找替代方式，预测这些情况可以是有用的。(最大的此罚分将是三个插入或三个缺失的罚分，并且由于机率，所以常常罚分将是较小的。)

将转录组表示为DAG准许将异构体自然地解释为穿过DAG的路径。与往常一样，当系统和方法表示数学结构中的物理项时，益处是大量的。例如，本发明的系统和方法允许序列以直观的方式可视化。例如，可以通过示出穿过DAG的两个路径来比较序列。

另外，本发明的方法有益于利用图表定理和算法。穿过DAG的路径的分析是在现代的数学和计算机科学中充满生机的研究区域。因此，利用该研究来改进转录组研究的能力是有价值的。例如，如果边缘被给定与由读数跨越的接合点的概率对应的权重，则路径的总权重可以携载关于实现不同异构体的先验可能性的信息。穿过DAG的路径最大化(或最小化)权重是已知的，并且最短路径算法可以被用于快速地找出那些最短路径。

图3提供了本发明的方法的图式。如上面所论述的，一旦DAG内置有以规范的顺序连接的外显子，则获得双端读数。然后，双端读数与DAG比对，并且在读数和包括表示RNA序列的节点和连接多对节点的边缘的DAG之间找出比对(具有满足预定指标的比对评分)。预定指标可以是例如最高评分比对。比对可以通过上述改良的史密斯-沃特曼算法进行。对于一对双端读数，比对可以包含基于穿过包含比对的有向非循环数据结构的路径，识别转录组中的转录本。通过使用插入大小作为比对的约束，可以有效地实行比对。双端读数可以被用于构建DAG(例如，添加边缘)。

该方法适合于分析由RNA测序获得的读数。在某些实施例中，有向非循环数据结构基本上表示至少一个染色体的所有已知的特征。可以在找出比对之前创建和存储有向非循环数据结构。在一些实施例中，通过将每个序列读数与穿过有向非循环数据结构的可能路径中的至少大部分进行比较来找出比对。方法可以包含基于找出的比对将多个序列读数装配成重叠群。

3.使用插入大小作为比对的约束

其中插入大小信息可以与DAG组合的一种方式是使用插入大小作为比对的约束。通过使用该方法，速度和精确性都有可能增加。因为该方法防止丢弃与比对相关的重要信息，所以与单独比对双端配偶相比，该方法可以更精确。应当注意的是，独立比对两个双端配偶近似地相当于在基因组中L个可能的位置(如果基因组具有长度L)之间选择两个位置。因此，比对包括在大约L^2种可能性之间进行选择。然而，第一比对被用于强有力地将另一比对约束到在第一比对周围的常量大小k的一些区域，然后仅有大约kL种可能性——在L而不是L^2的量级上。这将更加快速。

图4-图6例示使用插入大小信息以约束比对。

图4例示任意结构的假设的DAG和定义插入大小k的一对假设的双端读数。在如所绘制的DAG中，节点是方块，并且边缘是线。在所描绘的实施例中，从具有均值kbp的近似正态分布获取插入大小k。将该对的第一成员与DAG比对(例如，使用上述史密斯-沃特曼或改良的史密斯-沃特曼)。一般来说，双端包括片段，该片段的每端上具有衔接子，并且该片段具有从片段-衔接子边界朝向彼此朝内延伸的双端读数。在双端读数之间可存在空位。插入长度可以意味着片段的长度——即，从一个片段-衔接子边界到另一个片段-衔接子边界的距离。参见例如Lindgreen，2012，《AdapterRemoval：下一代序列读数的简单清理(AdapterRemoval:easycleaningofnext-generationsequencereads)》，《BMC研究纪要《(BMCResNotes)》5：337(例如，Lindgreen，图1，画面C)，其内容以引用方式并入。

在图5中，最左边黑色方块表示外显子，该对的第一成员在外显子处比对。确定了一组距离，它们占期望的插入大小值中的大部分(例如，99.95％)。例如，其中插入大小被正态分布，该组距离可以是k的3个标准差内的那些距离。所确定的该组距离被用于选择下游节点，且因此选择候选路径。所选择的下游节点表示该对双端读数的第二成员的潜在比对位点。在图5中，该组下游节点还被填充在黑色方块中。开始于上游节点且结束于下游节点中的一个的路径是候选路径，将该对双端读数与候选路径比对。

图6例示来自图4的DAG的候选路径。应当注意的是，图6还呈现包含在图4的DAG内的DAG。图6中例示了所有候选路径且仅例示了候选路径。可以将定义插入大小k(来自图4)的假设的双端读数对与图6中呈现的该DAG进行比对以节省计算量。在没有认识到对于比对的该约束的情况下，该对的每个读数可以沿着长度L的基因组并且在复杂的DAG(例如，图4中的DAG近似地为横越其长度4重的)中潜在地与任何位置比对，可存在多得多的位置(即，平均来说，如果DAG是横越任何基因位点n重的，则比对一对双端读数相当于从(nL)^2个选项中进行选择)。根据所例示的方法实行与DAG的比对可以使非常大的计划计算上成本小得多，产生具有来自RNA测序计划的高容量双端读数数据的大得多的吞吐量。

4.对可变剪接的变种的敏感度。

当读数跨越外显子的接合点时，当前技术使得难以或不可能在几种类型的情况之间进行区分。一种此情况是已发生突变的情况。即使在没有突变存在的情况下，当因为已经发生一种形式的可变剪接(其中密码子的一部分已被包含或跳过)，所以读数可以显得与接合点的任一侧都不匹配时，出现另一种此情况。此外，测序误差呈现另一种此情况。

考虑以下序列作为示例：

AUGCAUUUCCGUGGAUAGCGAUGCAUACUCACGAUGAUGAAAAUGCAUCAGAAAUAG(序列ID号：3)

其中明文区域表示外显子，第一带下划线的部分表示内含子，并且第二带下划线的部分是由于可变剪接所以可能或可能不包含在第二外显子中的区域。因此，成熟RNA可采取以下两种形式中的任一种：

AUGCAUUUCCGUGGAUAGAUGAAAAUGCAUCAGAAAUAG(序列ID号：4)

AUGCAUUUCCGUGGAUAGAUGCAUCAGAAAUAG(序列ID号：5)

现在假设我们获得以下读数：TGGATAGATGAA(序列ID号：6)。此读数可以具有在重要的方面不同的若干因果关系历史记录。

图7呈现包括这些序列的三个可替代的历史记录。在横越图7的顶部第三个所描绘的第一情况中，采取来自序列ID号：4的读数，并且没有突变或测序误差已发生。在第二情况中，采取来自序列ID号：5的读数，并且测序误差已发生。在第三情况中，采取来自转录序列ID号：5的读数，并且读数精确地反映突变。

如果线性参考序列仅由明文区域(图7中短划的下划线)组成，则很有可能比对评分将非常高——通过很多匹配升高且仅遭受一个不匹配罚分——使用中的算法将不使比对标记为可能错误。并非检查用于剪接信息的一些其它文件或数据库，而是将简单地假定突变或单个测序误差。与DAG参考进行比对的确较好。

图8描绘了序列ID号：3与DAG参考(由于没有了内含子，所以现在出现在图8中作为序列ID号：4)的比对。这里，因为精确的比对将具有与该DAG任何可能的比对的最高评分，所以精确的比对立即被恢复。由于DAG的本质，图8还描绘了：

AUGCAUUUCCGUGGAUAGAUGCAUCAGAAAUAG(序列ID号：5)

鉴于外显子的小部分遭受由可变剪接引起的包含或除去的频繁程度，这是很一般的现象并且将常常发生。当这些情况出现时，包括线性参考的系统将通常无法检测错误的特性描述，而防止这些错误的唯一方式是与如此多的线性参考进行比对，以至使过程过分地昂贵(在时间和金钱两方面)。

5.正确地使用双端数据：改良DAG。

本发明提供另外地利用插入大小和准确比对之间的关系。当处理插入大小作为对给定比对或一对比对的约束时，可以检视当前范例。这是对该情况的不完整的描述。将插入大小检视作为对于比对和转录本的组合的约束可以是更精确的。

本发明的方法包括不仅使用插入大小作为对比对的品质控制，还使用插入大小作为对于转录本的品质控制。如果两个双端配偶中的每个在转录DAG中具有高评分比对，并且如果该对比对对应于非常不合情理的插入大小，则与关于比对中的每个的正确性调整我们的置信相比，关于转录调整置信(例如，通过将边缘添加到图表)常常将是更合理的。

这可以通过——无论何时两个双端配偶与DAC比对——计算插入大小或由该比对暗示的大小来实施。(因为读数可能与若干特征一致，所以可能包括不止一个大小。)如果最小的此插入大小和插入大小的期望值的商大于一些常量(可能为2)，则检查DAG以找出一对未连接的节点，如果该对未连接的节点连接，则将必然伴有异构体的存在，对于该异构体，该异构体的所暗示的插入大小和插入大小的预期值的商在0.5和1.5之间。如果此对节点存在，则用边缘连接它们。如果若干个此对存在，则选择产生最接近于1的比率的那个对并且连接它。

在某些方面中，本发明提供一种分析转录组的方法，该方法包括从转录组获得至少一对双端读数、找出在该对的第一读数和在有向非循环数据结构(该数据结构具有表示RNA序列诸如外显子或转录本的节点和连接每对节点的边缘)中的节点之间的具有最佳评分的比对，识别包含节点的候选路径，该节点通过具有与该对双端读数的插入长度大体上类似的长度的路径连接到下游节点，和将双端读数与候选路径比对以确定最佳评分比对，以及改良DAG以在必要时包含新的边缘或者优选地创建由一对双端读数指定的候选路径。本领域中的技术人员将认识到，可以由任何合适的指标定义具有与该对双端读数的插入长度大体上类似的长度的路径。例如，大体上类似可以意指相同数量的核苷酸，或者可以优选地意指+10或+100核苷酸。在某些实施例中，大体上类似可以由用户定义以意指在长度上彼此差距在10％内，或优选地彼此差距在25％内，或彼此差距在50％内，或彼此差距在100％内。在优选的实施例中，长度大体上类似意指一个的长度不超过另一个的长度的两倍。该方法适合于许多双端读数，诸如来自cDNA片段的那些读数(例如，表示单个转录组的cDNA片段)。节点和下游节点可以表示一对外显子，该对双端读数从该对外显子获得。方法可以包含基于所确定的最佳评分比对来识别异构体。在一些实施例中，所识别的异构体是新颖的异构体。可以更新有向非循环数据结构以表示新颖的异构体(例如，通过添加新的边缘)。方法可以包含从包括该对双端读数的任何比对计算中排除不是候选路径的任何路径。本发明的其它方面通过例如可用于推断异构体频率的系统和方法提供了表达分析。

6.使用双端数据：推断异构体频率。

强有力的技术不是起因于使用双端结构以正确地比对，而是起因于使用比对以用从比对推断插入大小的中间步骤推断异构体频率。如果多对外显子出现在样本中的频率是已知的(参见例如Harrow，2012，《基因组研究(GenomeResearch)》22：1760)，则那些频率可以被用于构造异构体出现在样本中的频率。此外，不同的多对外显子通常将被不同的距离分隔开。因此，在双端配偶之间观测到的距离可以(稍微间接地)被用于跟踪外显子频率。

例示性示例如下。对于每对双端配偶，该对可以首先可选地与DAG参考进行比对。注意且跟踪每个对成员上的参考点(例如，最左边位置)之间的距离。这可以从与DAG的比对获得，虽然应当注意推断异构体频率的方法可以不需要将该对与DAG比对。

用于注意和跟踪在每个被比对的对的参考点之间的距离的一个合适的方法是通过使用散列表。在距离散列中，关键字可以都是在1和1000之间的距离，并且每个关键字的值被初始化为0。然后，可以通过查找对应于距离的关键字的值并递增该值来完成递增。应注意，可以极其快速地完成散列表的这类运算——实际上，在常数时间内。另外，应该注意的是，通过现有的生物信息学软件包诸如BioPerl很好地支持散列表。然后，每对外显子与对应于最接近外显子之间真实距离的距离的值的总和相关联。(也就是说：将为接近那些外显子之间距离的距离所观测到的所有“命中次数”加起来。在样本中，相对于的所有对的外显子，对的该数量和对的总数量的商将非常接近该对外显子发生的频率)。

可从该数据用代数方法确定异构体频率。方程式的优选的形式可以取决于相关异构体的结构，但一般技术获得如下：

假设有三个异构体，其中一个包含外显子A、B和C，其中另一个包含外显子A、B和D，并且其中第三个包含外显子B、C和D。外显子对的频率——其中[AB]是在所有的所比对的读数之间的外显子对AB的频率——从所比对的读数获得(即，在比对中，每当外显子i被连接到外显子j时，对应的频率递增1)。

然后，如果三个异构体的频率是(按顺序)[A]＝f，[B]＝g和[C]＝h，则外显子-对的频率应该如下：

[AB]＝f+g

[W]＝f+h

[X]＝g+h

当左侧值是从比对已知的时，可以解算异构体频率。一般来说，本发明的方法包含编写将外显子-对-频率与异构体-频率相关的方程式，并且使用前者以获得后者。本领域中的技术人员将认识到用于解算所描绘的方程式的代数方法。为了给出例示性示例，应指出上面的代数运算产生：

f＝[AB]-g，

以及：

h＝[BD]-g

因此：

[BC]＝f+[BD]-g

以及：

[BC]＝[AB]-g+[BD]-g

重排提供：

[BC]/2([AB]+[BD])＝g

由于[BC]、[AB]和[BD]都是从比对已知的，因此计算g。然后解算f和h是数学上最简单的。

在一些情况下，可以出现“自由度”问题，其中要解算的变量的数量不小于输入变量的数量。用于解决该显而易见的问题的任何合适的技术可以被采用，并且可以取决于特定的情况。可以找出的合适技术包含使变量标准化；解算超过一组异构体频率且事后选择一个；在外部获得至少一个输入值(例如，从生物样本获得)，或者提供输入信息(例如，从文献或假设提供)。

重要的是要注意关键步骤的图形解释：假设产生在x轴上具有距离和在y轴上具有频率的图表。然后，外显子-对将出现作为该图表上的峰值，并且经常那些峰值中的一些还将容易地可解释为异构体-频率。

因此，本发明提供通过从转录组获得多个双端读数对且可选地将那些读数与参考数据结构进行比对，推断转录组中异构体频率的方法。数据结构包含来自基因组的每个被表示为节点的多个外显子。多对节点或外显子在它们的近端处通过边缘连接。确定双端读数的所比对的对之间的距离和所比对的对之间的距离的频率。确定一组异构体路径和异构体频率，以使得表示以异构体频率穿过结构的异构体路径引起多对外显子以所比对的对之间的距离的频率被包含在该异构体路径中。所确定的该组异构体路径中的每个路径可以表示存在于生物体中的转录异构体。该比对的对之间的距离可以是从每对的第一比对成员的上游端到该对的第二比对成员的下游端的路径长度。在某些实施例中，确定该组异构体路径包含发现新颖的异构体路径。可以更新参考数据结构以包含新颖的异构体路径。

可以用代数方法例如通过解算一组线性方程式实行确定该组异构体路径。在一些实施例中，在解算该组线性方程式(例如，所以可以确定线性方程式中的其余的变量)之前，从生物学数据获得异构体频率中的一个。另外地或可替代地，在确定一组异构体路径和异构体频率之前，可以使在所比对的对之间的距离的频率标准化。

7.合并先验外显子计数分布和异构体计数分布。

在一些实施例中，本发明利用其中有用于异构体计数的相对频率的先验数据的情况。转录本的异构体分布可以被推断与蛋白质异构体的分布相关。即使不完全地，也已经彻底地研究了关于给定基因的外显子计数分布和异构体计数分布的问题，以及这些分布对关于序列的其它事实(例如，其是编码序列还是非编码序列)敏感的方式。

可以从文献或者通过实行分析获得异构体频率。例如，可以从如在上面第6节“使用双端数据：推断异构体频率”中所描述的双端读数RNA测序读数推断异构体频率。在一些实施例中，来自双端读数数据的插入大小用于推断异构体频率，如下面第8节“使用插入大小更有效地推断异构体概率”中所描述。在某些实施例中，从源诸如数据库或文献获得异构体频率。

图9表示由SwissProt氨基酸序列的研究获得的蛋白质异构体分布，并且是来自Nakao等人，2005，《人类可替代的蛋白质异构体的大规模分析：模式分类和与亚细胞定位信号的相关性(Large-scaleanalysisofhumanalternativeproteinisoforms:patternclassificationandcorrelationwithsubcellularlocalizationsignals)》，《核酸研究(NuclAcRes)33(8)：2355-2363的第三个图的再现。另外，Chang等人已经使用编程方法推导各种可变剪接形式的定量分布。Chang等人，2005，《来自EST数据库的可变剪接的定量分析中可变剪接图表的应用(TheapplicationofalternativesplicinggraphsinquantitativeanalysisofalternativesplicingformfromESTdatabase)》，《计算机技术应用国际期刊(IntJ.Comp.Appl.Tech)》22(1)：14。

由图9呈现的信息表示可以被用于辅助将读数映射到DAG的异构体计数分布。类似地，可以获得外显子计数分布信息。

例如，有每转录本的外显子的数量的显示的分布。参见例如Harrow等人，2012，《GENCODE：用于编码计划的参考人类基因组注释(GENCODE:ThereferencehumangenomeannotationforTheENCODEProject)》，《基因组研究(GenomeRes)》22：1760-1774。蛋白质-编码基因位点处的转录本在每转录本4个外显子处显示峰值，而IncRNA在两个外显子处示出不同的峰值(参见Harrow等人，2012的图2)。可以用给出相对频率而不是绝对频率的y轴概率性地解释来自Harrow等人的外显子计数分布或来自Nakao等人的异构体分布此信息可以被用于告知关于候选比对对应于实际核苷酸序列的似然性的判定。现有技术主要依赖于比对的评分。本发明提供包含异构体数量的先验知识的方式。

例如，对于大的读数组，75％可以明确地映射到一组8个异构体{I1、I2、...、I8}中的区域，而其另外25％可以被解释为或者映射I5或者映射到对应于另一异构体I9的外显子的一些其它组合。在该情况下，各种类型的事实都与如何比对读数的最后25％的决策相关。此类相关事实包含(i)与I5比对和与I9比对的评分；(ii)我们可能具有的关于调控机制的任何进一步的知识；以及(III)有8个或9个异构体是否更可能为先验的。

当将读数映射到DAG参考时，本发明合并异构体分布信息。可以以若干种不同的方式使用异构体频率分布信息。在映射读数中包含先验异构体数量的第一示例性方法包含：(1)为基因组的相关区域构造DAG(该DAG的边缘可以是加权的或未加权的)以及(2)将DAG的每个边缘与指示边缘的使用的值相关联。如果以跨越该边缘的方式所比对的读数的数量小于一些变量U，则设置使用“异常”。该变量直观地反映边缘是否已经被“使用”。估计该变量的一种方式可能是通过将其设置为等于(样本中读数的数量)*(外显子长度的合理估计和平均读数长度的商)*(.005)。

进一步地，该方法包含(3)将基因组的相关区域与对应于使用阈值T的变量相关联，可能设置T等于用于图表的边缘的数量的期望值加上用于图表的边缘的数量的标准差。这些数量可以通过以下步骤计算：估计DAG中边缘的数量，该DAG具有与为所讨论的基因而期望的经验外显子频率、异构体频率或这两者对应的节点和最大路径的数量。经验频率可以例如从Harrow等人，2012(用于外显子计数分布)或者从Nakao等人，2005(用于异构体分布)获得。

该方法包含(4)改良用于史密斯-沃特曼或其它算法的评分方案；如果“使用的”边缘(即，布尔型使用值是“正常的”或“真的”边缘)的总数N超过T，则通过一些罚分惩罚未使用的边缘的遍历。可以使该罚分随着N增加而增加。用于罚分的合理的公式可以是max(0，(N-T)*(插入罚分)*2)。0项出现以避免负罚分——即，以避免使用新边缘的加分——但此些加分在一些情形中可能是合适的，在该情形中，不需要采取0的最大值和后者表达的值。

为了避免呈比对对读数被处理的顺序敏感的形式的不被接受的路径依赖性，可以重复比对若干次且比较结果。这不影响算法的渐进运行时间。

在映射读数中包含先验异构体数量的可替代的示例性方法包含：

(i)比对每个读数与DAG；

(ii)将所观测的异构体数量与期望的异构体数量(例如，给定关于基因的背景信息，从异构体数量的先验概率分布获得的)比较；以及

(iii)如果所观测的数量和所预期的数量之间的差值超过某一阈值，则重新比对被分配给最罕见的候选异构体的读数。在重新比对中，用根据差值的量值的权重将与罕见的候选异构体的任何比对罚分。罚分可以是max(0，(N-T)*(插入罚分)*2)，其中N是“使用的”边缘(即，布尔型使用值是“正常的”或“真的”边缘)的总数，并且T是使用阈值。在一些实施例中，T是用于图表的边缘的数量的期望值加上用于图表的边缘的数量的标准差。

通过应用使用异构体数量的先验概率分布的方法，将推动读数被映射到DAG以产生接近先验分布的分布。

8.更有效地使用插入大小以推断异构体概率。

假设在不同的异构体中出现两个外显子E1和E2，并且在那两个异构体中外显子E1和E2通过距离L1和L2被分隔开。将偶尔发生获得一对双端读数，其一个配偶映射到E1，并且另一个映射到E2。然后，将该对归属于哪个异构体的选择可以被缩减为鉴于配偶的从外显子边缘的偏移是读数之间的真实距离，是否推断L1或L2或那些距离中的一个移位适当的量的问题。

因为我们通常具有与双端测序运行中插入大小的分布有关的良好的背景信息，所以将异构体-分配问题转换为关于插入大小的问题是有价值的。鉴于通常将没有其它理由相信更有可能从一个异构体或另一个异构体采集读数，我们可以实行简单的贝叶斯计算，以根据分布中插入大小L1和L2的概率的比率将读数概率性地分配给两个异构体。

这是对当前技术的显著改进，其该情况下总体上检视读数为携载用于两个异构体的等量证据权重。因此，以所描述的方式合并插入大小是不忽略相关证明的方式。

因此，在这些不明确的情况下，我们可以计数异构体，如下：

(1)计算插入大小的概率分布。(由P(s)指代具有插入大小s的读数的概率。)应注意这仅需要计算一次，并且这可以在预处理步骤中完成。

(2)如果两个双端配偶中的每个映射在外显子内部：

(2a)制作包含那些外显子的所有异构体的列表。

(2b)将每个异构体与将由已经从该异构体导出的读数暗示的插入大小相关联。

(2c)为(2b)中所计算的插入大小中的每个计算P(s)。

(2d)将每个异构体与P(s)除以(2c)中所计算的所有P(s)值的总和的相关值相关联。

应该注意的是，本文所描述的技术将以不影响使用“下游”工具(也就是说，将读数计数作为输入且给出其它有用的输出的工具)的能力的方式改进异构体的读数计数。

9.构造“最大DAG”。

在一些实施例中，构造DAG可以是有益的，其中，在此DAG中每对节点用从表示“较早”或“左侧”外显子到“较迟”或“右侧”外显子的节点的边缘连接。

使用Gencode或类似的经注释的转录组，创建节点以表示外显子和其它特征。Gencode提供足够信息以将这些节点与此信息相关联作为节点的序列、其在参考序列中占据的位置(并且因此其长度)等。因此，使用包括耦接到非暂时性存储器的处理器的计算机系统，来自基因组的多个特征中的每个被表示为节点。对于该多个特征中的每对，创建边缘，边缘通过该对的两个成员的近端连接该两个成员。此移动的生物学解释是猜测表示遵守外显子出现在基因组中的顺序的那些外显子的每个组合的异构体的存在，或者至少促进发现该异构体。对遵守外显子出现在基因组中的顺序敏感既有感应性的原因又有生物化学的原因。

例如，报告表明在多外显子前mRNA中的外显子总是按顺序保持。参见例如Black，2005，《剪接难题的简单答案(Asimpleanswerforasplicingconundrum)》，《美国国家科学院院刊(PNAS)》102：4927-8。没有被任何特定的机制束缚，从理论上可以讲，剪接机械设备辨识内含子的开始和结束，以及分支点腺嘌呤或接近剪接接合点的其它保守序列。剪接机械设备或剪接体沿着前mRNA从5'到3'移动，移除内含子且剪接外显子。因此看来存在用于以5'到3'顺序连接所有外显子的现象学以及文献支撑，在5'到3'顺序中，所有的外显子看起来沿着基因组DNA链。还参见Shao等人，2006，《人类中外显子重复接收、外显子加扰和反式剪接的生物信息学分析(Bioinformaticanalysisofexonrepetition,exonscramblingandtrans-splicinginhumans)》，《生物信息学(Bioinformatics)》22：692-698。对于某些假设，本发明的方法遵循人类前mRNA被剪接成保留由基因组序列定义的外显子顺序的线性分子的假设。

另外，在期望在有向非循环数据结构中包含外显子加扰的实例(参见例如Carrington等人，1985，《在伴刀豆球蛋白A的转译后修饰期间发生多肽接合(Polypeptideligationoccursduringpost-translationalmodificationofconcanavalinA)》，《自然(Nature)》313：64-67)的情况下，可以通过包含双重复外显子作为“虚构”大体上保持外显子以它们的原始顺序剪接的假设(或者，在可替代的方案中，可以废除该假设)。

因此，在例如基因组包含处于以下顺序：ABCDEFGHI的以下外显子的情况下

使用短划线作为边缘且字母作为节点，完整的该组边缘可以被表示为：

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。

果这足以证明发现所有的异构体而不是一种，并且该一种异构体——给出示例：ADEHG

则，借助仿真的或虚构的G——称它为G'，该非循环数据结构可以是“钉齿状的”。然后，可以由A-D-E-H-G'表示该异构体，并且结构不需要包含从H到G的任何边缘。然而，可以发现可以忽略外显子加扰和由所有包含的边缘保持的基因组外显子顺序。

可以通过使用来自经注释的转录组的输入数据，首先制备本发明的DAG。同样地，从头开始，DAG将包含所有已知的异构体。如果制成最大的DAG，则DAG将包含保持基因组外显子顺序的所有可能的边缘，并且然后，将因此包含既表示已知的异构体又表示新颖的异构体的节点和边缘。在存储器中节点和边缘被存储为有向非循环数据结构，并且方法包含将双端序列读数与包括由边缘的子集连接的特征的子集的路径进行比对，从而识别对应于路径的异构体。

当将新的数据与DAG(例如，双端RNA测序读数)比对时，那些数据可以包含在经注释的参考中未被识别的外显子。因为在该示例中外显子是新颖的，所以该外显子是异构体。使用DAG作为参考理想地适合于该情况，因为将新颖的异构体与现有的DAG进行比对还可以引导在DAG中创建新的节点的以表示新颖的外显子。因此，在所识别的异构体是新颖的情况下，比对转录序列可以包括在有向非循环数据结构中创建新节点。还可以为剪接变异体创建新节点，例如，其中一个外显子通过单个密码子不同于另一个外显子。

优选地，在分析双端读数之前且独立于分析双端读数，初始地全部或部分地构建最大DAG。因此，方法可以包含在比对转录本之前，为多个特征中的每对创建边缘。出于与有向非循环数据结构的本质相关的原因，数据结构中表示的参考可以任意大的，例如，数百或数千外显子、内含子或两者。

10.本发明的系统

图10例示用于实施本文所描述的方法的计算机系统401。本发明的系统可以包含图10中所示的任一个部件或任何数量的部件。一般来说，系统401可以包含能够通过网络415彼此通信的计算机433和服务器计算机409。另外，可以可选地从数据库405(例如，本地或远程)获得数据。在一些实施例中，系统包含用于获得测序数据的仪器455，仪器455可以被耦接到测序仪计算机451用于序列读数的初始处理。

在一些实施例中，通过并行处理实行方法，并且服务器409包含具有并行架构的多个处理器，即，能够将多个序列(例如，RNA测序读数)与参考序列构造(例如，DAG)进行比较的处理器和存储器的分布式网络。系统包括多个处理器，这些处理器同时执行多个读数与参考序列构造之间的多个比较。虽然可能有其它混合配置，但是并行计算机中的主存储器通常或者在单个地址空间中的所有处理元件之间共享，或者呈分布式，即，每个处理元件具有其自身的本地地址空间。(分布式存储器是指这样的事实：存储器以逻辑方式分布，但通常暗示它也以物理方式分布。)分布式共享存储器和存储器虚拟化组合这两种方法，其中处理元件具有其自身的本地存储器以及对非本地处理器上的存储器的存取权。对本地存储器的存取通常比对非本地存储器的存取更快。

其中可以用相等时延和带宽存取主存储器的每个元件的计算机架构被称为均匀存储器存取(UMA)系统。通常，UMA系统只能通过共享的存储器系统来实现，其中存储器并非以物理方式分布。不具有该特性的系统被称为不均匀存储器存取(NUMA)架构。分布式存储器系统具有不均匀存储器存取。

可以以若干种方式在硬件中实施处理器-处理器和处理器-存储器通信，这些方式包含经由共享的(或者多端口的或者多路复用的)存储器、纵横开关、共享的总线或无数拓扑的互连网络(包含星形、环形、树形、超立方体、胖超立方体(在一个节点处具有不止一个处理器的超立方体))或n维网格。

基于互连网络的并行计算机必须合并路由以实现在并非直接连接的节点之间的消息传递。用于处理器之间的通信的媒体很可能在大型多处理器机器中分层。此类资源在市面上可购买用于专用用途，或可以经由例如亚马逊云计算的“云”存取这些资源。

计算机总体上包含通过总线耦接到存储器和输入-输出(I/O)机构的处理器，存储器可以包含RAM或ROM，并且优选地包含存储指令的至少一个有形的非暂时性媒体，该指令可执行以致使系统实行本文所描述的功能。本领域中的技术人员将认识到，根据需要或者最适合于实行本发明的方法，本发明的系统包含一个或多个处理器(例如，中央处理单元(CPU)、图形处理单元(GPU)等)、计算机可读存储设备(例如，主存储器、静态存储器等)，或其组合，它们通过总线彼此通信。

处理器可以是本领域中已知的任何合适的处理器，诸如由英特尔(加利福尼亚州圣克拉拉(SantaClara,CA))以商标XEONE7出售的处理器，或由AMD(加利福尼亚州桑尼维尔(Sunnyvale,CA))以商标OPTERON6200出售的处理器。

根据本发明的输入/输出设备可以包含视频显示单元(例如，液晶显示器(LCD)或阴极射线管(CRT)监视器)、字母数字输入设备(例如，键盘)、光标控制设备(例如，鼠标或触控板)、磁盘驱动器单元、信号生成设备(例如，扬声器)、触摸屏、加速计、麦克风、蜂窝无线电频率天线，以及网络接口设备，网络接口设备可以是例如网络接口卡(NIC)、Wi-Fi卡或蜂窝调制解调器。

以引用的方式并入

贯穿本公开已经参考并且引用了其它文档，诸如专利、专利申请、专利公开案、期刊、书籍、论文、网络内容。所有此类文档在此以全文引用的方式并入本文中用于所有目的。

等效物

除本文示出且描述的那些之外，本领域的技术人员将从本文档的完整内容(包含对在本文中引用的科学和专利文献的参考)中明白本发明的各种改良及本发明的许多进一步的实施例。本文中的主题含有重要信息、范例和指南，它们可适于本发明在本发明的各种实施例及其等效物中的实践。

Claims

1.一种分析转录组的方法，所述方法包括：

从转录组获得一对双端读数；

找出在所述对的第一读数和有向非循环数据结构中的节点之间的具有最佳评分的比对，所述数据结构包括表示RNA序列的节点和连接多对所述节点的边缘；

识别包含所述节点的候选路径，所述节点通过具有与所述一对双端读数的插入长度大体上类似的长度的路径连接到下游节点；以及

将所述双端读数与所述候选路径比对以确定最佳评分比对。

2.根据权利要求1所述的方法，进一步包括从所述转录组获得多个双端读数。

3.根据权利要求2所述的方法，进一步包括为所述多个双端读数对中的每个确定最佳评分比对。

4.根据权利要求3所述的方法，其中所述多个双端读数是从来自所述转录组的cDNA片段获得的。

5.根据权利要求1所述的方法，进一步包括基于所确定的最佳评分比对来识别异构体。

6.根据权利要求5所述的方法，其中所述所识别的异构体是新颖的异构体。

7.根据权利要求6所述的方法，进一步包括更新所述有向非循环数据结构以表示所述新颖的异构体。

8.根据权利要求7所述的方法，其中更新所述有向非循环数据结构以表示所述新颖的异构体包括添加至少一个新节点。

9.根据权利要求8所述的方法，其中所述节点和所述下游节点表示一对外显子，且所述一对双端读数是从所述一对外显子获得的。

10.根据权利要求1所述的方法，进一步包括从包括所述一对双端读数的任何比对计算中排除不是候选路径的任何路径。

11.根据权利要求1所述的方法，进一步包括：

将多个双端读数对与所述有向非循环数据结构进行比对；

确定所述多个双端读数对的所述所比对的对之间的距离和所述所比对的对之间的所述距离的频率；以及

确定一组异构体路径和异构体频率，以使得：表示以所述异构体频率穿过所述结构的所述异构体路径引起多对特征以所述所比对的对之间的所述距离的所述频率被包含在所述异构体路径中。

12.根据权利要求11所述的方法，其中所述特征包含外显子。

13.根据权利要求12所述的方法，其中所确定的所述组异构体路径中的每个路径表示存在于生物体中的转录异构体。

14.根据权利要求11所述的方法，其中所述所比对的对之间的所述距离包括从每对的第一比对成员的上游端到所述对的第二比对成员的下游端的路径长度。

15.一种用于转录组学的系统，所述系统包括：

包括耦接到具有存储于其中的有向非循环数据结构的非暂时性存储器的处理器的计算机系统，所述数据结构包括表示RNA序列的节点和连接多对所述节点的边缘，其中所述系统可操作以：

从cDNA片段获得一对双端读数；

找出在所述对的第一读数和所述数据结构中的节点之间的具有最佳评分的比对；

识别包含所述节点的候选路径，所述节点通过具有与所述双端读数的插入长度大体上类似的长度的路径连接到下游节点；以及

将所述双端读数与所述候选路径比对以确定最佳评分比对。

16.根据权利要求15所述的系统，进一步可操作以从cDNA片段获得多个双端读数。

17.根据权利要求16所述的系统，进一步可操作以基于所述所确定的最佳评分比对来识别新颖的异构体。

18.根据权利要求17所述的系统，进一步可操作以更新所述有向非循环数据结构以表示所述新颖的异构体。

19.根据权利要求18所述的系统，其中更新所述有向非循环数据结构以表示所述新颖的异构体包括添加至少一个新节点。

20.根据权利要求15所述的系统，进一步可操作以从包括所述对双端读数的任何比对计算中排除不是候选路径的任何路径。

21.一种推断转录组中异构体频率的方法，所述方法包括：

将来自cDNA的多个双端读数对与参考数据结构进行比对，所述参考数据结构包括来自基因组的每个被表示为节点的多个特征，其中多对的所述多个特征在它们的近端处通过边缘连接；

确定所述双端读数的所述所比对的对之间的距离和所述所比对的对之间的所述距离的频率；以及

22.根据权利要求21所述的方法，其中所述cDNA大体上表示生物体的整个转录组，并且所述特征包含外显子。

23.根据权利要求22所述的方法，其中所确定的所述组异构体路径中的每个路径表示存在于所述生物体中的转录异构体。

24.根据权利要求23所述的方法，其中确定所述组异构体路径包含发现新颖的异构体路径。

25.根据权利要求24所述的方法，进一步包括更新所述参考数据结构以包含所述新颖的异构体路径。

26.根据权利要求21所述的方法，其中所述所比对的对之间的所述距离包括从每对的第一比对成员的上游端到所述对的第二比对成员的下游端的路径长度。

27.根据权利要求21所述的方法，其中确定所述组异构体路径包括解算一组线性方程式。

28.根据权利要求27所述的方法，其中确定所述组异构体路径包括在解算所述组线性方程式之前从生物学数据获得所述异构体频率中的一个。

29.根据权利要求21所述的方法，进一步包括在确定一组异构体路径和异构体频率之前使在所述所比对的对之间的所述距离的所述频率标准化。

30.根据权利要求21所述的方法，其中所述数据结构包括表示基因特征数据库中的每个经注释的特征的节点。

31.一种推断转录组中异构体频率的系统，所述系统包括：

包括耦接到具有存储于其中的参考数据结构的非暂时性存储器的处理器的计算机系统，所述参考数据结构包括来自基因组的每个被表示为节点的多个特征，其中多对的所述多个特征在它们的近端处通过边缘连接，其中所述系统可操作以：

确定来自cDNA的多个双端读数对之间的距离和所述对之间的所述距离的频率；以及

32.根据权利要求31所述的系统，其中所述cDNA大体上表示生物体的整个转录组，所述特征包含外显子，并且所述系统进一步可操作以将来自cDNA的多个双端读数对与所述参考数据结构进行比对。

33.根据权利要求32所述的系统，其中所确定的所述组异构体路径中的每个路径表示存在于所述生物体中的转录异构体。

34.根据权利要求33所述的系统，其中确定所述组异构体路径包含发现新颖的异构体路径。

35.根据权利要求34所述的系统，进一步可操作以更新所述参考数据结构以包含所述新颖的异构体路径。

36.根据权利要求31所述的系统，其中所述所比对的对之间的所述距离包括从每对的第一比对成员的上游端到所述对的第二比对成员的下游端的路径长度。

37.根据权利要求31所述的系统，其中确定所述组异构体路径包括解算一组线性方程式。

38.根据权利要求37所述的系统，其中确定所述组异构体路径包括在解算所述组线性方程式之前从生物学数据获得所述异构体频率中的一个。

39.根据权利要求31所述的系统，进一步可操作以在确定所述组异构体路径和异构体频率之前使在所述所比对的对之间的所述距离的所述频率标准化。

40.根据权利要求31所述的系统，其中所述数据结构包括表示远程基因特征数据库中的每个经注释的特征的节点。