用于在有向非循环结构中使用双端数据的系统和方法
【专利摘要】本申请涉及一种分析转录组的方法，该方法包括从来自生物体的转录组获得至少一对双端读数，找出在该对的第一读数和有向非循环数据结构(该数据结构具有表示RNA序列诸如外显子或转录本的节点和连接多对节点的边缘)中的节点之间的具有最佳评分的比对，识别包含节点的候选路径，该节点通过具有与该对双端读数的插入长度大体上类似的长度的路径连接到下游节点，并且将双端读数与候选路径比对以确定最佳评分比对。
【专利说明】
用于在有向非循环结构中使用双端数据的系统和方法
[0001] 相关申请的交叉参考
[0002] 本申请要求2014年1月17日提交的美国专利申请序列号14/157,979的权益和优先 权，该美国专利申请要求2013年10月21日提交的美国临时专利申请序列号61/893，467的权 益和优先权，该临时专利申请的内容W引用方式并入本文。
[000；3] 序列表
[0004] 本申请包含已通过网络电子提交系统化FS-Web)WASCII格式提交的序列表，并且 W其全部内容W引用方式并入本文。2014年10月14日创建的ASCII格式化的序列表被命名 为 SBG-009-01W0-2-Sequences_ST25,并且大小为 1,217 字节。
技术领域
[0005] 本发明大体上设及生物信息学，并且尤其设及使用有向非循环数据结构进行分 析。
【背景技术】
[0006] 对于特定的发育阶段或生理状态，转录组是细胞中完整的一组转录本和它们的数 量。转录组对生物体的功能、发育和疾病极为重要。生物体的本质和活力来源于所有的转录 本(包含mRNA、非编码RNA和小RNA)、它们的表达水平、基因剪接模式和转录后修饰。事实上， 人类基因组中与期望相比大得多的部分现在已知为被转录的。参见Bedone等人，2004,《具 有基因组嵌合阵列的人类转录序列的全局识别（Global identification of human transcribed sequences with genome tiling arrays)》，《科学（Science)》306:2242- 2246 〇
[0007] 已经报道基于下一代测序(NGS)技术的用于转录组分析的方法。参见，例如，Wang 等人，2009，《RNA测序：用于转录组学的革命性的工具（RNA-Seq: a revolutionary tool for transcriptomics)》，《遗传学自然评论（化t Rev Genet)》10(l):57-63 2009eRNA测序 (将RNA测序)方法承诺快速生成高容量的转录组数据。
[000引然而，RNA测序面临信息挑战，诸如存储和分析大量数据，运必须被解决W更好地 使用读数。比对和分析RNA测序读数不仅在设及的信息的本质上呈现出问题，还在容量上呈 现出问题。具有大数据集的读映射是计算上的挑战，并且用于差异表达的分析方法才开始 出现。参见Garber等人，2011，《使用RNA测序的用于转录组注释和量化的计算方法 (Comput曰tion曰1 methods for tr曰nscriptome 曰nnot曰tion 曰nd qu曰ntific曰tion using RNA-Seq)》，《自然方法(化t Meth)》8(6):469-477。映射来自RNA测序的大量读数的计算需 求被将来可获得的大量参考数据大大倍增。
[0009]例如，GENC0DE联合体力图识别已经被报道的人类基因组中的所有基因特征。 Harrow等人，2012,《GENC0DE:用于编码计划的参考人类基因组注释（GENC0DE: The reference human genome annotation for The ENCODE Project)》，《基因组石开究（Genome Res)》22:1760-1774。由GENCODE计划包含的材料的量是庞大的。不仅仅要继续添加新的蛋 白质编码基因位点，而且经注释的可变剪接转录本的数量稳定地增加 。GENCODE 7版本包含 超过20,000个蛋白质-编码和几乎10,000个长的非编码RNA基因位点（IncRNA)，W及未在其 它源中展现的超过33，000个编码转录本。GENCODE还包含其它特征诸如非转录区域(UTR)、 长基因间非编码RNA(lincRNA)基因、短非编码RNA和可变剪接模式。即使RNA测序研究开始 仅具有有限量的新数据，来自源诸如GENCODE的潜在参考数据的容量也是如此大，使得理解 新数据在计算上是具有挑战性的。
【发明内容】
[0010] 本发明大体上提供了用于使用组织为有向非循环图(DAG)或类似的数据结构的参 考对双端读数进行分析的系统和方法。参考中的特征诸如外显子和内含子在DAG中提供节 点，并且那些特征通过边缘W它们的规范的基因顺序彼此链接。DAG可W按比例绘制为任何 大小，并且事实上，在第一实例中可W被填入来自外来经注释的参考的导入序列。在新数据 包含双端读数的情况下，一对中的每个读数可W潜在地沿长度L的基因组且在复杂的DAG中 与任何位置比对，可存在多得多的位置（即，平均来说，如果DAG是横越任何基因位点η重的， 则比对一对双端读数相当于从(nL)~2个选项中进行选择）。然而，插入长度信息提供了用于 将潜在的比对约束在DAG内的手段一一不仅是在路径内远离的距离中，还通过排除路径、将 选项的数量限制为相对于L是线性的。双端读数还可W提供其它信息。例如，双端读数可W 提供用于修饰DAG的基础，可W被用于推断异构体频率、或者其它。事实上，具有插入长度信 息的双端读数集携载信息的品质，该品质非常适合于聚焦基于DAG的比对和更新DAGW表示 可能没有被发现或注意到的实际上出现的固有异构体。因此，双端读数特别适合于RNA测序 计划，允许相对于很大的参考集分析很大的RNA测序数据集。使用双端RNA测序读数发现转 录组的所有的异构体并不是计算上难处理的，并且转录组学可先前不可能的步伐和吞 吐量移动。可W从一个或任何数量的生物体研究一个或任何数量的转录组W增加对生物体 功能、发育和疾病的理解。
[0011] 本发明的某方面设及创建包含参考的特征(例如，外显子、内含子等）的DAG，参考 的特征W基因组顺序通过它们的近端边缘连接。该DAG数据结构具有转录组学中的应用，包 含例如识别异构体。此DAG参考被用于设及包含配对读数的RNA测序数据的转录组学中。可 W使用包含作为节点的特征(例如从经注释的源基因组或转录组创建的）的DAG参考来分析 RNA测序数据。本发明的系统和方法提供了使用双端读数的改进的比对、使用双端读数信息 构建DAG、推断异构体概率或频率，或其组合。
[0012] 在某些方面中，本发明提供了分析转录组的方法，该方法包括从来自生物体的转 录组获得至少一对双端读数，找出在该对的第一读数和有向非循环数据结构中的节点之间 的具有最佳评分的比对(该数据结构具有表示RNA序列诸如外显子或转录本的节点和连接 多对节点的边缘），识别包含节点的候选路径，该节点通过具有与该对双端读数的插入长度 大体上类似的长度的路径连接到下游节点，W及将双端读数与候选路径比对W确定最佳评 分比对。方法可W包含从包括该对双端读数的任何比对计算中排除不是候选路径的任何路 径。
[0013] 该方法适合于任何数量的双端读数。节点和下游节点可W表示一对外显子，该对 双端读数从该对外显子获得。该方法可W包含基于所确定的最佳评分比对来识别异构体。 在一些实施例中，所识别的异构体是新颖的异构体。有向非循环数据结构可w被更新w表 示新颖的异构体(例如，通过添加至少一个新节点或边缘）。
[0014] 本发明的相关方面提供用于转录组学的计算机系统。该系统具有禪接到非暂时性 存储器的处理器，并且被用于储存有向非循环数据结构，该有向非循环数据结构具有表示 RNA序列的节点和连接多对节点的边缘。系统可操作W从转录组获得一对双端读数，在该对 的第一读数和数据结构中的节点之间找出具有最佳评分的比对，识别包含节点的候选路 径，该节点通过具有与该对双端读数的插入长度大体上类似的长度的路径连接到下游节 点，W及将该双端读数与该候选路径比对W确定最佳评分比对。优选地，系统从随后的比对 尝试中排除所有非候选路径。在某些实施例中，该系统进一步可操作W从cDNA片段中获得 多个双端读数。优选地，系统可操作W为多个双端读数对中的每个确定最佳评分比对。cDNA 片段可W表示单个转录组。系统可W基于所确定的最佳评分比对来识别异构体，诸如新颖 的异构体。系统可操作W更新有向非循环数据结构W表示新颖的异构体。在一些实施例中， 更新有向非循环数据结构W表示新颖的异构体包括添加至少一个新节点。节点和下游节点 可W表示一对外显子，双端读数对从该外显子对获得。系统可操作W从包括该对双端读数 的任何比对计算中排除不是候选路径的任何路径。
[0015] 本发明的其它方面通过例如可用于推断异构体频率的系统和方法提供了表达分 析。
[0016] 在某些方面中，本发明提供推断转录组中异构体频率的方法。该方法包含从生物 体的转录组获得多个双端读数对，并且可选地将那些读数与参考数据结构进行比对。数据 结构包含来自基因组的每个被表示为节点的多个特征。多对节点或特征在它们的近端处通 过边缘连接（即，在对内，两个对成员在基因组上在它们的近端处通过边缘彼此连接；由于 任何对成员还可W是其它对的成员，所W特征的任一端或两端可W通过一个或任何数量的 边缘被连接到其它特征）。确定双端读数的所比对的对之间的距离和该所比对的对之间的 距离的频率。确定一组异构体路径和异构体频率，W使得表示W异构体频率穿过结构的异 构体路径引起多对特征W所比对的对之间的距离的频率被包含在该异构体路径中。所确定 的该组异构体路径中的每个路径可W表示存在于生物体中的转录异构体。所比对的对之间 的距离可W是从每对的第一比对成员的上游端到该对的第二比对成员的下游端的路径长 度。在某些实施例中，确定该组异构体路径包含发现新颖的异构体路径。可W更新参考数据 结构W包含新颖的异构体路径。
[0017] 可W用代数方法例如通过解算一组线性方程式实行确定该组异构体路径。在一些 实施例中，在解算该组线性方程式(例如，所W可W确定线性方程式中的其余的变量)之前， 从生物学数据获得异构体频率中的一个。另外地或可替代地，在确定一组异构体路径和异 构体频率之前，可W使所比对的对之间的距离的频率标准化。
[0018] 可W从现有的基因特征数据库构建数据结构。例如，每个节点可W表示来自此数 据库的特征。
[0019] 在相关的方面中，本发明提供推断转录组中异构体频率的计算机系统。系统包含 禪接到非暂时性存储器的处理器，并且存储参考数据结构，该参考数据结构包括来自基因 组的每个被表示为节点的多个特征，其中该多个特征的对在它们的近端处通过边缘连接。 系统可W被用于获得或确定来自转录组的多个双端读数对之间的距离和该对之间的距离 的频率。然后可w确定一组异构体路径和异构体频率，w使得表示w异构体频率穿过结构 的异构体路径引起多对特征W所比对的对之间的距离的频率被包含在该异构体路径中。特 征可W包含外显子，并且系统进一步可操作W将多个双端读数对与参考数据结构进行比 对。所确定的该组异构体路径中的每个路径可W表示存在于生物体中的转录异构体。确定 该组异构体路径可W包含发现新颖的异构体路径。系统可W更新参考数据结构W包含新颖 的异构体路径。所比对的对之间的距离可W包括从每对的第一比对成员的上游端到该对的 第二比对成员的下游端的路径长度。在一些实施例中，确定该组异构体路径包括解算一组 线性方程式。在某些实施例中，确定该组异构体路径包括在解算该组线性方程式之前，从生 物学数据获得异构体频率中的一个。系统可操作W在确定该组异构体路径和异构体频率之 前，使在所比对的对之间的距离的频率标准化。数据结构可W包含表示远程基因特征数据 库中的每个经注释的特征的节点。
【附图说明】
[0020] 图1示出与DAG比较的读数的图形表示。
[0021] 图2示出对应于比较的实际矩阵。
[0022] 图3提供了本发明的方法的图式。
[0023] 图4图示了假设的DAG和双端读数对。
[0024] 图5示出穿过假设的DAG的候选路径。
[0025] 图6示出来自假设的DAG的候选路径。
[00%]图7呈现包括序列的Ξ个可替代的历史记录。
[0027]图8根据本发明描绘了比对。
[002引图9示出异构体频率分布。
[0029] 图10例示了用于实施本发明的方法的系统。
【具体实施方式】
[0030] 应当理解，一般来说，基因不是用于代码单个、特定蛋白质。代替地，在不同的时间 W不同的比例转录不同子组的给定的基因的外显子。由运些不同的转录本创建的不同的蛋 白质被称为异构体，并且相关地，不同的转录本本身也被称为异构体。如上面所论述的，现 代生物信息学已经提供了潜在地用作参考的至少部分地经注释的转录组。例如，Gencode计 划已经不仅编录基因，还编录人类基因组中的外显子(连同其它区域）。
[0031] 本文所描述的系统和方法被用于有效地利用那些经注释的转录组中的信息。实际 上，本发明包含将双端读数与经注释的转录组进行比较理想地适合于分析方法的桐见，在 该分析方法中序列数据被存储为基于DAG的数据结构。
[0032] 本发明的一些实施例提供用经注释的转录组中的数据创建有向非循环数据结构 的方法。在此数据结构中，来自经注释的转录组的特征(例如，外显子和内含子)被表示为节 点，该节点通过边缘连接。因此，由穿过基于DAG的数据结构的对应路径表示来自现实世界 转录组的异构体。
[0033] 本发明的某方面设及创建包含来自一个或多个已知参考的特征诸如内含子和外 显子的DAG。在本领域中DAG被理解为是指可W被呈现为图表的数据W及呈现数据的图表。 本发明提供用于将DAG存储为可W由计算机系统读取W用于生物信息学处理或用于呈现为 图表的数据的方法。可任何合适的格式保存DAG，该合适格式包含例如节点和边缘的列 表、矩阵或表示矩阵的表、一批阵列或表示矩阵的类似的变量结构，W内置有图表语法的语 言保存，W用于图表目的的通用标记语言保存，或其它合适格式。
[0034] 在一些实施例中，DAG被存储为节点和边缘的列表。此方式的其中一种是创建包含 所有节点和所有边缘的文本文件，该所有节点具有分配给每个节点的ID，该所有边缘中的 每个具有开始和结束节点的节点ID。因此，例如，如果为两个语句"See化ne run"和"Run， Jane run"创建DAG,则可W创建不分大小写的文本文件。可W使用任何合适的格式。例如， 文本文件可包含逗号分隔值。将该DAG命名为"Jane"用于将来的参考，W该格式，DAG"Jane" 可W读作如下：Isee、化1111、3^'日]16、41'11]1，1-3、2-3、3-4。本领域的技术人员将理解，该结构可 适用于图7和图8中所表示的内含子和外显子，W及下面的随附论述。
[0035] 在某些实施例中，DAG被存储为表示矩阵(或一批阵列或表示矩阵的类似的变量结 构）的表，其中NX N矩阵中的（i，j)项指示节点i和节点j被连接(其中N是W基因组顺序包含 节点的矢量）。为了 DAG是非循环的，仅需要所有的非零项在对角线上面(假定节点W基因组 顺序表示）。在二进制情况下，0项表示从节点i到节点j不存在边缘，并且1项表示从i至Ijj有 一个边缘。本领域的技术人员将理解矩阵结构允许除0到1之外的值与边缘相关联。例如，任 何项可W是指示权重或所使用的次数，反映世界中所观测的数据的一些固有品质的数值。 矩阵可W被写入到文本文件作为表或一系列线性的行（例如，首先是行1，接着是分隔符 等），因此提供简单的串行化结构。
[0036] 在定义节点之后，使矩阵DAG串行化的一种有用的方式将是为项使用逗号分隔值。 使用该格式，DAG "Jane""将包含与上述相同的节点定义，接着是矩阵项。该格式可W读作：
[0037] lseej;run、3jane、4;run
[0038] ，，1，\，，1，\，，，1\，，，
[0039] 其中简单地省略了零(0)的项，并且7"是换行符。
[0040] 本发明的实施例包含W内置有图表语法的语言存储DAG。例如，具有被称为图表可 视化(Gra地viz)的图表可视化软件包的DOT语言提供了可W被用于存储具有辅助信息且可 W使用从Graphviz网站商购获得的多个工具转化成图表文件格式的数据结构。Gra地viz是 开源图表可视化软件。图表可视化是一种将结构信息表示为抽象的图表和网络的图式的方 式。其具有在联网、生物信息学、软件工程、数据库和网页设计、机器学习中的应用，W及在 用于其它技术领域的可视界面中的应用。Graphviz布局规划采取W简单文本语言为形式的 图表描述，并且W有用的格式(诸如用于网页的图像和SVG;用于包含在其它文档中的PDF或 附言;或者在交互式图表浏览器中的显示)制作图式。
[0041] 在相关的实施例中，DAGW用于图表目的的通用标记语言或其它格式被存储。根据 上面线性文本文件或逗号分隔矩阵的描述，本领域的技术人员将认识到，语言诸如XML可W 被用于(扩展到)创建定义节点和它们的报头或ID、边缘、权重等的标签(标记）。然而，DAG被 结构化且被存储，本发明的实施例包括使用节点表示特征诸如外显子和内含子。运为分析 RNA测序读数与发现、识别和表示异构体提供了有用的工具。
[0042] 如果节点表示特征或特征的片段，则可W由穿过运些片段的路径表示异构体。由 将一些先前的外显子连接到一些之后的外显子的边缘表示外显子是"被跳过的"。本文所呈 现的是用于构造 DAGW表示基因的可变剪接或异构体的技术。在Lee和Wang, 2005,《可变剪 接的生物信息学分析（Bioinformatics analysis of alternative splicing)》，《简明生 物信息学（Brief Bioinf)》6( 1):23-33 ;Heber等人，2002,《剪接图表和EST装配问题 (Splicing graphs and EST assembly problems)》，《生物信息学（Bioinformatics)》 18Suppl:sl81-188;Leipzig等人，2004,《可变剪接巷道(ASG):桥接基因组和转录组之间的 空j立（The alternative splicing 邑allery(ASG):Bridging the gap between genome and transcriptome)》，《核酸研究（Nucl Ac 1?63)》23(13):3977-2983;^及1^66日1111和 Dewey ,2013,《来自具有概率性的剪接图表的RNA测序数据的可变剪接的推断（Inference of 曰Iterative splicing from RNA-Seq d曰t曰 with probabilistic splice graphs)》， 《生物信息学(Bioinformatics)》29(18) :2300-2310中论述了可变剪接，运些文献中的每个 文献的内容W引用方式并入。可W在Florea等人，2005,《具有AIR的基因和可变剪接注释 (Gene and alternative splicing annotation with AIR)》，《基因组石开究（Genome Research)》15:54-66 ;Kim等人，2005，《ECgene:用于可变剪接的基于基因组的EST群集和基 因建模（EC邑ene :Genome-based EST clustering and 邑ene modeling for alternative splicing)》，《基因组研究（Genome Research)》15:566-576; W及Xing等人，2006,《用于来 自剪接图表的全长异构体的概率性重构的期望-最大化算法（An expectation- maximization algorithm for probabilistic reconstructions of full-length isoforms from splice gra地s)》，《核酸研究（Nucleic Acids Resea;rch)》，34,3150-3160 中找出另外的论述，运些文献中的每个文献的内容W引用方式并入。
[0043] 1.使用双端读数
[0044] 本发明设及包括RNA测序数据的转录组学，该RNA测序数据包含双端读数。本发明 的系统和方法提供使用双端读数的改进的比对。在某些实施例中，双端读数被用于约束比 对，由此大大降低RNA测序计划的计算成本。
[0045] 很多测序数据是双端数据。存在其中此类数据的特征应该在比对中体现他们自身 的方式:双端配偶在比对中被分配的位置之间的距离应该大致对应于由测序运行给定的插 入大小。一般来说，由插入大小给定的信息设及比对算法试图掲示的信息。
[0046] 用于产生成对的核酸读数的若干方法在本领域中是已知的。通常，该方法包括选 择性地将核酸样本片段化从而得到已知大小分布的序列的一些形式。通常用官能团或通过 自环化来保护片段化的样本的端部，于是剩余的核酸样本材料被移除。然后，端部被扩增且 测序，导致被认为彼此之间产生一些距离的核酸读数的群体。
[0047] 配偶对样本制备成套装置是例如从化xtera Illumina(加利福尼亚州圣迭戈(San Diego, CA))商购的。为了制备成对的核酸用于在此平台上测序，核酸样本被片段化成长度 在化b和化b之间的段。然后，使用亲和柱分隔，片段端被生物素标记W促进之后的恢复。生 物素标记的端部被接合W创建环状DNA段。此时，来自样本的非环化的核酸被消化且随后被 移除。然后，剩余的环化的DNA被再片段化W产生适合于群集和测序的片段。可W用亲和柱 或链霉亲和素包被的小珠从环化的DNA的剩余片段移除被生物素标记的原始端。然后，恢复 的端进行端修复、聚腺巧酸加尾，并且在扩增和测序之前添加测序衔接子。
[004引 2.比对
[0049]本发明提供用于将序列读数与表示经注释的参考的有向非循环数据结构进行比 对的方法和系统。使用本发明的比对算法，即使存在与RNA测序结果相关联的大数量和大体 上全面的参考，读数也可W被快速地映射。与在其它情况下相同，DAG是比对的理想目标。与 之后在某个比对的工作流程中使用经注释的转录组相比，实行与由经注释的转录组构造的 DAG进行比对是有效得多的。
[0050] 通过使用DAG作为参考获得大量益处。例如，和与一个参考进行比对，并且然后试 图根据经注释的转录组调整某个比对的结果相比，与DAG进行比对更精确。运主要是因为前 一方法在用于初始比对的序列和转录组中所表示的另一序列之间强加了不自然的不对称。 和试图与每种物理可能性(此类可能性的数量将在数个接合点中总体上指数增长）的线性 序列进行比对相比，与潜在地表示所有的相关物理可能性的目标进行比对是计算上有效得 多的。
[0051] 本发明的实施例包含将一个或多个读数与参考数据结构进行比对。该比对可W包 含基于DAG的比对，其找出最优地表示读数的数据结构的边缘和节点。该比对还可W包含W 成对的方式将读数与来自DAG的节点的外显子比对的实例。
[0052] 成对比对总体上包括沿目标的一部分放置一个序列，根据算法引入空位，对两个 序列匹配的程度评分，和优选地沿着参考对不同的位置重复进行该成对比对。最优评分匹 配被认为是比对，并且表示关于序列数据表示的内容的推断。在一些实施例中，对一对核酸 序列的比对进行评分包括为置换和插入缺失的概率设置值。当比对了单个碱基时，匹配或 不匹配W置换概率贡献于比对评分，置换概率可W是例如对于匹配为及对于不匹配为- 0.33。插入缺失W将空位罚分从比对评分中扣除，空位罚分可W是例如-1。空位罚分和置换 概率可W基于关于序列如何演化的经验知识或先验假设。它们的值影响得到的比对。特别 地，空位罚分和置换概率之间的关系影响置换或插入缺失是否在得到的比对中将是有利 的。
[0053] 正式地说，比对表示两个序列X和y之间的推断关系。例如，在一些实施例中，序列X 和y的比对A将X和y分别映射到可W包含空隙的另外两个字符串X'和y'，使得：（i)lx'l = ly'l;(ii)从x'和y'移除空隙应该分别回到x和y;W及(III)对于任何i，x'[i]和y'[i]无法 两个都是空隙。
[0054] 空位是X'或y'中的任一者中的连续空隙的最大子串。比对A可W包含W下Ξ种区 域：（i)匹配的对(例如，义'[。=7'[。）；。。不匹配的对，（例如，义'[。辛7'[。，并且两者 均不是空隙）；或（III)空位(例如，或者或者是空位）。在某些实施例中， 仅匹配的对具有高的正评分a。在一些实施例中，不匹配的对总体上具有负评分b，并且长度 r的空位还具有负评分g+rs，其中g，s<0。对于DNA，一个通用评分方案(例如，由BLAST所使用 的)使得评分3=1、评分b = -3、g = -5和s = -2。比对A的评分是所有的匹配的对、不匹配的对 和空位的评分的总和。X和y的比对评分可W被定义为在X和y的所有可能的比对之中的最大 评分。
[0055] 在一些实施例中，任何对具有由置换概率的4X4矩阵B定义的评分a。例如，B(i，i) =1和0<B(i，j)i<〉j<l是一个可能的评分系统。例如，在与颠换相比变换被认为是更加生物 学上可能的情况下，矩阵B可包含B(C，T) = 7和B(A，T) = 3,或者期望的或由本领域中已知的 方法确定的任何其它组值。
[0056] 根据本发明的一些实施例，比对包含成对比对。一般来说，成对比对包括一一具有 m个字符的序列Q(查询)和η个字符的参考基因组Τ(目标)一一找出和评估Q和T之间的可能 的局部比对。对于任何l<i<n和l<j<m，计算T比..i]和Q[k. . j]的最大可能比对评分（即，在 位置i处结束的T的任何子串和在位置j处结束的Q的任何子串的最优比对评分），其中h<i且 k<j。运可W包含检查所有具有cm个字符的子串，其中C是根据相似模型的常量，并且将每个 子串单独与化k对。每个比对被评分，并且具有优选的评分的比对被接受为比对。本领域中 的技术人员将理解，存在序列比对的精确算法和近似算法。精确算法将找出最高评分的比 对，但是在计算上会昂贵的。两种众所周知的精确算法是尼德曼-翁施算法(Needleman- Wunsch)(《分子生物学杂志（J Mol Biol)》，48(3) :443-453,1970)和史密斯-沃特曼算法 (Smith-Waterman)(《分子生物学杂志（J Mol Biol )》，147( 1): 195-197，1981;《数学进展 (八(1￥.111]\&1地.）》20(3)，367-387，1976)。后藤(6〇1〇}1)(《分子生物学杂志〇1〇161〇1)》， 162(3) ,705-708,1982)对史密斯-沃特曼算法的进一步改进将计算时间从0(m2n)减少到0 (mn)，其中m和η是被比较的序列大小，该改进更可修正为并行处理。在生物信息学领域，正 是后藤的改良算法通常被称为史密斯-沃特曼算法。史密斯-沃特曼方法用于将较大序列集 与较大参考序列进行比对，因为可更普遍并且更便宜地获得并行计算资源。参见例如亚马 逊(Amazon)的云计算资源。本文所参考的所有期刊文章 W其全部内容W引用的方式并入。
[0057] 史密斯-沃特曼(SW)算法通过奖励序列中的碱基之间的重叠并且对序列之间的空 位罚分来比对线性序列。史密斯-沃特曼算法还与尼德曼-翁施算法不同，不同之处在于SW 不要求较短序列跨越描述较长序列的字母字符串。也就是说，SW不假定一个序列是另一个 序列的全部内容的读数。此外，因为SW并不一定找出横跨字符串全长的比对，所W局部比对 可W在两个序列内的任何地方开始和结束。
[0058] 在一些实施例中，根据点矩阵法、动态规划法或单字法，进行成对比对。动态规划 法一般实施史密斯-沃特曼(SW)算法或尼德曼-翁施(NW)算法。根据NW算法的比对总体上根 据具有线性空位罚分d的相似矩阵S(a，b)(例如，诸如前述矩阵B)对比对的字符进行评分。 矩阵S(a，b)总体上提供置换概率。SW算法类似于NW算法，但是任何负评分矩阵单元被设置 为0。在美国专利5,701，256和美国公布2009/0119313中更详细地描述了 SW算法和NW算法及 其实施方式，该专利和公布W其全部内容W引用方式并入本文。
[0059] 实施史密斯-沃特曼算法的一个版本的比对程序是MUMmer，MUMme;r可W从由 Geeknet (弗吉尼亚州费尔法克斯(化irfax，VA))维护的Source化巧e网站商购获得。MUMmer 是用于快速比对基因组规模序列的系统（KurtZ、S等人，《基因组生物学（Genome Biology)》，5 :R12(2004) ;Delcher、A.L.等人，《核酸研究（Nucl. Acids Res.)》，27:11 (1999))。例如，MUMmer 3.0可W在2.4G化Linux台式计算机上使用78MB存储器在13.7秒内 找出在一对5-兆碱基基因组之间的所有20-碱基对或更长的精确匹配。MUMmer可W处理来 自鸟枪法测序计划的100个或1000个重叠群，并且将使用系统包含的NUCmer程序将100个或 1000个重叠群与另一组重叠群或参考进行比对。如果对于DNA序列比对来说物种太相异而 不能检测相似性，贝化ROmer程序可W基于两者输入序列的六框翻译生成比对。
[0060] 其它示例性比对程序包含:核巧酸数据库的高效大规模比对化fficient Large- scale Ali 即 ment of Nucleotide Databases 化 LAND)) 或序列和变异的同感评估 (ELANDv2component of the Consensus Assessment of Sequence and Variation (CASAVA))软件(加利福尼亚州圣迭戈的Illumina(Illumina，San Diego，CA))的ELANDv2部 件；Real Time Genomics公司的实施基因组学（RTG)调查器（加利福尼亚州旧金山（San FYancisco , CA)) ;Novocraft的Novoalign(马来西亚雪兰義州（Selangor .Malaysia)); Exonerate,欧洲生物信息研究所化uropean Bioinformatics Institute)(英国辛克斯顿 (Hinxton ,υΚ)) (Slater，G.和Birney ,Ε.，《BMC生物信息学（BMC Bioinfo;rmatics)》6:31 (2005))，都柏林大学的Clus1:al Omega,(爱尔兰都柏林(Dublin, Ireland)) (Sievers F等 人，《分子系统生物学(Mol Syst Biol)》7，adicle 539(2011));都柏林大学的Clus化IW或 011131日1乂（爱尔兰都柏林（〇11131;[]1,^61日]1(1))化日'10]11.4等人，《生物信息学 (Bioinformatics)》，23,2947-2948(2007);和FASTA，欧洲生物信息研究所化uropean Bioinformatics Insti1:ute)(英国辛克斯顿巧inxton,UK)) (Pearson W.R等人，《美国国家 科学院院刊（PNAS)》85(8) :2444-8( 1988) ;Lipman，D.J.，《科学（Science)》227(4693): 1435-41(1985))。
[0061] 如上面所论述的，当将RNA测序读数与有向非循环经注释的参考基因组进行比对 时，实施SW比对算法或（下面进一步更详细地论述的）的改良的版本可W是优选的或期望 的。
[0062] 根据W下方程式（1)，对于表示长度η和m的两个字符串的nXm矩阵H，易于表达SW 算法：
[006；3] Hk〇 = H〇i = 0()(ffO<k<nS〇<l<m) (1)
[0064]出j = max {出-1, j-1+s (ai，bj)，出-广Win，出,j-广Wdei，0}
[00化](对于1 < i < η?? < j < m)
[0066] 在W上方程式中，s(ai，bj)表示匹配奖分（当ai = b川寸）或不匹配罚分（当ai辛bj 时），并且对插入和缺失分别给出罚分Win和Wdel。在大多数情况下，所得矩阵具有为零的许多 元素。运种表示使得更容易在矩阵中从高到低、从右到左回溯，因此识别比对。
[0067] 在已经用评分完全填充矩阵后，SW算法实行回溯W确定比对。开始于矩阵中的最 大值，算法将基于Ξ个值化1-U-1、出或出J-1)被用于计算每个单元的最终最大值来进行 回溯。当达到零时回溯停止。最佳评分比对可W包含大于插入和缺失的最小可能数量的可 能数量，同时包含远小于置换的最大可能数量的可能数量。
[0068] 当作为SW或SW-后藤应用时，该技术使用动态规划算法来实行分别具有大小m和η 的两个字符串S和A的局部序列比对。此动态规划技术采用表或矩阵来保存匹配评分并避免 对连续单元的重新计算。可W根据序列的字母索引字符串的每个元素，也就是说，如果S是 字符串 ATCGAA，则 S [ 1 ] = A。
[0069] 代替将最佳比对表示为出,^上面），在W下方程式（2)中将最佳比对表示为B[j， k]:
[0070] B[ j，k]=max(p[ j，k]，i[ j，k]，d[ j，k]，0)(对于0<j <m，0<k< η) (2)
[0071] 在W下方程式（3)-(5)中概述最大值函数B[ j，k]的变量参数，其中MI5141'畑_ 阳 NALTY、MATCH_B0NUS、IN沈RTI0N_阳NALTY、DELETI0N_阳NALTY和0PENING_阳NALTY都是常 量，并且除MATCH_B0NUSW外都是负数。通过W下方程式(3)给出匹配变量参数p[j，k]:
[0072] 如果S[ j]辛A[k]，则p[ j，k] = max(p[ [ j-1，k-l]，i[ j-1，k-l]，d[ j-1，k-l]) + MISMATCH_PENALTY(3)
[0073] 如果S[ j]=A[k]，则p[ j，k]=max(p[[j-l，k-]，i[ j-l，k-l]，d[j-l，k-l])+MATCH_ BONUS,
[0074] 通过W下方程式(4)给出插入变量参数i[j，k]:
[0075] i[j，k]=max(p[j-l，k]+(PENING_PENALTY，i[j-l，k]，d[j-l，k] + (4)
[0076] 0 阳 NINGJENALTY)+INS 邸 TIONJENALTY
[0077] 并且缺失变量参数d[j，k]通过W下方程式(5)给出：
[007引 d[j，k]=max(p[j，k-l]+0阳NING_PENALTY，i[j，k-l] + (5)
[0079] 0阳NING_PENALTY，d[j，k-1])+DELETION_PENALTY
[0080] 对于所有Ξ个变量参数，将[0,0]元素设置为零W确保回溯完成，即，p[0,0] = i [0,0] =d[0,0] =0〇
[0081] 评分参数是稍微任意的，并可经调整W实现计算的行为。DM的评分参数设置的一 个示例（Huang,《第3章：生物序列比较和比对（Bio-Sequence Comparison and Alignment)》，当前顶级比较分子生物学（Curr Top Comp Mol Biol.)丛书，马萨诸塞州剑 桥市(Cambridge,Mass.):麻省理工学院出版社(The Μ口 Press) ,2002)将是：
[0082] MATCH-BONUS:10
[0083] MISMATCH_PENALTY：-20
[0084] IN沈RTI0N_PENALTY:-40
[0085] 0阳NING_PENALTY:-10
[0086] DELETI0N_PENALTY：-5
[0087] W上空位罚分（INS邸TI0N_阳NALTY、0PENING_阳NALTY)之间的关系有助于限制空 位开放的数目，即，支持通过设置高于空位开放成本的空位插入罚分来将空位集合在一起。 当然，MISMATCH_PENALTY、MATCH_BONUS、INSERTION_PENALTY、OPENING_PENALTY 和 DELETIONJENALTY之间可能存在可替代的关系。
[0088] 在一些实施例中，本发明的方法及系统并入多维比对算法。本发明的多维算法提 供了序列信息的"向后看"类型分析(如在史密斯-沃特曼中），其中通过包含多个通路和多 个节点的多维空间进行向后看。多维算法可W被用于将RNA测序读数与DAG类型参考进行比 对。该比对算法通过关于包含在DAG(例如，参考序列构造)上的位置处的每个序列识别最大 评分为Cl,J识别最大值。事实上，通过在先前位置处"向后"看，有可能横越多个可能的路径 识别最佳比对。
[0089] 在上面所论述的读数(亦称"字符串"）和有向非循环图表化AG)上执行本发明的算 法。出于定义该算法的目的，假设S是要比对的字符串，并且假设D是S将进行比对的有向非 循环图表。W从1开始的索引对字符串S的元素加括号。因此，如果S是字符串ATCGAA，那么S [1] =A、S[4] =G 等。
[0090] 在某些实施例中，对于DAG，节点的序列的每个字母将被表示为独立元素 cLd的前 趋被定义为：
[0091] (i)如果d不是其节点的序列的首字母，则其节点中在d之前的字母是其（唯一）前 趋；
[0092] (ii)如果d是其节点的序列的首字母，则为d的节点的父节点的任何节点（例如，基 因组中上游的所有外显子)的序列的最后一个字母是d的前趋。
[0093] 所有前趋的组继而表示为P[d]。
[0094] 为了找出"最优"比对，算法寻求M[j，d](S的前j个元素与在d之前(且包含d)的DAG 的一部分的最佳比对的评分）的值。该步骤类似于在W上的方程式1中发现出,具体来说， 确定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* ] +DELETEJENALTY}
[009引 i=M[ j-1，d]+IN 沈 RT_PENALTY
[0099] 对于在 P[d]中的 p*，如果 S[j]=d，则 a=max{M[j-l，p*]+MATCH_SCORE};
[0100] 对于在P [ d ]中的 p*，如果 S [ j ]辛 d，则max {M [ j -1，p* ] +MI SMATCH_PENALTY}
[0101] 如上面所描述，e是S的前j个字符与DAG直到（但不包含)d的部分的比对加上另外 的DELETE_PENALTY的最高值。因此，如果d不是节点的序列的首字母，则仅存在一个前趋p， 并且S的前j个字符与DAG(直到且包含P)的比对评分等效于1。'，9]+061^了6_阳魁1；1￥。在其 中d是其节点的序列的首字母的实例中，可W存在多个可能的前趋，并且因为DELETE, PENALTY是常量，所W求^。'柳]+061^了6_阳魁1；1'門的最大值和选择具有与8的前^'个字符 比对的最高比对评分的前趋相同。
[0102] 在方程式(6)中，i是字符串S的前j-1个字符与DAG的直到并且包含d的部分的比对 加上INSERTJENALTY，其类似于SW中的插入变量参数的定义(参看方程式1)。
[0103] 另外，a是S的前j个字符与DAG的直到但不包含d的部分比对的最高值，加上或者 MTCH_SC0RE(如果S的第j个字符与字符村目同）或者MISMATCHJENALTY(如果S的第j个字符 与字符d不同）。如同e-样，运意味着如果d不是其节点的序列的首字母，则仅存在一个前 趋，即P。运意味着a是S的前j-1个字符与DAG(直到并且包含P)的比对评分，即，M[ j-1，p](再 加上或者MI SMATCH_PENALTY或者MATCH_SC0RE)取决于d与S的第j个字符是否匹配。在其中d 是其节点的序列的首字母的实例中，可W存在多个可能的前趋。在该情况下，求{M[j，p*] + MI SMATCHJENALTY或MATCH_SC0RE}的最大值和W下相同：选择具有与S的前j-1个字符的最 高比对评分(即，候选M[ j-1，p*]变量参数的最高值)的前趋并且取决于d与S的第j个字符是 否匹配而加上或者MISMATCHJENALTY或者MATCH_SC0RE。
[0104] 再次，与在SW算法中相同，罚分例如DELETE_PENALTY、INSERT_PENALTY、MATCH_ SCORE和MI SMATCH_PENALTY可W被调整W促进与较少空位等的比对。
[0105] 如W上方程式中所描述，该算法通过不仅计算该元素的插入、缺失和匹配评分，而 且向后看(逆着DAG的方向巧化AG上的任何之前节点W找出最大评分，来找出每个读数的最 佳(即，最大)值。因此，该算法能够详细研究穿过DAG的含有已知突变的不同路径。因为图表 是有向的，所W逆着图表方向移动的回溯跟随朝向图表起点的优选异构体，并且最佳比对 评分识别高度确定性的最可能比对。
[0106] 图1示出与DAG比较的读数的图形表示。图1的顶部区域呈现序列读数"ATCGAA"W 及参考序列TTGGATATGGG(序列ID号：1)和已知的插入事件TTGGATCGAATTATGGG(序列ID号： 2)，其中插入是带下划线的。图1的中间示出读数、序列ID号：1和序列ID号:2之间的关系，保 持每个序列线性W辅助观察关系。图1的底部区域示出使用DAG构造的比对。在所描绘的DAG 中，即使沿不同路径，但可W通过从DAG的5'端到DAG的3'端读取来读取序列ID号1和序列ID 号2两者。如所描绘的，序列读数被示出为与上部路径比对。
[0107] 图2示出对应于比较的实际矩阵。相似于史密斯-沃特曼技术，本发明的所例示算 法识别最高评分，并且实行回溯W识别读数的正确位置。图1和图2还突出本发明产生字符 串与该构造的实际匹配。在其中序列读数包含未被包含在DAG中的变异体的情况下，将通过 空位、插入等报告经比对的序列。
[0108] 将RNA测序读数与经注释的参考的DAG进行比对允许发现转录组中的异构体。另 夕h已经发生的不同形式的可变剪接的可能性(例如，一个密码子已在相关外显子中的一个 中被跳过的可能性何W被"免费"预期。
[0109] 外显子可自动地分成两个节点，一个对应于第一密码子，并且另一个对应于外显 子的其余部分。运允许在DAG中表示两种类型的剪接，即使还未观测到它们。鉴于一个密码 子的与比对的差值的价值将往往不足W惩罚候选比对W诱导我们寻找替代方式，预测运些 情况可W是有用的。（最大的此罚分将是Ξ个插入或Ξ个缺失的罚分，并且由于机率，所W 常常罚分将是较小的。）
[0110] 将转录组表示为MG准许将异构体自然地解释为穿过DAG的路径。与往常一样，当 系统和方法表示数学结构中的物理项时，益处是大量的。例如，本发明的系统和方法允许序 列W直观的方式可视化。例如，可W通过示出穿过DAG的两个路径来比较序列。
[0111] 另外，本发明的方法有益于利用图表定理和算法。穿过DAG的路径的分析是在现代 的数学和计算机科学中充满生机的研究区域。因此，利用该研究来改进转录组研究的能力 是有价值的。例如，如果边缘被给定与由读数跨越的接合点的概率对应的权重，则路径的总 权重可W携载关于实现不同异构体的先验可能性的信息。穿过DAG的路径最大化(或最小 化)权重是已知的，并且最短路径算法可W被用于快速地找出那些最短路径。
[0112] 图3提供了本发明的方法的图式。如上面所论述的，一旦DAG内置有W规范的顺序 连接的外显子，则获得双端读数。然后，双端读数与DAG比对，并且在读数和包括表示RNA序 列的节点和连接多对节点的边缘的DAG之间找出比对(具有满足预定指标的比对评分）。预 定指标可W是例如最高评分比对。比对可W通过上述改良的史密斯-沃特曼算法进行。对于 一对双端读数，比对可W包含基于穿过包含比对的有向非循环数据结构的路径，识别转录 组中的转录本。通过使用插入大小作为比对的约束，可W有效地实行比对。双端读数可W被 用于构建DAG (例如，添加边缘）。
[0113] 该方法适合于分析由RNA测序获得的读数。在某些实施例中，有向非循环数据结构 基本上表示至少一个染色体的所有已知的特征。可W在找出比对之前创建和存储有向非循 环数据结构。在一些实施例中，通过将每个序列读数与穿过有向非循环数据结构的可能路 径中的至少大部分进行比较来找出比对。方法可W包含基于找出的比对将多个序列读数装 配成重叠群。
[0114] 3.使用插入大小作为比对的约束
[0115] 其中插入大小信息可W与DAG组合的一种方式是使用插入大小作为比对的约束。 通过使用该方法，速度和精确性都有可能增加。因为该方法防止丢弃与比对相关的重要信 息，所W与单独比对双端配偶相比，该方法可W更精确。应当注意的是，独立比对两个双端 配偶近似地相当于在基因组中L个可能的位置(如果基因组具有长度L)之间选择两个位置。 因此，比对包括在大约L~2种可能性之间进行选择。然而，第一比对被用于强有力地将另一 比对约束到在第一比对周围的常量大小k的一些区域，然后仅有大约吐种可能性一一在L而 不是L~2的量级上。运将更加快速。
[0116] 图4-图6例示使用插入大小信息W约束比对。
[0117] 图4例示任意结构的假设的DAG和定义插入大小k的一对假设的双端读数。在如所 绘制的DAG中，节点是方块，并且边缘是线。在所描绘的实施例中，从具有均值k bp的近似正 态分布获取插入大小k。将该对的第一成员与DAG比对(例如，使用上述史密斯-沃特曼或改 良的史密斯-沃特曼）。一般来说，双端包括片段，该片段的每端上具有衔接子，并且该片段 具有从片段-衔接子边界朝向彼此朝内延伸的双端读数。在双端读数之间可存在空位。插入 长度可W意味着片段的长度一一即，从一个片段-衔接子边界到另一个片段-衔接子边界的 距离。参见例如Lindgreen，2012,《AdapterRemoval :下一代序列读数的简单清理 (AdapterRemoval:easy cleaning of next-邑eneration sequence reads)》，《BMC石开究参己 要《(BMC Res Notes)》5:337(例如，Lindgreen，图1，画面C)，其内容W引用方式并入。
[0118] 在图5中，最左边黑色方块表示外显子，该对的第一成员在外显子处比对。确定了 一组距离，它们占期望的插入大小值中的大部分(例如，99.95%)。例如，其中插入大小被正 态分布，该组距离可W是k的3个标准差内的那些距离。所确定的该组距离被用于选择下游 节点，且因此选择候选路径。所选择的下游节点表示该对双端读数的第二成员的潜在比对 位点。在图5中，该组下游节点还被填充在黑色方块中。开始于上游节点且结束于下游节点 中的一个的路径是候选路径，将该对双端读数与候选路径比对。
[0119] 图6例示来自图4的DAG的候选路径。应当注意的是，图6还呈现包含在图4的DAG内 的DAG。图6中例示了所有候选路径且仅例示了候选路径。可W将定义插入大小k(来自图4) 的假设的双端读数对与图6中呈现的该DAG进行比对W节省计算量。在没有认识到对于比对 的该约束的情况下，该对的每个读数可W沿着长度L的基因组并且在复杂的DAG(例如，图4 中的DAG近似地为横越其长度4重的）中潜在地与任何位置比对，可存在多得多的位置（即， 平均来说，如果DAG是横越任何基因位点η重的，则比对一对双端读数相当于从(nL)~2个选 项中进行选择）。根据所例示的方法实行与DAG的比对可W使非常大的计划计算上成本小得 多，产生具有来自RNA测序计划的高容量双端读数数据的大得多的吞吐量。
[0120] 4.对可变剪接的变种的敏感度。
[0121] 当读数跨越外显子的接合点时，当前技术使得难W或不可能在几种类型的情况之 间进行区分。一种此情况是已发生突变的情况。即使在没有突变存在的情况下，当因为已经 发生一种形式的可变剪接(其中密码子的一部分已被包含或跳过），所W读数可W显得与接 合点的任一侧都不匹配时，出现另一种此情况。此外，测序误差呈现另一种此情况。
[0122] 考虑W下序列作为示例：
[01U] AUGCAUUUCCGUGGAUAGCGAUGCAUACUCACGAUG AUGAAAAUGCAUCAGAAAUAG(序列ID号： 3)
[0124]其中明文区域表示外显子，第一带下划线的部分表示内含子，并且第二带下划线 的部分是由于可变剪接所W可能或可能不包含在第二外显子中的区域。因此，成熟RNA可采 取W下两种形式中的任一种：
[01 巧]AUGCAUUUCCGUGGAUAGAUGAAAAUGCAUCAGAAAUAG (序列 ID 号:4)
[0126] AUGCAUUUCCGUGGAUAGAUGCAUCAGAAAUAG (序列 ID号:5)
[0127] 现在假设我们获得W下读数:TGGATAGATGAA(序列ID号：6)。此读数可W具有在重 要的方面不同的若干因果关系历史记录。
[0128] 图7呈现包括运些序列的Ξ个可替代的历史记录。在横越图7的顶部第Ξ个所描绘 的第一情况中，采取来自序列ID号:4的读数，并且没有突变或测序误差已发生。在第二情况 中，采取来自序列ID号：5的读数，并且测序误差已发生。在第Ξ情况中，采取来自转录序列 ID号:5的读数，并且读数精确地反映突变。
[0129] 如果线性参考序列仅由明文区域（图7中短划的下划线）组成，则很有可能比对评 分将非常高一一通过很多匹配升高且仅遭受一个不匹配罚分一一使用中的算法将不使比 对标记为可能错误。并非检查用于剪接信息的一些其它文件或数据库，而是将简单地假定 突变或单个测序误差。与DAG参考进行比对的确较好。
[0130] 图8描绘了序列ID号:3与DAG参考（由于没有了内含子，所W现在出现在图8中作为 序列ID号:4)的比对。运里，因为精确的比对将具有与该DAG任何可能的比对的最高评分，所 W精确的比对立即被恢复。由于DAG的本质，图8还描绘了：
[0131] AUGCAUUUCCGUGGAUAGAUGCAUCAGAAAUAG (序列 ID号:5)
[0132] 鉴于外显子的小部分遭受由可变剪接引起的包含或除去的频繁程度，运是很一般 的现象并且将常常发生。当运些情况出现时，包括线性参考的系统将通常无法检测错误的 特性描述，而防止运些错误的唯一方式是与如此多的线性参考进行比对，W至使过程过分 地昂贵(在时间和金钱两方面）。
[0133] 5.正确地使用双端数据:改良DAG。
[0134] 本发明提供另外地利用插入大小和准确比对之间的关系。当处理插入大小作为对 给定比对或一对比对的约束时，可W检视当前范例。运是对该情况的不完整的描述。将插入 大小检视作为对于比对和转录本的组合的约束可W是更精确的。
[0135] 本发明的方法包括不仅使用插入大小作为对比对的品质控制，还使用插入大小作 为对于转录本的品质控制。如果两个双端配偶中的每个在转录DAG中具有高评分比对，并且 如果该对比对对应于非常不合情理的插入大小，则与关于比对中的每个的正确性调整我们 的置信相比，关于转录调整置信(例如，通过将边缘添加到图表)常常将是更合理的。
[0136] 运可W通过一-无论何时两个双端配偶与DAC比对一一计算插入大小或由该比对 暗示的大小来实施。（因为读数可能与若干特征一致，所W可能包括不止一个大小。）如果最 小的此插入大小和插入大小的期望值的商大于一些常量(可能为2)，则检查DAGW找出一对 未连接的节点，如果该对未连接的节点连接，则将必然伴有异构体的存在，对于该异构体， 该异构体的所暗示的插入大小和插入大小的预期值的商在0.5和1.5之间。如果此对节点存 在，则用边缘连接它们。如果若干个此对存在，则选择产生最接近于1的比率的那个对并且 连接它。
[0137] 在某些方面中，本发明提供一种分析转录组的方法，该方法包括从转录组获得至 少一对双端读数、找出在该对的第一读数和在有向非循环数据结构（该数据结构具有表示 RNA序列诸如外显子或转录本的节点和连接每对节点的边缘）中的节点之间的具有最佳评 分的比对，识别包含节点的候选路径，该节点通过具有与该对双端读数的插入长度大体上 类似的长度的路径连接到下游节点，和将双端读数与候选路径比对W确定最佳评分比对， W及改良DAGW在必要时包含新的边缘或者优选地创建由一对双端读数指定的候选路径。 本领域中的技术人员将认识到，可W由任何合适的指标定义具有与该对双端读数的插入长 度大体上类似的长度的路径。例如，大体上类似可w意指相同数量的核巧酸，或者可w优选 地意指+10或+100核巧酸。在某些实施例中，大体上类似可W由用户定义W意指在长度上彼 此差距在10%内，或优选地彼此差距在25%内，或彼此差距在50%内，或彼此差距在100% 内。在优选的实施例中，长度大体上类似意指一个的长度不超过另一个的长度的两倍。该方 法适合于许多双端读数，诸如来自cDNA片段的那些读数（例如，表示单个转录组的cDNA片 段）。节点和下游节点可W表示一对外显子，该对双端读数从该对外显子获得。方法可W包 含基于所确定的最佳评分比对来识别异构体。在一些实施例中，所识别的异构体是新颖的 异构体。可W更新有向非循环数据结构W表示新颖的异构体(例如，通过添加新的边缘）。方 法可W包含从包括该对双端读数的任何比对计算中排除不是候选路径的任何路径。本发明 的其它方面通过例如可用于推断异构体频率的系统和方法提供了表达分析。
[0138] 6.使用双端数据:推断异构体频率。
[0139] 强有力的技术不是起因于使用双端结构W正确地比对，而是起因于使用比对W用 从比对推断插入大小的中间步骤推断异构体频率。如果多对外显子出现在样本中的频率是 已知的（参见例如化rrow，2012,《基因组研究(Genome Research)》22:1760)，则那些频率可 W被用于构造异构体出现在样本中的频率。此外，不同的多对外显子通常将被不同的距离 分隔开。因此，在双端配偶之间观测到的距离可W(稍微间接地)被用于跟踪外显子频率。
[0140] 例示性示例如下。对于每对双端配偶，该对可W首先可选地与DAG参考进行比对。 注意且跟踪每个对成员上的参考点（例如，最左边位置)之间的距离。运可W从与DAG的比对 获得，虽然应当注意推断异构体频率的方法可W不需要将该对与DAG比对。
[0141] 用于注意和跟踪在每个被比对的对的参考点之间的距离的一个合适的方法是通 过使用散列表。在距离散列中，关键字可W都是在1和1000之间的距离，并且每个关键字的 值被初始化为0。然后，可W通过查找对应于距离的关键字的值并递增该值来完成递增。应 注意，可W极其快速地完成散列表的运类运算一一实际上，在常数时间内。另外，应该注意 的是，通过现有的生物信息学软件包诸如BioPerl很好地支持散列表。然后，每对外显子与 对应于最接近外显子之间真实距离的距离的值的总和相关联。（也就是说:将为接近那些外 显子之间距离的距离所观测到的所有"命中次数"加起来。在样本中，相对于的所有对的外 显子，对的该数量和对的总数量的商将非常接近该对外显子发生的频率）。
[0142] 可从该数据用代数方法确定异构体频率。方程式的优选的形式可W取决于相关异 构体的结构，但一般技术获得如下：
[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] 当左侧值是从比对已知的时，可W解算异构体频率。一般来说，本发明的方法包含 编写将外显子-对-频率与异构体-频率相关的方程式，并且使用前者w获得后者。本领域中 的技术人员将认识到用于解算所描绘的方程式的代数方法。为了给出例示性示例，应指出 上面的代数运算产生：
[0149] f=[AB]-g，
[0150] W及：
[0151] h=[BD]-g
[0152] 因此；
[0153] [BC]=f+[BD]-g
[0154] W及；
[0155] [BC] = [AB]-g+[BD]-g
[0156] 重排提供：
[0157] [BC]/2([AB] + [BD])=g
[0158] 由于[BC]、[AB]和[BD]都是从比对已知的，因此计算g。然后解算f和h是数学上最 简单的。
[0159] 在一些情况下，可W出现"自由度"问题，其中要解算的变量的数量不小于输入变 量的数量。用于解决该显而易见的问题的任何合适的技术可W被采用，并且可W取决于特 定的情况。可W找出的合适技术包含使变量标准化;解算超过一组异构体频率且事后选择 一个;在外部获得至少一个输入值(例如，从生物样本获得），或者提供输入信息(例如，从文 献或假设提供）。
[0160] 重要的是要注意关键步骤的图形解释:假设产生在X轴上具有距离和在y轴上具有 频率的图表。然后，外显子-对将出现作为该图表上的峰值，并且经常那些峰值中的一些还 将容易地可解释为异构体-频率。
[0161] 因此，本发明提供通过从转录组获得多个双端读数对且可选地将那些读数与参考 数据结构进行比对，推断转录组中异构体频率的方法。数据结构包含来自基因组的每个被 表示为节点的多个外显子。多对节点或外显子在它们的近端处通过边缘连接。确定双端读 数的所比对的对之间的距离和所比对的对之间的距离的频率。确定一组异构体路径和异构 体频率，W使得表示W异构体频率穿过结构的异构体路径引起多对外显子W所比对的对之 间的距离的频率被包含在该异构体路径中。所确定的该组异构体路径中的每个路径可W表 示存在于生物体中的转录异构体。该比对的对之间的距离可W是从每对的第一比对成员的 上游端到该对的第二比对成员的下游端的路径长度。在某些实施例中，确定该组异构体路 径包含发现新颖的异构体路径。可W更新参考数据结构W包含新颖的异构体路径。
[0162] 可W用代数方法例如通过解算一组线性方程式实行确定该组异构体路径。在一些 实施例中，在解算该组线性方程式(例如，所W可W确定线性方程式中的其余的变量)之前， 从生物学数据获得异构体频率中的一个。另外地或可替代地，在确定一组异构体路径和异 构体频率之前，可W使在所比对的对之间的距离的频率标准化。
[0163] 7.合并先验外显子计数分布和异构体计数分布。
[0164] 在一些实施例中，本发明利用其中有用于异构体计数的相对频率的先验数据的情 况。转录本的异构体分布可W被推断与蛋白质异构体的分布相关。即使不完全地，也已经彻 底地研究了关于给定基因的外显子计数分布和异构体计数分布的问题，W及运些分布对关 于序列的其它事实(例如，其是编码序列还是非编码序列)敏感的方式。
[0165] 可W从文献或者通过实行分析获得异构体频率。例如，可W从如在上面第6节"使 用双端数据:推断异构体频率"中所描述的双端读数RNA测序读数推断异构体频率。在一些 实施例中，来自双端读数数据的插入大小用于推断异构体频率，如下面第8节"使用插入大 小更有效地推断异构体概率"中所描述。在某些实施例中，从源诸如数据库或文献获得异构 体频率。
[0166] 图9表示由SwissProt氨基酸序列的研究获得的蛋白质异构体分布，并且是来自 Nakao等人，2005，《人类可替代的蛋白质异构体的大规模分析:模式分类和与亚细胞定位信 号的相关性（L曰r邑e-sc曰le 曰η曰lysis of hum曰η 曰Iterative protein isoforms:p曰ttern classification and correlation with subcellular localization signals)》，《核酉畫 研究(Nucl Ac Res)33(8):2355-2363的第Ξ个图的再现。另外，化ang等人已经使用编程方 法推导各种可变剪接形式的定量分布。化ang等人，2005,《来自EST数据库的可变剪接的定 量分析中可变剪接图表的应用（The a卵lication of alternative splicing胖日地S in qu曰ntit曰tive analysis of 曰Iterative splicing form from EST d曰t曰b曰se)》，《计算 机技术应用国际期刊（Int J.Comp.App 1.Tech)》22(1): 14。
[0167] 由图9呈现的信息表示可W被用于辅助将读数映射到DAG的异构体计数分布。类似 地，可W获得外显子计数分布信息。
[0168] 例如，有每转录本的外显子的数量的显示的分布。参见例如化rrow等人，2012， 《GENC0DE:用于编码计划的参考人类基因组注释（GENC0DE: The reference human genome anno1:ation for The ENCODE Project)》，《基因组研究（Genome Res)》22:1760-1774。蛋白 质-编码基因位点处的转录本在每转录本4个外显子处显示峰值，而IncRNA在两个外显子处 示出不同的峰值(参见HarrOW等人，2012的图2)。可W用给出相对频率而不是绝对频率的y 轴概率性地解释来自化rrow等人的外显子计数分布或来自化kao等人的异构体分布此信息 可W被用于告知关于候选比对对应于实际核巧酸序列的似然性的判定。现有技术主要依赖 于比对的评分。本发明提供包含异构体数量的先验知识的方式。
[0169] 例如，对于大的读数组，75%可W明确地映射到一组8个异构体{II、12.....18}中 的区域，而其另外25%可W被解释为或者映射15或者映射到对应于另一异构体19的外显子 的一些其它组合。在该情况下，各种类型的事实都与如何比对读数的最后25%的决策相关。 此类相关事实包含（i)与15比对和与19比对的评分；（ii)我们可能具有的关于调控机制的 任何进一步的知识；W及(III)有8个或9个异构体是否更可能为先验的。
[0170] 当将读数映射到DAG参考时，本发明合并异构体分布信息。可若干种不同的方 式使用异构体频率分布信息。在映射读数中包含先验异构体数量的第一示例性方法包含： (1)为基因组的相关区域构造 DAG(该DAG的边缘可W是加权的或未加权的及(2)将DAG的 每个边缘与指示边缘的使用的值相关联。如果W跨越该边缘的方式所比对的读数的数量小 于一些变量U，则设置使用"异常"。该变量直观地反映边缘是否已经被"使用"。估计该变量 的一种方式可能是通过将其设置为等于(样本中读数的数量)*(外显子长度的合理估计和 平均读数长度的商)*(.〇〇5)。
[0171] 进一步地，该方法包含（3)将基因组的相关区域与对应于使用阔值T的变量相关 联，可能设置T等于用于图表的边缘的数量的期望值加上用于图表的边缘的数量的标准差。 运些数量可W通过W下步骤计算:估计DAG中边缘的数量，该DAG具有与为所讨论的基因而 期望的经验外显子频率、异构体频率或运两者对应的节点和最大路径的数量。经验频率可 W例如从化rrow等人，2012(用于外显子计数分布)或者从化kao等人，2005(用于异构体分 布)获得。
[0172] 该方法包含(4)改良用于史密斯-沃特曼或其它算法的评分方案;如果"使用的"边 缘（即，布尔型使用值是"正常的"或"真的"边缘）的总数N超过T，则通过一些罚分惩罚未使 用的边缘的遍历。可W使该罚分随着N增加而增加。用于罚分的合理的公式可W是max(0， (N-T)*(插入罚分)*2)。0项出现W避免负罚分一一即，W避免使用新边缘的加分一一但此 些加分在一些情形中可能是合适的，在该情形中，不需要采取0的最大值和后者表达的值。
[0173] 为了避免呈比对对读数被处理的顺序敏感的形式的不被接受的路径依赖性，可W 重复比对若干次且比较结果。运不影响算法的渐进运行时间。
[0174] 在映射读数中包含先验异构体数量的可替代的示例性方法包含：
[01巧](i)比对每个读数与DAG;
[0176] (ii)将所观测的异构体数量与期望的异构体数量(例如，给定关于基因的背景信 息，从异构体数量的先验概率分布获得的）比较;W及
[0177] (iii)如果所观测的数量和所预期的数量之间的差值超过某一阔值，则重新比对 被分配给最罕见的候选异构体的读数。在重新比对中，用根据差值的量值的权重将与罕见 的候选异构体的任何比对罚分。罚分可W是max(0，（N-T)*(插入罚分)*2)，其中N是"使用 的"边缘（即，布尔型使用值是"正常的"或"真的"边缘)的总数，并且T是使用阔值。在一些实 施例中，T是用于图表的边缘的数量的期望值加上用于图表的边缘的数量的标准差。
[0178] 通过应用使用异构体数量的先验概率分布的方法，将推动读数被映射到DAGW产 生接近先验分布的分布。
[0179] 8.更有效地使用插入大小W推断异构体概率。
[0180] 假设在不同的异构体中出现两个外显子E1和E2,并且在那两个异构体中外显子E1 和E2通过距离L1和L2被分隔开。将偶尔发生获得一对双端读数，其一个配偶映射到E1，并且 另一个映射到E2。然后，将该对归属于哪个异构体的选择可W被缩减为鉴于配偶的从外显 子边缘的偏移是读数之间的真实距离，是否推断L1或L2或那些距离中的一个移位适当的量 的问题。
[0181] 因为我们通常具有与双端测序运行中插入大小的分布有关的良好的背景信息，所 W将异构体-分配问题转换为关于插入大小的问题是有价值的。鉴于通常将没有其它理由 相信更有可能从一个异构体或另一个异构体采集读数，我们可W实行简单的贝叶斯计算， W根据分布中插入大小L1和L2的概率的比率将读数概率性地分配给两个异构体。
[0182] 运是对当前技术的显著改进，其该情况下总体上检视读数为携载用于两个异构体 的等量证据权重。因此，W所描述的方式合并插入大小是不忽略相关证明的方式。
[0183] 因此，在运些不明确的情况下，我们可W计数异构体，如下：
[0184] (1)计算插入大小的概率分布。（由P(s)指代具有插入大小S的读数的概率。）应注 意运仅需要计算一次，并且运可W在预处理步骤中完成。
[0185] (2)如果两个双端配偶中的每个映射在外显子内部：
[0186] (2a)制作包含那些外显子的所有异构体的列表。
[0187] (2b)将每个异构体与将由已经从该异构体导出的读数暗示的插入大小相关联。
[0188] (2c)为(2b)中所计算的插入大小中的每个计算P(s)。
[0189] (2d)将每个异构体与P(s)除W(2c)中所计算的所有P(s)值的总和的相关值相关 联。
[0190] 应该注意的是，本文所描述的技术将W不影响使用"下游"工具(也就是说，将读数 计数作为输入且给出其它有用的输出的工具)的能力的方式改进异构体的读数计数。
[0191] 9.构造"最大 DAG"。
[OW] 在一些实施例中，构造 DAG可W是有益的，其中，在此DAG中每对节点用从表示"较 早"或"左侧"外显子到"较迟"或"右侧"外显子的节点的边缘连接。
[0193] 使用Gencode或类似的经注释的转录组，创建节点W表示外显子和其它特征。 Gencode提供足够信息W将运些节点与此信息相关联作为节点的序列、其在参考序列中占 据的位置(并且因此其长度)等。因此，使用包括禪接到非暂时性存储器的处理器的计算机 系统，来自基因组的多个特征中的每个被表示为节点。对于该多个特征中的每对，创建边 缘，边缘通过该对的两个成员的近端连接该两个成员。此移动的生物学解释是猜测表示遵 守外显子出现在基因组中的顺序的那些外显子的每个组合的异构体的存在，或者至少促进 发现该异构体。对遵守外显子出现在基因组中的顺序敏感既有感应性的原因又有生物化学 的原因。
[0194] 例如，报告表明在多外显子前mRNA中的外显子总是按顺序保持。参见例如Black, 2005,《剪接难题的简单答案(A simple answer for a splicing conun化um)》，《美国国家 科学院院刊(PNAS)》102:4927-8。没有被任何特定的机制束缚，从理论上可W讲，剪接机械 设备辨识内含子的开始和结束，W及分支点腺嚷岭或接近剪接接合点的其它保守序列。剪 接机械设备或剪接体沿着前mRNA从5 '到3 '移动，移除内含子且剪接外显子。因此看来存在 用于W5'到3'顺序连接所有外显子的现象学W及文献支撑，在5'到3'顺序中，所有的外显 子看起来沿着基因组DNA链。还参见化ao等人，2006,《人类中外显子重复接收、外显子加扰 和反式剪接的生物信息学分析（Bioinformatic analysis of exon repetition,exon scrambling and trans-splicing in humans)》，《生物信息学(Bioinformatics)》22:692- 698。对于某些假设，本发明的方法遵循人类前mRNA被剪接成保留由基因组序列定义的外显 子顺序的线性分子的假设。
[01M]另外，在期望在有向非循环数据结构中包含外显子加扰的实例（参见例如 化；rrington等人，1985,《在伴刀豆球蛋白A的转译后修饰期间发生多肤接合(Polypeptide ligation occurs during post-translational modification of concanavalin A)》， 《自然(Nature)》313:64-67)的情况下，可W通过包含双重复外显子作为"虚构"大体上保持 外显子W它们的原始顺序剪接的假设(或者，在可替代的方案中，可W废除该假设）。
[0196] 因此，在例如基因组包含处于W下顺序:ABCDEFGHI的W下外显子的情况下
[0197] 使用短划线作为边缘且字母作为节点，完整的该组边缘可W被表示为：
[019 引 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]果运足W证明发现所有的异构体而不是一种，并且该一种异构体一一给出示例： ADE服
[0200] 贝IJ，借助仿真的或虚构的G-一称它为G'，该非循环数据结构可W是"钉齿状的"。 然后，可W由A-D-E-H-G'表示该异构体，并且结构不需要包含从昭化的任何边缘。然而，可 W发现可W忽略外显子加扰和由所有包含的边缘保持的基因组外显子顺序。
[0201] 可W通过使用来自经注释的转录组的输入数据，首先制备本发明的DAG。同样地， 从头开始，DAG将包含所有已知的异构体。如果制成最大的DAG，贝化AG将包含保持基因组外 显子顺序的所有可能的边缘，并且然后，将因此包含既表示已知的异构体又表示新颖的异 构体的节点和边缘。在存储器中节点和边缘被存储为有向非循环数据结构，并且方法包含 将双端序列读数与包括由边缘的子集连接的特征的子集的路径进行比对，从而识别对应于 路径的异构体。
[0202] 当将新的数据与DAG(例如，双端RNA测序读数）比对时，那些数据可W包含在经注 释的参考中未被识别的外显子。因为在该示例中外显子是新颖的，所W该外显子是异构体。 使用DAG作为参考理想地适合于该情况，因为将新颖的异构体与现有的DAG进行比对还可W 引导在DAG中创建新的节点的W表示新颖的外显子。因此，在所识别的异构体是新颖的情况 下，比对转录序列可W包括在有向非循环数据结构中创建新节点。还可W为剪接变异体创 建新节点，例如，其中一个外显子通过单个密码子不同于另一个外显子。
[0203] 优选地，在分析双端读数之前且独立于分析双端读数，初始地全部或部分地构建 最大DAG。因此，方法可W包含在比对转录本之前，为多个特征中的每对创建边缘。出于与有 向非循环数据结构的本质相关的原因，数据结构中表示的参考可W任意大的，例如，数百或 数千外显子、内含子或两者。
[0204] 10.本发明的系统
[0205] 图10例示用于实施本文所描述的方法的计算机系统401。本发明的系统可W包含 图10中所示的任一个部件或任何数量的部件。一般来说，系统401可W包含能够通过网络 415彼此通信的计算机433和服务器计算机409。另外，可W可选地从数据库405(例如，本地 或远程)获得数据。在一些实施例中，系统包含用于获得测序数据的仪器455,仪器455可W 被禪接到测序仪计算机451用于序列读数的初始处理。
[0206] 在一些实施例中，通过并行处理实行方法，并且服务器409包含具有并行架构的多 个处理器，即，能够将多个序列（例如，RNA测序读数)与参考序列构造(例如，DAG)进行比较 的处理器和存储器的分布式网络。系统包括多个处理器，运些处理器同时执行多个读数与 参考序列构造之间的多个比较。虽然可能有其它混合配置，但是并行计算机中的主存储器 通常或者在单个地址空间中的所有处理元件之间共享，或者呈分布式，即，每个处理元件具 有其自身的本地地址空间。（分布式存储器是指运样的事实:存储器W逻辑方式分布，但通 常暗示它也W物理方式分布。）分布式共享存储器和存储器虚拟化组合运两种方法，其中处 理元件具有其自身的本地存储器W及对非本地处理器上的存储器的存取权。对本地存储器 的存取通常比对非本地存储器的存取更快。
[0207] 其中可W用相等时延和带宽存取主存储器的每个元件的计算机架构被称为均匀 存储器存取(UMA)系统。通常，UMA系统只能通过共享的存储器系统来实现，其中存储器并非 W物理方式分布。不具有该特性的系统被称为不均匀存储器存取(NUMA)架构。分布式存储 器系统具有不均匀存储器存取。
[0208] 可W W若干种方式在硬件中实施处理器-处理器和处理器-存储器通信，运些方式 包含经由共享的（或者多端口的或者多路复用的）存储器、纵横开关、共享的总线或无数拓 扑的互连网络(包含星形、环形、树形、超立方体、胖超立方体(在一个节点处具有不止一个 处理器的超立方体））或η维网格。
[0209] 基于互连网络的并行计算机必须合并路由W实现在并非直接连接的节点之间的 消息传递。用于处理器之间的通信的媒体很可能在大型多处理器机器中分层。此类资源在 市面上可购买用于专用用途，或可W经由例如亚马逊云计算的"云"存取运些资源。
[0210] 计算机总体上包含通过总线禪接到存储器和输入-输出（I/O)机构的处理器，存储 器可W包含RAM或ROM,并且优选地包含存储指令的至少一个有形的非暂时性媒体，该指令 可执行W致使系统实行本文所描述的功能。本领域中的技术人员将认识到，根据需要或者 最适合于实行本发明的方法，本发明的系统包含一个或多个处理器（例如，中央处理单元 (CPU)、图形处理单元(GPU)等）、计算机可读存储设备(例如，主存储器、静态存储器等），或 其组合，它们通过总线彼此通信。
[0211] 处理器可W是本领域中已知的任何合适的处理器，诸如由英特尔（加利福尼亚州 圣克拉拉(San化Clara，CA))W商标XE0NE7出售的处理器，或由AMD(加利福尼亚州桑尼维 尔（Sunnyva 1 e，CA)) W商标0PTER0N 6200出售的处理器。
[0212] 根据本发明的输入/输出设备可W包含视频显示单元（例如，液晶显示器化CD)或 阴极射线管(CRT)监视器）、字母数字输入设备(例如，键盘）、光标控制设备(例如，鼠标或触 控板）、磁盘驱动器单元、信号生成设备(例如，扬声器）、触摸屏、加速计、麦克风、蜂窝无线 电频率天线，W及网络接口设备，网络接口设备可W是例如网络接口卡(NIC)、Wi-Fi卡或蜂 窝调制解调器。
[0213] ^引用的方式并入
[0214] 贯穿本公开已经参考并且引用了其它文档，诸如专利、专利申请、专利公开案、期 刊、书籍、论文、网络内容。所有此类文档在此W全文引用的方式并入本文中用于所有目的。 脚引等效物
[0216]除本文示出且描述的那些之外，本领域的技术人员将从本文档的完整内容(包含 对在本文中引用的科学和专利文献的参考）中明白本发明的各种改良及本发明的许多进一 步的实施例。本文中的主题含有重要信息、范例和指南，它们可适于本发明在本发明的各种 实施例及其等效物中的实践。
【主权项】
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所述的系统，其中所述数据结构包括表示远程基因特征数据库中 的每个经注释的特征的节点。
【文档编号】G06F19/22GK105830078SQ201480067827
【公开日】2016年8月3日
【申请日】2014年10月15日
【发明人】丹尼斯·库拉尔, 南森·梅维斯
【申请人】七桥基因公司