专利名称:用于关于参考多核苷酸序列标注样本多核苷酸序列中的变异的方法和系统的制作方法
用于关于参考多核苷酸序列标注样本多核苷酸序列中的变
异的方法和系统
背景技术:
在以下讨论中将描述某些论文和方法用于背景和介绍性目的。在此包含的任何事物不解释为承认现有技术。申请人明确地保留在适合的情况中、按照可应用的法定规定展示在此引用的论文和方法不构成现有技术的权利。基因研究近年来已经看到快速的进展。已经测序包括一些个体人类的特定有机体的完整基因组,并且其作为参考变得可用,即用作用于研究相同物种的成员的参考的基因序列。在基因研究、基因测试、个体化医药和多个其它应用中,通常有用的是获取基因材料的样本、确定该样本的序列、以及关于一个或多个现有参考分析该样本以鉴定序列变异或获取关于样本的其它有用信息。基因测试的现有方法通常定位或映射长、邻接样本序列到参考中的位置。然而,一些用于获得样本序列的技术产生包括具有预定空间关系的多个较短序列(有时称为寡聚体)的多核苷酸序列,在一些情况下包括具有可变距离的空间关系的多个较短序列的多核苷酸序列。在后者的情况下,在这些较短序列中碱基的相对基因组位置仅是近似已知,并且通常是以具有可变但受限数量的间隔或重叠(称为缺口距离)的短邻接读出(read)的形式。对于序列组装可用的传统技术不足以提供高速低成本的、包括可变有缺口邻接读出的短序列的从头(de novo)组装或重组装。因此,存在对于用于变异标注和组装的改进方法和系统的需要。本发明解决该需要。发明概述提供该发明概述从而以简化的形式介绍以下在详述中进一步描述的概念的节选。 该发明概述不旨在鉴定请求保护的主题的关键或实质特征,也不旨在用于限制请求保护的主题的范围。请求保护的主题的其它特征、细节、功用和优点将从以下撰写的、包括所附权利要求中定义的那些方面的详述变得显然。示例性实施方案提供用于与参考参考多核苷酸序列相比,标注样本多核苷酸序列中的变异的方法和系统。映射配对读出可包括与参考比对的两臂,或者与参考比对的一臂。 实施方案的各方面包括在至少一个计算机上执行定位参考多核苷酸序列中的局部区域的应用,该局部区域存在样本多核苷酸序列的一个或多个碱基从参考多核苷酸序列中的对应碱基变化的似然,其中至少部分基于样本多核苷酸序列的映射配对读出确定该似然;对每个局部区域产生至少一个序列假设,并且对局部区域的至少一部分优化该至少一个序列假设,以对局部区域得到高概率的一个或多个优化的序列假设;以及分析优化的序列假设以鉴定样本多核苷酸序列中的一系列变异标注(call)。根据示例性实施方案,在样本多核苷酸序列中进行的变异标注受到数据的有力支持。根据某些示例性实施方案,统计分析可以用于变异标注和/或用于很大程度基于已知参考的邻接序列的组装,但是邻接序列包括不同于参考的序列的改变和变异,这样的变化包括但不限于一个或多个碱基的删除、插入、突变、多态性和重复或重排。尽管可以描述示例性实施方案主要用于包括可变有缺口读出的配合读出。但是示例性实施方案也可以配置用于与具有其它预测或限定的空间关系的配合读出(例如,具有无缺口读出的配合读出或配合读出的邻近端是无缺口的情况)一起使用。附图简述
图1是图示根据一个示例性实施方案用于标注样本多核苷酸序列中的变异的系统的框图。图2是图示一个示例性配合读出的框图。图3是图示映射配合读出(即,映射到参考多核苷酸序列中的位置的配合读出) 的框图。图4A是图示根据一个示例性实施方案用于标注与参考多核苷酸序列相比从样本多核苷酸序列获得的映射配合读出中变异的处理的图。图4B是图形地描绘对每个局部区域生成一个或多个序列假设的集合的框图。图4C是图形地描绘对每个局部区域300优化序列假设412的框图。图4D是图形地图示变异标注的图。图5是图示根据一个示例性实施方案由用于标注样本多核苷酸序列中的变异的变异标注器执行的处理的细节的流程图。图6是图示使用贝叶斯(Bayesian)公式在计算参考分期间执行的处理的流程图。图7A和7B是图示用于使用部分德布鲁因(de Bruijn)图计算局部从头区间的处理的图。图8是图示根据一个示例性实施方案用于得到优化区间的处理的流程图。图9A和9B是图示用于使用贝叶斯公式和德布鲁因图执行优化的处理的流程图。发明详述除非另外指明,在此描述的技术的实践可以采用有机化学、聚合物技术、分子生物学(包括重组技术)、细胞生物学、生物化学和测序技术的传统技术和描述,这些在本领域的技术人员的范围内。这样的传统技术包括聚合物阵列合成、多核苷酸的杂交和连接、以及使用标记物的杂交的检测。通过参考这里的示例可以获得适当的技术的具体例示。然而,当然也可以使用其它等效的传统过程。这样的传统技术和描述可以在标准实验室手册中找到,诸如 Green 等编、1999 年出版的 Genome Analysis :A Laboratory Manual Series (Vols. I—IV),Weiner、Gabriel、St印hens 编的 0007 年)的 Genetic Variation :A Laboratory Manual, Dieffenbach、Dveksler 编(2003 年)的 PCR Primer :A Laboratory Manual, Bowtell 禾口 Sambrook(2003 年)的 DNA Microarrays :A Molecular Cloning Manual, Mount (2004 年)的 Bioinformatics :Sequence and Genome Analysis, Sambrook 禾口 Russell (2006 年)白勺 Condensed Protocols from Molecular Cloning :A Laboratory Manual, VX R Sambrook 禾口尺ussell(2002 ^ )的 Molecular Cloning :A Laboratory Manual (所有来自 Cold Spring Harbor Laboratory 出版社),Stryer,L. (1995 年)的、纽约州纽约的W. H. Freeman出版的Biochemistry (第4版),Gait的、伦敦的IRL出版社1984 年出版的“Oligonucleotide Synthesis :A Practical Approach”,Nelson禾口Cox(2000年) 的、纽约州纽约的 W. H. Freeman 出版的 Lehninger,Principles of Biochemistry (第 3版), 以及Berg等(2002年)的、纽约州纽约的W. H. Freeman出版的Biochemistry (第5版),所有这些通过引用将其全文合并在此用于所有的目的。注意,如在此和所附权利要求中使用的,单数形式“一个/ 一种”和“这个/这种” 包括复数指代物,除非上下文另外清楚地指示。因此,例如,谈及“(一个)配合读出(read),, 是指一个或多个多核苷酸配合读出,并且谈及“概率分析”包括提及本领域技术人员已知的等效步骤和方法等。除非另外定义,在此使用的所有技术和科学术语具有与实施方案所属技术领域的普通技术人员通常理解的相同的含义。在此提及的所有公开通过引用结合,用于描述和公开可以与目前描述的实施方案有关地使用的设备、公式和方法的目的。在提供值的范围的地方,可理解该范围的上限和下限和该声明范围中的任何其它声明的或居间的值之间的每个居间值包括在实施方案之中。这些较小范围的上限和下限可以独立地包含在较小的范围内,并且也可以包含在实施方案之中,服从声明的范围中具体排除的任何极限。常数值仅用于图示,并且用于执行示例性实施方案的系统可以配置为与其它值一起工作。定义以下定义有助于提供用于理解本发明的背景。在此使用的“多核苷酸”、“核酸”、“寡核苷酸”、“寡”或语法等同物通常是指以线性方式共价连接到一起的至少两个核苷酸。核酸通常包含磷酸二酯键、尽管在一些情况下可以包括具有可替代主链,诸如亚磷酰胺(phosphoramidite)、二硫代磷酸酯 (phosphorodithioate)、或甲基亚磷酰胺(methylphophoroamidite)连接,或者肽核酸主链和连接的核酸类似物。其它类似物核酸包括具有双环结构的核酸,包括锁定核酸、正主链、非离子主链和非核糖主链。术语“参考多核苷酸序列”或者简单的“参考”是指参考有机体的已知核苷酸序列。 参考可以是参考有机体的完整基因组序列、参考基因组的部分、多个参考有机体的共有序列、基于不同有机体的不同组分的汇编序列、从有机体群体提取的基因组序列的集合、或者任何其它合适的序列。参考还可以包括关于已知要在有机体群体中得到的参考的变异的信息。参考有机体还可以是对分别提取的正在测序的样本特异性的,可能是有关个体或相同个体(可能对补充肿瘤序列正常的)。“样本多核苷酸序列”是指自基因,调节元件,基因组DNA,cDNA,包括mRNA、rRNA、 siRNA.miRNA的RNA等及其片段衍生的样本或目标有机体的核酸序列。样本多核苷酸序列可以是来自样本的核酸、或者诸如扩增反应的产物的二代核酸。对于样本多核苷酸序列、或要从样本多核苷酸(或任何多核苷酸)“衍生”的多核苷酸片段,可以意味着通过物理地、 化学地、和/或酶促的样本多核苷酸(或者任何多核苷酸)片段化形成样本序列/多核苷酸片段。要从多核苷酸“衍生”还可以意味着片段是源多核苷酸的核苷酸序列的特定子集的复制或扩增的结果。“配合读出(mated read) ”通常是指源自相距几百个或几千个碱基定位的基因组序列的两个不同区域(臂)的一组个体核苷酸读出。在从由要变异标注(call)和/或重组装(reassemble)的样本有机体得到的较大邻接多核苷酸(例如,DNA)的片段测序期间, 可以生成配合读出。“映射(map) ”是指例如通过匹配举例的配合读出与对应于参考中的位置的索引中的一个或多个码(key),使配合读出与同配合读出类似的参考中的零个、一个或多个位置相关的处理。“贝叶斯(Bayes)定理”是指用于计算条件概率的数学公式,其中在数据的给定体 E的条件下的假设H的概率是假设和数据的联合的非条件概率与单独数据的非条件概率的比。以E为条件的H的概率定义为(H) = P (H&E) /P (E)。“德布鲁因(De bruijn)图”是指其顶点或节点是来自某一字母表的符号的序列, 并且其边缘指示可能重叠的序列的图。“部分德布鲁因图处理”是指不计算完整的德布鲁因图,但是基本遵循顶点的初始图(优选从参考得出)、并且确定是否存在分支的处理。可以使用超过某一阈值的任何分支,以指示在局部区域中存在变异的似然。“重组装”和“重测序(resequence) ”是指对照参考组装映射配合读出,以建立类似于原始参考但不必与其相同的序列的方法。这不同于其中配合读出组装在一起以形成新的之前未知的序列的从头组装。在以下描述中,阐述许多具体细节以提供本实施方案的更彻底理解。然而,实践本实施方案不需一个或多个这些具体细节将对本领域的技术人员显而易见。在其它示例中, 不再描述公知特征和对本领域技术人员公知的过程,以便避免实施方案不明显。图1是图示根据一个示例性实施方案用于标注样本多核苷酸序列中的变异的系统的框图。在该实施方案中,系统100可以包括1-n个计算机12的计算机集群、以及数据存储库14。例如,在一个特定实施方案中,系统可以包括32个计算机。计算机12可以经由高速局域网(LAN)切换结构16连接到数据存储库14。至少部分计算机12可以并行执行变异标注器(caller) 18的实例。变异标注器18可以包括贝叶斯公式20组件和德布鲁因图 22组件。数据存储库14可以存储包括一个或多个数据库的若干数据库,该一个或多个数据库存储参考多核苷酸序列M、在生化处理期间从样本多核苷酸序列获得的配合读出观、 以及从配合读出观生成的映射配合读出观。参考多核苷酸序列24 (以下简称为参考)是指参考有机体的已知核苷酸序列(例如,已知基因组)。这包括包含具有在基因组内的一个或多个位置的已知变异的序列的参考。多核苷酸分子是由链中共价键合的核苷酸单体构成的有机聚合物分子。脱氧核糖核酸 (DNA)和核糖核酸(RNA)是具有不同生物功能的多核苷酸的示例。有机体的基因组是编码为DNA或RNA的有机体的遗传信息的总和。单倍体基因组包含有机体的每个遗传单元的一个拷贝。在诸如哺乳动物的二倍体动物中,基因组是一系列互补多核苷酸,其包括组织为具有离散基因单元或等位基因的染色体的集合的大部分遗传信息的两个拷贝。等位基因的每个拷贝提供在单个染色体上的具体位置,并且基因组中每个等位基因的基因型包括确定具体特征或性状的同源染色体上的特定位置处存在的等位基因对。基因组包括等位基因的两个相同拷贝的情况中,这对该等位基因是纯合的,并且当基因组包括两个不同的等位基因时,这对该基因座是杂合的。DNA本身组织为互补的两条多核苷酸链。在以下描述中,对应于参考的链可以称为“正”,并且与参考互补的链可以称为“负”。参考M可以是完整的基因组序列、参考基因组的部分、多个参考有机体的共有序列、基于不同有机体的不同组分的汇编序列、或任何其它合适的序列。序列M还可以包括关于已知要在有机体群体中得到的参考的变异的信息。在对要分析的样本有机体的多核苷酸序列(例如,来自基因、基因组DNA、RNA或其片段的核酸序列)执行测序处理期间,可以获得配合读出26。配合读出沈可以是得自包含完整基因组,诸如完整哺乳动物基因组,更特别的完整人基因组的样本。在另一实施方案中,配合读出沈可以是来自完整基因组的具体片段。在一实施方案中,通过对扩增的核酸结构(诸如,使用聚合酶链反应(PCR)或滚环复制创建的扩增物(amplimer))执行测序,可以获得配合读出26。例如在美国专利公布No. 20090111705、No. 20090111706和 No. 20090075343中描述可以使用的扩增物的示例,通过引用并入其全部内容。映射配合读出28是指已经映射到参考M中的位置的配合读出26,如以下映射部分中进一步说明的。在一实施方案中,变异标注器18包括程序指令,其至少部分基于贝叶斯公式20执行的证据推理、以及德布鲁因图22组件和部分德布鲁因图22’组件二者执行的基于德布鲁因图的算法的组合,对参考M和映射配合读出观执行诸如概率分析或准许计分可替代假设的其它方法的统计分析。使用统计分析,变异标注器18为了鉴定和标注与参考M相比在映射配合读出观的序列中检测到的变异或差别的目的,生成和计分序列假设。在一实施方案中,变异标注器 18能够鉴定和标注包括一个或多个碱基的删除、插入、突变、多态性和重复或重排等的序列的大规模结构变异。变异标注器18可以输出包含鉴定的变异的变异标注32文件、列表或其它数据结构,每个描述观察映射配合读出观的序列的部分在或靠近具体位置区别于参考M的方式。 在一实施方案中,变异标注器18也可以输出由于计算的不确定性不能标注的无标注区间的列表。在另一实施方案中,变异标注器18可以配置为使用概率分析或计分可替代假设的其它方法,以重组装或重测序来自映射配合读出26的样本有机体的多核苷酸序列,其中组装的多核苷酸序列实质上是基于参考对,但是包括鉴定的变异。在一实施方案中,变异标注器18可以实施为链接到库(未示出)以执行具体任务的单个可执行。在另一实施方案中,变异标注器18可以实施为相互通信的独立应用模块。 此外,变异标注器18的功能可以跨比图1所示数量多或少的软件模块/组件分布。为了加速处理,可以配置计算机集群10,使得在不同的计算机12上执行的变异标注器18的示例并行地对参考M和映射配合读出沈的不同部分进行操作。任务调度器30 负责在计算机集群10中跨各种计算机12分配任务或工作包。计算机12可以包括典型的硬件组件(未示出),包括一个或多个处理器、输入设备 (例如,键盘、指向设备等)、以及输出设备(例如,显示设备等)。计算机12可以包括计算机可读/可写介质,例如包含当由处理器执行时实施公开的功能的计算机指令的存储器和存储器设备(例如,快闪存储器、硬盘驱动器、光盘驱动器、磁盘驱动器等)。计算机12可以还包括用于实施数据存储库14和用于存储变异标注32的计算机可写介质。计算机12可以还包括用于通信的有线或无线网络通信接口。在一实施方案中,计算机集群10可以实施为以另外的处理器水平衡量的商品Intel/Linux集群。数据生成
在一些实施中,生化技术可以用于生成从要分析的有机体的样本多核苷酸获得的配合读出沈。在一实施方案中,生化技术提供在离散但是有关的集合中的数据,使得配合读出26的内容可以包括预测的空间关系和/或分离变异。可以基于关于用于生成配合读出 26的生化处理的现有知识(即,基于如果生化处理应用于样本期望获得的序列)、根据配合读出26的序列数据或其子集的初步分析的经验估计、专家的估计、或其它合适的技术,确定关系。许多生化处理可以用于生成供当前组装方法使用的配合读出沈。这些包括但不限于如美国专利No. 6,864,052、No. 6,309,824,6, 401,267中公开的杂交方法,如美国专利 No. 6,210,891、No. 6,828,100、No. 6,833,246、No. 6,911,345、No. 7,329,496 和 Margulies 等 2005 年在 Nature 杂志 437 期 376-380 页和 Ronaghi 等 1996 年在 Anal. Biochem 杂志242期84-89页中公开的边合成边测序(sequencing-by-synthesis)方法,如美国专利No. 6,306,597、W02006073504、W02007120208中公开的基于连接的方法,如美国专利 No. 5,795,782、No. 6,015,714、No. 6,627,067、No. 7,238,485 和 No. 7,258,838 以及美国专利申请No. 2006003171和No. 20090029477中公开的纳米孔测序技术,以及如美国专利申请No. 20090111115中公开的纳米通道测序技术,所有这些通过引用并入其全文。在一具体的实施中,在一些实施方案中使用复合探针锚连接(Combinatorial Probe Anchor Ligation, cPAL)处理(参见美国专利申请公开No. 200802;34136和No. 20070099208,在此通过引用并入其全文)。图2是图示一个示例性配合读出200的框图。在一实施方案中,在测序期间从较大邻接多核苷酸(例如,DNA)的片段生成配合读出200,较大邻接多核苷酸从要变异标注和 /重组装的样本有机体获得。仅原始片段的部分碱基在每个配合读出200中读取。产生配合读出200的生化操作之后,来自片段的碱基称为的碱基读出202。这些来自片段的碱基读出202由此组成原始较大多核苷酸中的碱基的小部分。单个测序实验可以从样本多核苷酸序列(例如,样本基因组G)生成Nd配合读出 26,并且每个配合读出200可以仅包括原始片段的nD碱基,称为碱基读出202。这导致总数 nDND个碱基。假设存在样本基因组G的Ne个碱基,每等位基因的平均覆盖是c = nDND/Ne。 这意味着忽略覆盖偏倚,样本基因组G的Ne个碱基的每个(在二倍体区域中计数二个等位基因)接收碱基读出的平均数c。在具体的测序实验中,c 25,并且原始片段的碱基读出的数量nD = 70,因此配合读出200的数量Nd = 2xl09o在其最简单的形式中,每个单个碱基读出202分别包括与DNA或RNA的核苷酸中得到的4个碱基对应的四符号值A,C,G,T (或U)之一。关于一些碱基读出202的碱基值可以偶然地丢失,在这样的情况下碱基值由符号N代表,并且称为无标注。在一实施方案中, 每个碱基读出202与指示碱基读出202的估计质量的质量分214 (在例如0和1、或0和99 之间的数)关联。变异标注器18可以使用质量分确定来自读出204中单个碱基读出202 的证据贡献。在一实施方案中,在生化处理期间读取的碱基读出202的相对位置仅是近似地已知,并且以md邻接读出204的形式从0到md-l编号。用于生成配合读出沈的生化处理可以导致各种长度的读出204,包括单个配合读出200中的不同长度。第i个读出204具有Ii 碱基读出202的长度,在一实施方案中Ii可以是对给定测序项目中的所有有关的配合读出200相同。该限制以某一另外的记法复杂性为代价,可以可选地放松。在一实施方案中,配合读出200包括数据中的多个缺口。例如,md读出204可以由标号从0到md-2的md-l个缺口 206 (gl到g5)分隔。关于具体配合读出200的每个缺口 206的值可以是可变的和未知的。然而缺口值的统计分布是可得到的,先验的或来自通过映射的经验采样。应当注意,图2所示配合读出200仅用于说明的目的,以及不同的库(S卩,读出)架构可以具有变化数量的读出204,并且相应地读出之间的读出内缺口可以具有可变的长度, 以及缺口可以出现在任何读出之间。在一实施方案中,配合读出200可以包括将配合读出200划分为称为左臂2IOA和右臂210B(统称为臂210)的两个子集要素的配对缺口 208。配对缺口 208可以包含相对大数量的值,诸如例如峰值在约400-500个碱基的分布和100个碱基的宽度。在一实施方案中,左臂210A和右臂210B可以分别包括md,1204。在一非常具体的实施方案中,生化处理产生每个具有35个碱基读出202的左和右臂210,并且在该实施方案中,左和右臂210每个可以包括4个邻接读出204011^ = md,2 = 4)。继而4个邻接读出204每个包括三个10碱基读出和一个5碱基读出。配合读出200 中的5碱基读出可能可以是相对于配对缺口 108的最远的一个(例如,rl和r8)。臂内缺口 206每个可以仅具有小数量的值,其是小的正或负整数。负缺口值对应于单个配合读出200中读取多于一次的相同基因组碱基。在另一实施方案中,可以配合多于两个臂,这可以是合适的,如果存在可以采用大范围值的另外的缺口的话。在一实施方案中,一特定配合读出200的md-l个缺口值g可以分割为三组,分别为 Hidjl-IU和md,2-l的gpgjP &,对应于左臂210A、配对缺口和右臂210B。关于三组缺口值的概率分布假设为不相关,因此关于缺口值g的概率分布可以写为乘积P(g) = P(g1)P(gm)P(g2)尽管整体配合读出200的分布可以写为乘积,但是第一和第三缺口值概率分布 p(gl)和p(g2)不必需解耦为关于各个缺口值的概率的乘积。此外,这是基于缺口概率独立于其它附近序列的假设。阅读本公开时用某一另外的标记复杂性完成消除该假设对本领域的技术人员是显而易见的。基因组模型在整体基因组的重测序组装中,假设测序的样本(例如,基因组G)非常类似于已知的参考对(例如,基因组&)。参考基因组(^1可以提供为毗连群(contigs,来自 contiguous)的集合,其是自基因源衍生的非重叠多核苷酸段的集合。在更复杂的基因组 (例如,哺乳动物的基因组)中,这些毗连群可以对应于染色体的子集。然而,因为所有的毗连群比配合读出沈源自的每个片段的长度长很多,为了该讨论的目的将好像包括碱基的单个邻接序列一样处理(V Gchi或者(^(i)应当用于指示中在位置i处的碱基。的每个部分具有关联的倍数性p,其是基因组中存在的特定基因单元的等位基因的数量。倍数性在二倍体有机体的基因组的多数区域中是2,包括人类基因。更大或更小的倍数性值在正常基因组也是可能的,例如,人类的多数“Y”染色体是单倍体,并且包括来自癌细胞或来自具有基因异常性的个体的基因组的异常基因组通常具有带有不常见的倍数性的区域。基因组部分的P等位基因不必彼此相等,而是全部非常类似于对应的的部分。在具体的方面,为了分析的缘故处理每个区域的倍数性好像其预先已知,尽管这样的假设对于所述实施方案的方法不是严格必需的。术语用于指示遍及所有等位基因合计的基因组中碱基的总数,其用于重测序数据的组装。对于正常的二倍体人类基因组,该数字约为6X109。配合读出概率模型为了变异标注和重组装中使用统计计算的目的,如下建模配合读出200的产生过程。每个配合读出200起源于正在测序的样本基因组G的两个等位基因之一上的一随机位置,两个互补链之一上。由于覆盖偏倚,起源基因组位置、链和等位基因不是相等地可能的。配合读出200其第一碱基在等位基因a的链s的给定位置χ的概率采用P(x, s, a) ——b(xfsf a)
2ng其中b(x,s, a)是具有平均数1的未知覆盖偏倚。在以下多数描述中以及在变异标注过程中,假设b(x,s,a) = 1,并且不考虑覆盖偏倚。在该情况下,P(X,s,a)是常数,且 P(x, s,a) = 1八2队),并且等于基因组的DNA的两条链中的碱基的总数的倒数。碱基读出碱基读出202的质量分214例如使用配合读出200或其要素映射到参考M期间获得的统计(以下说明)、或直接从图像分析得出的估计,可以转换为估计的误差率Q。作为最简单的可能假设,来自参考M的碱基(例如,基因组碱基)以概率1-6’等于配合读出 200中的碱基读出202,并且以概率6’/3等于三个剩余读出的任一。这部分基于单个碱基读出202彼此不相关的暗含假设。如从本公开的方法显然的,在另一实施方案中有利于不同的误差模型可以撤销该假设。更正式地,如果bg是真实的正在读取的基因组碱基,并且bd是作为配合读出200 中读取的碱基,以前者为条件的后者的概率是由下式给出的4X4误差矩阵P(bd\bg) = (1- 6)S(bg,bd) + [1 - S(bg,bd)]6/3其中如果克罗内克(Kronecker) δ符号的自变量相同,则其是1,否则是0。更复杂的误差模型是可能的,其中每个碱基读出带有作为P (bg |bd)的估计的4个权数。以下可以进一步扩展贝叶斯公式,以考虑这样的估计。在另一实施方案中,可以以在序列的两侧为条件建模误差概率,或者以其它碱基是否有误为条件的误差概率,等等。映射映射配合读出观可以从配合读出沈的数据集合创建。高速映射软件(例如,映射到参考(map-to-reference)组件(未示出))可以用于将原始配合读出沈映射到配合读出沈可能出现在的参考M中的碱基位置。映射到参考组件可以输出以参考M中配合读出沈的位置的形式的映射配合读出观。映射可以容忍来自参考M的小变异,诸如读出 202中个别基因组变异、读出误差或未知碱基导致的变异。图3是图示映射配合读出28的框图。配合读出是已经映射到参考M中的位置的配合读出。在一实施方案中,映射配合读出观可以包括仅具有与参考比对的一臂210的配合读出,即半映射。更具体地,当对于要通过映射处理(诸如以下段中列出的专利申请中描述的)识别的相似性,臂210的读出204足够地类似于在一位置和链开始的参考M的序列时,考虑配合读出以映射到参考M中的该位置和链。从这样映射可以推断配合读出可以从关注的样本的对应位置衍生。在另一实施方案中,映射配合读出观可以包括具有与参考上的规定位置比对的两臂210的配合读出,即全映射。为了支持包括大规模结果改变或密集变异的区域的较大变异的组装,可以单独映射配合读出26的左和右臂210,在比对之后估计配对约束。在一实施方案中,映射到参考组件可以执行如以下临时专利申请中公开的映射处理2009年2月3日提交的、申请号为61/149,670的OLIGOMER SEQUENCES MAPPING (CGI002PRV), 2009 年2 月 3 日提交的、申请号为 61/149,665 的 INDEXING A REFERENCE SEQUENCE FOR OLIGOMER SEQUENCE MAPPING (CGI003PRV),以及 2009 年 2 月 3 日提交的、申请号为61/149,689 的OLIGOMER SEQUENCES MAPPING(CGI004PRV),所有这些申请转让给本申请的受让人,并且在此通过引用并入其全文。总结该处理,配对臂读出可以使用二阶段处理比对到参考对。首先,左和右臂210 可以使用参考基因组的索引独立地比对。在一实施方案中,该初始搜索可以得到至多以两个单碱基取代匹配左臂210A的参考M中的所有位置,但是可以得到具有达到五个错配的一些位置。其次,对于在第一阶段中鉴定的左臂的每个位置,右臂210B经历可以约束为由配合距离(例如,离开0到700个碱基)的分布告知的基因组区间的局部比对处理。随后可以使用臂210B的位置,对臂210A重复相同的局部比对处理。在两个阶段,可以通过尝试缺口值的多个组合执行有缺口臂读出210的比对。如果配合读出200具有臂210的任何一致位置(即,左和右臂按适合的顺序并在期望配合距离分布中在相同的链上),则只有这些位置可以保留。否则,可以保留臂210的所有位置。在任一情况下,出于性能的原因,在一实施方案中,对每个臂210报告的位置的数量可以限制为例如50个位置。在一实施方案中,对映射配合读出观的存储记录可以包括对按参考顺序分类的参考多核苷酸序列(Gtl) 24的臂映射列表,以鉴定可以可能对L(G)的贝叶斯计算有贡献的序列,以下进一步说明。该记录可以不仅包含映射信息,还可以包含对应的配合读出200的碱基读出204和质量分214。变异标注处理图4是图示根据一个示例性实施方案用于关于参考多核苷酸序列M标注样本多核苷酸序列中的变异的处理的图。该处理可以包括在至少一计算机12上执行定位参考多核苷酸序列中的局部区域的应用,该局部区域中存在样本多核苷酸序列的一个或多个碱基从参考多核苷酸序列中的对应碱基变化的似然,其中至少部分基于样本多核苷酸序列的映射配合读出确定该似然(块400)。通过横跨参考M和配合读出28中的对应位置的垂直条,图3示出局部区域300 的示例。该局部区域300鉴定映射配合读出观有可能具有来自参考多核苷酸序列M中的对应碱基的一个或多个碱基中的变异或变化的位置。局部区域300可以包括一个或一序列的碱基位置。在一实施方案中,至少部分在计算参考分过程期间使用贝叶斯公式20组件、以及使用部分德布鲁因图22’组件以计算局部从头区间,由变异标注器18执行定位有可能具有来自参考的变化的局部区域300,如以下参考图5进一步所述。在一实施方案中,可能变化的小或单独局部区域300可以组合以在得到优化区间过程期间形成称作优化区间的较大局部区域,如图5所述。变异标注器18对每个局部区域300生成至少一个序列假设,并且对至少部分局部区域优化该至少一个序列假设以对该局部区域得到高概率的一个或多个优化的序列假设 (块40幻。在一实施方案中,对局部区域生成序列假设可以包括遍历局部区域中的每个碱基位置和用每个可能的可替代碱基值迭代地改变碱基,包括在该位置的插入或删除,其可以导致对每个局部区域300生成的多个序列假设。图4B是图形地描绘对每个局部区域生成一个或多个序列假设的集合的框图。在该示例中,已经在参考多核苷酸序列M中鉴定可能变化的两个局部区域300A和300B (统称为局部区域300)。对局部区域300A和300B的每个,通过在各自的局部区域300中改变碱基值,分别生成序列假设412A和412B (统称为序列假设41 的集合。图4C是图形地描绘对每个局部区域300优化序列假设412的框图。在一实施方案中,优化序列假设412可以包括对其中每个单独序列假设计算概率比。对局部区域300的每个导致高概率比的序列假设412可以用于形成优化的序列假设414A和414B(统称为优化的序列假设414)的集合。在一实施方案中,来自最大化概率比的优化序列假设414A的序列假设应用于局部区域300A。根据一实施方案,在利用贝叶斯公式20组件和德布鲁因图组件22的序列假设产生和优化过程期间,执行优化,如参考图5进一步说明。变异标注器18分析优化的序列假设的每个集合,以鉴定样本多核苷酸序列中的一系列变异标注,其中变异标注通过映射配合读出有力地支持004)。在一实施方案中,可以在提取标注处理期间鉴定一系列变异标注,该处理可以包括指示在或靠近具体位置与参考M相比样本多核苷酸序列沈的映射配合读出观的碱基中检测的变异。在一具体实施方案中,可以至少部分基于统计概率分析鉴定变异标注。图4D是图形地图示变异标注的图。在该示例中,单独检查优化的序列假设414的全部或子集,并且关于最正确地描述参考M中的对应序列的优化假设的一或多部分,做出变异标注32。在一实施方案中,由提取标注组件做出变异标注,如图5中进一步说明。图5是图示根据一个示例性实施方案由用于标注样本多核苷酸序列中的变异的变异标注器18执行的处理的细节的流程图。变异标注器18可以包括五个主要处理使用贝叶斯公式计算参考分500,使用部分德布鲁因图计算局部从头区间501,得到优化区间502, 使用贝叶斯公式和德布鲁因图生成假设并优化504,以及提取标注506。使用这些处理,变异标注器18可以通过迭代地最大化基因组的后验概率P (G/MtdRds),获得重组装的样本多核苷酸序列,例如基因组(G),该后验概率解释所有映射配合读出(MtdRdsUS。该处理可以通过变异标注器18接收参考多核苷酸序列(Gtl) 24和映射配合读出观作为对全部五个处理的只读输入开始。变异标注器18也可以接收作为只读输入的元数据 508,其可以描述任务,即鉴定参考M、输入文件的位置、定义映射配合读出观的结构等。如前所述,参考对包括参考有机体的已知核苷酸序列,例如国家生物技术信息中心(NCBI)建立的参考人类基因组36。映射配合读出观每个包括由可变缺口分隔的多个邻接读出204。 构成映射配合读出观的读出204的每个包括位于映射配合读出中各个碱基位置处的多个碱基。也将这些碱基映射到参考M中的碱基位置。使用该输入,使用贝叶斯公式计算参考分500处理计算参考分510,并且使用部分德布鲁因图计算局部从头区间501处理计算局部从头区间512。参考分510在一实施方案中,使用贝叶斯公式计算参考分500(以下称为计算参考分500)处理有助于鉴定可能变化的局部区域300。该处理对参考M的每个碱基计算在该碱基位置的所有可能的1碱基变异与在该碱基位置的参考残基之间的概率比。在一实施方案中,每个 1碱基变异可以称为初始假设。图6是图示在使用贝叶斯公式计算参考分期间执行的处理的流程图。计算参考分 500处理可以包括通过以所有可能的1碱基变异修改ρ等位基因中在该位置的碱基值,对参考中的每个碱基位置生成初始假设的集合(块600)。修改可以包括单个碱基变化(以反映序列中潜在的单核苷酸多态性(SNP))、以及插入和删除(统称为插入和删除)。除了单碱基变化之外,当在包括串联重复的基因组区域进行修改时,初始假设可以包括达到10 个碱基的插入和删除,例如在三联体碱基串联重复的区域中可以插入或删除三个碱基。计算参考分500处理确定接近(并可能重叠)参考的当前碱基位置的映射配合读出的集合(块60幻。在一实施方案中,鉴定映射配合读出28,其具有映射到近似地在合适方向和链远离当前碱基位置的配对缺口期望长度的位置的一臂,以满足期望的配对关系。贝叶斯公式20用于通过对对应的集合中的每个初始假设计算概率比Pv/PKef,对参考的每个碱基位置计算参考分,其中Pv是1碱基初始假设的概率,且Pltef是参考中碱基值的概率,并且其中靠近每个碱基位置的映射配合读出的集合在计算每个碱基位置的概率比期间使用(块604)。在又一实施方案中,假设样本多核苷酸序列包括基因组G,则每个参考分可以包括
对数似然比L(G),其中L(G) = l0g(Pv/Pltef)。在又一实施方案中,对数似然比可以表达为
(以下更充分地描述)
Γ π "“、 , P(GlMtdRds)0107 L(G) = loa—~---—.
、J y P(G0\MtdRds)在多数情况下,L(G)计算为大和负的,表达在该位置不存在来自的变异的事实。 然而,在存在1碱基变异的位置,L(G)将计算为正数。该计算将足以标注和组装1碱基变异,但更长的变异的似然使情况复杂。在存在更长的变异的局部区域300,对1碱基变异的 L(G)可以仍是正的,但是比存在1碱基变异的局部区域中的值低很多。L(G)中的低值正数可以用于识别变异存在的似然,并且鉴定局部区域300。其中一个或多个概率比超过阈值的参考M中的位置可以输出为有可能已经相对于参考改变的局部区域300。计算参考分500的输出是包含计算参考分510的表格的文件。该表格的每一行可以引用参考基因组的位置,并且包含关于该位置的参考分,其中参考分是在该位置纯合或杂合形式的所有可能的1碱基变异假设的计算概率。在一实施方案中,对每个可能的假设, L(G) = log(Pv/PEef)可以存储为参考分。在一实施方案中,参考分510的表格的每一行可以包含鉴别关于该行的参考基因组的位置的字段。这可以包括染色体标识和位置,或等效信息;三个字段,包含关于在感兴趣的位置的所有三个可能的纯合SNP的L(G)=log(Pv/PEef);一个字段,包含关于在感兴趣的位置的纯合删除的L(G) = log(Pv/PEef);四个字段,包含关于紧接在感兴趣的位置之前的所有四个可能的纯合1碱基插入的 L (G) = log(Pv/PEef);以及复制之前八个字段的另外的八个字段,但是针对杂合变异而不是纯合变异。在一些实施方案中,报告具有任一评估的替代的最高参考分的初始假设的分作为关于给定位置的参考分是有用的,包括评分多碱基修饰(诸如关于由局部从头处理建议的序列或短串联重复的串联重复拷贝变化)的似然。局部从头区间512根据示例性实施方案,使用部分德布鲁因图计算局部从头区间501的处理(以下称作计算局部从头区间501)补充计算参考分500处理,构成用于定位有可能具有变化的样本多核苷酸序列中的局部区域300的另一部件。该处理可以通过得到映射到参考位置但是也延伸越过该位置的映射配合读出28,生成超过单个碱基变化的变异。部分德布鲁因图 501处理是指没有计算完整德布鲁因图的处理,但是基本上遵循参考顶点的初始图,并且确定是否存在分支。任何超过某一阈值的分支可以认为指示局部区域300中存在变异的似然。这样的局部区域300可以经历进一步的分析(例如,以下描述的优化),使得确定如果有存在什么变异。作为总的看法,局部从头处理用从碱基序列(例如每30个碱基长度)创建的参考顶点,从参考M初始化局部德布鲁因图。随后迭代的扩大该图。对图中已经存在的每个顶点,得到良好的映射到对应的参考顶点、但是包括以任何可能的1碱基值(A、C、G或T/U)延伸超过参考顶点的任一端的映射配合读出的集合。在一实施方案中,当至少读出的14个碱基匹配该序列,并且至多读出的1个碱基错配该序列时,可以认为将映射配合读出良好地映射到对应的碱基序列。在一实施方案中,该过程对每个碱基延伸计算延伸强度,延伸强度代表关于至少部分基于具有相同延伸的映射配合读出的数量、以及那些映射配合读出与正处理的顶点的序列匹配和错配的数量,以每个1碱基值延伸参考顶点的支持量。每个这样的1碱基延伸对应于图中的一个可能的顶点。增加具有最高延伸强度和还没有在图中的碱基延伸。在一实施方案中,用于构造部分德布鲁因图的图建立过程是递归的,但是也可以使用正确地重复必需的步骤的任何过程。在图建立过程的一实施方案中,在一个方向以深度优先的方式计算每个顶端碱基延伸的延伸强度,并且在对具有高于阈值的延伸强度的碱基延伸的每个计算之后,创建新边沿和新顶点。如果在路径中不存在具有高于阈值的的延伸强度的碱基延伸,则对该路径返回失败。如果图建立过程的计算创建等于参考顶点之一的碱基序列、以及与SNP或短插入删除(约50个碱基或更短)一致的新分支顶点,则计算结束并返回该路径。图7A和7B是图示用于使用部分德布鲁因图计算局部从头区间的处理的细节的图。该处理可以包括用参考顶点初始化德布鲁因图,其中每个参考顶点代表来自参考基因组的N个邻接碱基的碱基序列(块700)。在一实施方案中,例如,每个碱基序列可以重叠紧邻的碱基序列,例如共享差不多1个碱基。在一实施方案中,参考顶点可以包括例如30个碱基的碱基序列(N = 30)。对每个参考顶点,确定映射到参考顶点、并且包括以任何可能的碱基值(A、C、G或Τ/υ)在左或右方向延伸超过参考顶点的任一端的碱基延伸的映射读出的集合,并且对每个碱基延伸计算代表关于以每个1碱基值延伸参考顶点的支持量的碱基延伸强度,并且指示碱基延伸的左或右方向(块702)。在一实施方案中,可以至少部分基于具有相同延伸的映射配合读出的数量、以及那些映射配合读出与正处理的顶点的序列匹配和错配的数量,计算延伸强度。该处理可以将参考多核苷酸序列分割为参考段(块704),以将全部计算划分为较短的段。在一实施方案中,参考多核苷酸序列可以分割为50个碱基的参考重叠或相邻段。通过得到与参考段的每个中的参考顶点不相容、以及其碱基延伸强度大于第一阈值的达到M个顶端碱基延伸,得到来自参考的可能变异(块706)。在一实施方案中,可以使用顶端二(M = 2)碱基延伸。在一实施方案中,例如,每个碱基延伸可以越过参考顶点延伸 56个碱基。M个顶端碱基延伸可以用作部分德布鲁因图中的分支顶点,并且每个分支顶点和延伸方向输入到建立部分德布鲁因图的递归过程中(块708)。在深度优先递归处理中,确定映射到分支顶点、以及包括以每个可能的碱基值(Α、 C、G或T/U)在对应的方向延伸超过分支顶点的碱基延伸的映射读出的集合,并且对每个碱基延伸计算代表关于以每个碱基值延伸分支顶点的支持量的分支延伸强度(块710)。该处理从具有高于分支延伸阈值的分支延伸强度的碱基延伸创建新分支顶点 (块 712)。然后在分支顶点与每个新的分支顶点之间构造新的边沿,其中由一个或多个边沿连接的分支顶点的链形成通过该图的路径(块714)。对路径中的边沿计算边沿延伸强度,并且找出具有高于边沿强度阈值的最大边沿延伸强度的顶端N个边沿,以及使用顶端N个边沿的每个执行块722 (块716)。如果块722无法找出新分支顶点和参考顶点之间的匹配,则顶端N个新边沿输入到递归过程(710)作为分支顶点(块718)。如果没有边沿具有高于强度阈值的延伸强度,则结束递归过程,并对该路径返回失败(块720)。每个新分支顶点与代表来自参考的附近碱基序列的部分德布鲁因图中的参考顶点相比,并且如果新分支顶点等于参考顶点的碱基序列的任一,结束递归过程,并且返回该路径(块722)。部分德布鲁因图过程的输出可以是包含局部从头区间512的文件、列表或其它数据结构。局部从头区间512可以包括在该处部分德布鲁因图自参考发散分支的碱基位置的列表。在一实施方案中,计算参考分500处理和计算局部从头区间501处理可以同时并且并行地跨计算机集群10执行,以便得到映射配合读出的处理不需要重复。由各批计算的参考分510和局部从头区间512处理的结果可以合并为单个的参考分文件以及单个的局部从头区间文件。得到优化区间502以上处理得到有可能具有从参考变化的碱基的局部区域300。如下进一步处理这些局部区域300。在得到优化区间502处理期间,可以组合由参考分510和局部从头区间512代表的可能的变化的各个局部区域300,以形成较大的优化区间。在概念的水平,所得优化区间鉴定为可能包含变异的区域,并且如果有,要求进一步的分析以确定可能的变异。得到优化区间502处理考虑作为候选优化区间的局部从头区间512和与高概率比 (例如,高L(G))关联的参考分510。可以组合重叠或相隔小于阈值碱基距离(例如,小于臂210的长度)的候选优化区间,以形成优化区间。然而,一些优化区间可以变得过长,并且优化这样的优化区间会变得太昂贵。在一实施方案中,这些优化区间因此是先验的无标注。图8是图示根据一个示例性实施方案用于得到优化区间的处理的流程图。得到优化区间502处理可以包括分析与每个参考碱基位置关联的参考分,以及鉴定在该位置参考分高于预定数量阈值(即,高L(G))的连续位置的每个范围作为峰(块800)。在一实施方案中,用于高L(G)的数量阈值可以是例如-10db。对每个鉴定的峰,峰区间可以定义为由峰的碱基位置以及在该碱基位置的任一侧的预定数量的碱基形成的连续碱基位置的集合(块80幻。在一实施方案中,峰区间可以包括峰本身外加在每一侧的3个额外碱基。得到优化区间502处理可以从每个峰区间和在德布鲁因图的分支中得到的每个局部从头区间,形成各自的候选优化区间(块804)。得到优化区间502处理也可以创建重叠的候选优化区间(即,各个峰区间和/或局部从头区间)的各自并集(块806)。可以合并重叠或相隔小于规定距离的邻近候选优化区间(块808)。在一实施方案中,如果邻近候选优化区间相隔小于序列臂210的长度,合并它们。在一具体示例中,如果邻近候选优化区间在20个碱基的距离内,可以合并它们。无标注分配给比预定碱基长度长的候选优化区间,并且这样的候选优化区间增加到无标注区间列表(块810)。所有剩余的候选优化区间增加到优化区间列表(块812)。将无标注分配给候选优化区间的处理是降低计算成本的额外步骤,但是不是数学地要求。在一可替代实施方案中,无标注的分配可以跳过,或者等效地,预定碱基长度可以设置为很高的值,使得没有候选优化区间是无标注的。得到优化区间502处理可以输出包含优化区间514的文件、列表或其它数据结构。 文件的每个行可以包含优化区间514的开始和结束的位置。这可以包括染色体鉴定和开始 /结束位置,或者等效的信息。优化区间也可以包括计算参考分500处理期间生成的每个1 碱基变异的初始假设、或者计算局部从头区间502处理期间生成的对应于德布鲁因图路径的可替代序列。得到优化区间502处理也可以输出包含无标注区间516的文件、列表或其它数据结构,无标注区间是比预定碱基长度长和已经鉴定为处理太困难(典型地因为低读出覆盖)的优化区间。不对这些区间执行优化,并且这些区间最终在最后的变异文件32中是无标注的。每个无标注区间516可以以类似于对优化区间进行的方式鉴定。在一实施方案中, 无标注区间可以最终以单独的方式分析。在一实施方案中,得到优化区间502处理可以并行地跨计算机集群10执行,并且由各批计算的优化区间514可以组合为优化区间的单个文件和无标注区间516的单个文件。优化
优化区间514对应于参考中感兴趣的样本基因组可以具有可替代(多个)序列 (称为变异)的部分。前述步骤不仅提供区间,而且提供关于每个区间中的具体变异的试验性初始假设,这些初始假设可以是计算参考分500处理期间鉴定的碱基的最佳组合,或它们可以是对应于计算局部从头区间501的延伸期间得到的新颖的(多个)路径的(多个) 序列。使用贝叶斯公式和德布鲁因图生成假设和优化504(以下称作生成假设和优化) 处理分析每个优化区间514,并且尝试修订这些初始假设以改进它们对数据的拟合(S卩,映射配合读出观)。生成假设和优化504处理可以对每个优化区间输出高概率的优化序列假设的集合。在一实施方案中,优化504阶段使用贝叶斯公式20,其中对每个优化区间,优化区间中的初始假设的碱基是迭代地改变的,一次一个碱基,包括插入和删除的碱基,直到概率增大到最大值。贝叶斯公式20处理用于其意图目的的工作,但是在一些情况下,使用1碱基变异不能增大这些修订的序列假设的概率。在一实施方案中,在这些情况下通过使用更大的变异,例如五个碱基来增大概率。根据一实施方案,德布鲁因图22可以用于产生额外的开始序列假设作为用于优化504阶段的贝叶斯公式的变异方向。贝叶斯计分包括对每个序列假设计算概率比L(G)。在一实施方案中,优化处理计算L(G)为L(G) = log(PH/PEef)其中,Ph是有效区间中修订的序列假设的概率,且Pltef是参考的概率,以及其中在每个碱基位置在概率的计算期间使用有效区间中映射配合读出的集合。保持导致最大的概率比的修订的序列假设。该修订的序列假设被认为是对该特定的有效区间的最优假设。从当前假设得到最佳1碱基替代和选择所得更新的假设作为下一轮的基础的处理对初始序列假设的每个或子集重复。在又一实施方案中,对数似然比L(G)可以表达(以下更充分地描述)为概念上,以上略述的贝叶斯公式允许多核苷酸序列变异和重组装表示为得到给出最大可能的L(G)的基因组G的优化问题。简单的过程可以使用可以算法地描述为以下的贪心优化以设置G =(;。开始,并且计算L(G)。对G应用小变化,并且重计算L(G)。如果L(G)相对于前一值增大,保持应用的小变化。否则,取消改变。保持迭代,直到应用的小变化无一导致L(G)中的进一步增大。以该方式获得的最终G是在该样本基因组的最佳猜测。然而,不进行改进这样的过程会陷入大量的计算问题。一个问题是每次更新G时, 要求得到G上具有良好映射的配合读出200的高效方法。该映射要求配合读出标引的某种形式,这在存储器要求方面可能是有问题的。另一问题是计算要求随机访问Nd配合读出的参考分和碱基读出,这又导致了大的存储器要求。在每次迭代时重计算所有映射也是昂贵的,因此需要在小的变化应用于G时可以递增更新的数据结构。这再一次造成存储器问题,并可以引起重计算的大成本。由于这些原因,示例性实施方案提供更简单的方法,其中在优化过程期间相互地单独分析优化区间形式的样本M的小的局部区域。当优化视为有效区间的当前的优化区间时,假设正在被优化的有效区间之外任何地方G = G00使用该约束,只在有效区间优化样本多核苷酸序列G。在一实施方案中,有效区间的大小保持小于最小可能配对长度。利用该假设,对于具有对G的良好映射的配合读出200,配合读出200的一臂210必须具有近似一个配对远离有效区间的对的良好映射。该方案使得实现有效区间贪心优化过程,其中在每个迭代,仅影响有效区间的所有可能1碱基变化应用于G的当前ρ等位基因。这包括将各个碱基位置改变为不同的值, 删除其,或对其右侧或左侧插入新碱基。如果存在多于一个等位基因,依次对每个等位基因进行该优化(使剩余的等位基因不变)。对每个以该方式获得的修饰的G,计算L(G),如上所述。导致最大L(G)的1碱基变化应用于G,并且迭代直到不能获得L(G)中的进一步的改进。对关于G的所有假设,存储具有得到的最可能的假设的给定因数中的计算概率的细节, 直到最大值。模拟表明在大部分的情况中,有效区间贪心优化过程可以对L(G)的真实最优收敛,换句话说,优化过程正确地重构要组装的序列,例如基因组中局部区间的序列。然而,在很多情况中,其对L(G)的局部最优收敛。由于该原因,示例性实施方案用修改的德布鲁因图处理补充有效区间贪心优化,以为优化过程提供额外的开始序列假设(种子)。这极大地增大了优化过程的成功率。在可替代实施方案中,可以使用降低局部最优的影响的许多其它用于优化的方法,诸如模拟退火或基因算法。为了建立优化过程,鉴定如下的配合读出的集合,其具有至少一臂映射到一定距离的有效区间、臂和链附近(和可能地重叠),以便配合臂可以对有效区间中的组装有贡献。在一实施方案中,配合读出的该集合用于执行有效区间中的序列的从头组装。图9A和9B是图示用于使用贝叶斯公式和德布鲁因图生成假设和优化的处理的流程图。该处理可以包括从优化区间的集合取回优化区间,并且设置优化区间作为有效区间 (块 900)。确定靠近有效区间的映射配合读出的集合(块90 。根据一实施方案,当配合读出臂210之一的读出204可以从与臂210—致的取向的有效区间中的序列得到,配合读出视为映射到有效区间。不存在映射第二臂的要求,或者如果映射,不存在在有效区间中或靠近有效区间映射第二臂的要求。使用包含有效区间的参考的部分和关于有效区间的映射配合读出的集合,构造德布鲁因图(块904)。在可替代实施方案中,可以使用包含有效区间加某一预定长度的有效区间任一侧的序列的参考(^l的部分,构造德布鲁因图。在一实施方案中,该处理可以包括用参考顶点初始化德布鲁因图,参考顶点每个代表来自参考基因组的N个邻接碱基的碱基序列,并且其中每个碱基序列重叠紧邻的碱基序列,例如共享差不多1个碱基。通过该图的参考顶点的路径称为参考路径。映射配合读出可以用于从参考路径创建分支。取回重连接参考路径的德布鲁因图中的一个或多个序列的路径,用作用于序列假设的开始点(块906)。通过仅在有效区间中对G的当前的ρ等位基因应用所有可能的1碱基变化,优化序列假设(块90 。对每个修饰的G,使用靠近有效区间的映射配合读出,计算L(G)(块910)。导致最大L(G)的所有1碱基变化,例如大于阈值,应用于G(块912)。随后确定L(G)是否被改进(块914)。如果L(G)被改进,则该处理继续优化序列假设(块908)。如果确定L(G)没有被改进(块914),则以图9B继续,对关于G的所有优化序列假设,存储具有得到的最可能的假设的给定因数中的计算概率的细节,直到优化序列假设的最大数(块916)。在一实施方案中,保存至少给定数量的优化序列假设(或者所有优化序列假设,如果该数量没有达到预定阈值),加上概率的给定因数内的所有另外的假设,以及比该因数更差的最佳假设。随后确定德布鲁因图中是否存在更多路径(块918)。如果是,则处理通过取回另一路径继续(块906)。如果德布鲁因图中不存在更多路径,则确定是否存在更多优化区间 (块920)。如果是,处理以取回另一优化区间继续(块900)。否则,处理结束。因此,优化处理补充有其中在优化区间中计算局部从头组装的过程,以生成多个似乎可能的种子序列,而其又用于向全局最优驱动优化过程。本实施方案的局部从头组装已经实质上从在邻接的读出上使用的现有的德布鲁因图方法修改,使得适应可变有缺口读出,如下所述。用于每个局部从头组装的配合读出臂的集合选自将远离有效区间的一个配对映射到参考的配合读出。该播种过程可以导致对P (GI MtdRds) /P (G01 MtdRds)环境中局部最优的存在更有弹性的优化过程。而且,尽管一次处理一个基因组区域,但是计算远距离变异对的联合概率,导致参考中片段重复区域或其它重复要素中的伪正值的实质上的减少, 如下所述。优化处理504的输出是假设文件518,其包含对已经优化的每个选择的优化区间, 关于有效区域中G的ρ等位基因的最有可能的优化假设的列表,以及对应的计算的L (G),相对于参考多核苷酸序列假设的计算的概率比。在一实施方案中,假设文件518可以包含用于在其上执行优化的每个优化区间的块,其中该块包含1)该块引用的优化区间的开始和结束位置,其可以包括染色体标识和开始/结束位置,或者等效信息;2)用于该优化区间的优化处理期间遇到的最有可能的假设的列表。该文件可以对每个假设存储代替关于每个等位基因的参考序列的序列、以及 L(G) = log(PH/PKef)。可以并行地跨计算机集群10执行优化处理。来自计算机集群10的批处理结果可以组合为假设518的文件、列表或其它数据结构。在一实施方案中,配合读出支持文件包含要用于计算远距离变异对之间的相关性的△(配合读出)信息,以下进一步描述。提取标注提取标注506处理通过对每个优化区间分析假设文件518中存储的优化假设的列表、以及做出关于顶端假设的哪些部分有被断定的充分置信度的决定,关于在具体位置的参考多核苷酸序列推断样本多核苷酸序列的配合读出的碱基中的变异。讨论的区域是与用于优化区间的所有其它顶端假设一致的区域,其中顶端假设是具有比最可能假设的概率除以阈值高的计算概率的假设,阈值可以设置约1000的因数,即30dB。更具体地,得到在顶端假设中一致得到的变异,并且基于顶端假设对与该变异不一致的最佳假设的似然比计分每个变异。其中给出矛盾结果的最可能假设的区域是无标注的。在许多情况下,优化处理504可以导致具有单个优化假设的许多优化区间,单个优化假设计算为比所有其它假设有更多的似然,并且随后可以有置信度地标注以描述有效区间中的正确序列。然而,经常发生优化区间可以具有全部在最可能假设的小概率因数中的许多类似假设。在这样的情况下,在优化区间中G的序列上做出标注更困难。例如,考虑存储的最可能假设作为表1所示示例的示例表 权利要求
1.一种用于关于参考多核苷酸序列标注样本多核苷酸序列中的变异的计算机执行方法,该方法包括在至少一个计算机上执行定位参考多核苷酸序列中的局部区域的应用,该局部区域中存在样本多核苷酸序列的一个或多个碱基从参考多核苷酸序列中的对应碱基变化的似然, 其中至少部分基于样本多核苷酸序列的映射配对读出确定该似然;对每个局部区域生成至少一个序列假设,并且对局部区域的至少部分优化该至少一个序列假设,以对局部区域得到高概率的一个或多个优化序列假设;以及分析优化序列假设以鉴定样本多核苷酸序列中的一系列变异标注。
2.根据权利要求1的方法,其中每个配合读出包括可变有缺口读出。
3.根据权利要求1的方法,其中每个配合读出包括无缺口读出。
4.根据权利要求1的方法,还包括在数据存储库中存储参考多核苷酸序列、以及从映射到参考多核苷酸序列中的位置的样本多核苷酸序列获得的映射配合读出,并且通过在计算机集群中的多个计算机上并行执行的变异标注器来执行该方法,多个计算机的每个通过网络耦接到数据存储库。
5.根据权利要求1的方法,其中定位可能已经从参考变化的局部区域还包括使用贝叶斯公式处理计算参考分。
6.根据权利要求5的方法,其中使用贝叶斯公式处理计算参考分还包括通过以所有可能的ι碱基变异修改在Ρ等位基因中的碱基位置的碱基值,对参考多核苷酸序列中的每个碱基位置生成初始假设的集合;确定靠近参考多核苷酸序列的当前碱基位置的映射配合读出的集合;及通过对对应的集合中的每个初始假设计算概率比Pv/PMf,计算关于每个碱基位置的参考分,其中Pv是1碱基变异假设的概率,且Pm是参考多核苷酸序列中碱基值的概率,以及其中在计算在每个碱基位置的概率比期间,使用靠近每个碱基位置的映射配合读出的集合 O
7.根据权利要求6的方法,其中样本多核苷酸序列包括基因组G,以及其中每个参考分包括关于每个假设作为对数似然比的L(G),其中L(G) = log(Pv/Pref)。
8.根据权利要求6的方法,其中相互独立地生成映射配合读出,并且通过下式计算考虑所有映射配合读出的概率估计
9.根据权利要求6的方法,还包括用插入罚分近似代表(勞广,使得G的等位基因中的每个额外碱基以因数exp(-C/nD)导致P(GlMtdRds)的减小,其中nD代表每个映射配合读出中的碱基数量,使得额外的碱基不增加到G,除非额外的碱基具有充分的映射配合读出支持。
10.根据权利要求6的方法,还包括通过在计算L(G)期间仅考虑在配合读出的总和中具有小量的错配的配合读出的映射,控制单个配合读出可以向L(G)给出的贡献的量,以及用常数代表具有大量不配合的映射,该常数假设为对所有多核苷酸序列相同并且独立于G。
11.根据权利要求1的方法,其中定位可能已经从参考变化的局部区域还包括使用部分德布鲁因图计算局部从头区间以得到超过单碱基变化的变异。
12.根据权利要求11的方法,还包括用参考顶点初始化德布鲁因图,从来自参考多核苷酸序列的碱基序列创建参考顶点;对每个参考顶点,确定映射到该参考顶点、并且包括以任何可能的1碱基值延伸超过参考顶点的任一端的碱基延伸的映射配合读出的集合;至少部分基于一些具有相同延伸的映射配合读出、和具有正被处理的顶点的序列的那些映射配合读出的匹配和错配的数量,对每个碱基延伸计算代表关于以每个1碱基值延伸参考顶点的支持量的延伸强度;将与参考顶点不相容的、具有最高延伸强度的碱基延伸用作部分德布鲁因图中的分支顶点;在一个方向以深度优先方式计算关于每个分支顶点的在延伸方向的延伸强度,并且在来自具有高于阈值的延伸强度的碱基延伸的每个计算之后,创建新边沿和分支新顶点;如果路径中不存在具有高于阈值的延伸强度的碱基延伸,返回关于该路径的失败;以及如果创建等于参考顶点之一的碱基序列、并且与SNP或短插入删除一致的新分支顶点,结束计算并返回该路径。
13.根据权利要求1的方法,其中定位可能已经从参考变化的局部区域还包括通过组合由参考分和局部从头区间表示的各个可能变化的局部区域,得到优化区间。
14.根据权利要求13的方法,还包括与高概率比Pv/PMf关联的局部从头区间和参考分视为候选优化区间,其中Pv是1碱基变异假设的概率,并且PMf是参考多核苷酸序列中碱基值的概率;以及组合重叠或者相隔小于阈值碱基距离的候选优化区间为优化区间。
15.根据权利要求14的方法,其中样本多核苷酸序列包括基因组G,以及其中每个参考分包括关于每个假设作为对数似然比的L(G),其中L(G) = log(Pv/Pref)。
16.根据权利要求1的方法,其中对局部区域优化序列假设还包括遍历局部区域中初始假设中的每个碱基位置,并且用包括插入和删除碱基的每个可能的可替代碱基值迭代地变化碱基,对每个变化计算第二概率比;以及将最大化第二概率比的变化应用于局部区域。
17.根据权利要求16的方法,还包括用在正被优化的当前局部区域之外样本多核苷酸序列等于参考多核苷酸序列的假设,对相互分离的局部区域执行优化,以便仅在当前局部区域优化样本多核苷酸序列。
18.根据权利要求17的方法,还包括计算第二概率比为IV^rf,其中I3h是当前局部区域中序列假设的概率,且PMf是参考的概率。
19.根据权利要求18的方法,其中样本多核苷酸序列包括基因组G,以及其中每个参考分包括关于每个假设作为对数似然比的L(G),其中L(G) = log(PH/Pref)。
20.根据权利要求17的方法,还包括使用贝叶斯公式和德布鲁因图的组合执行优化, 其中德布鲁因图用于生成用于贝叶斯公式的另外的开始序列假设,以驱动优化处理朝向全局最优。
21.根据权利要求20的方法,还包括使用修改的德布鲁因图以处理可变有缺口读出。
22.根据权利要求17的方法,还包括要求当前局部区域的大小保持小于最小可能配对长度。
23.根据权利要求1的方法,其中分析每个序列假设的集合还包括推断在具体位置相对于的参考多核苷酸序列,样本多核苷酸序列的映射配合读出的碱基中变异。
24.根据权利要求23的方法,还包括对每个局部区域检查在对该局部区域列出的最可能假设的变异分阈值中的所有假设, 产生最可能假设的集合;得到最可能假设的集合中的每个假设中存在的公共特征;以及存储该公共特征作为各自的变异标注。
25.根据权利要求23的方法,还包括得到最可能假设的集合中的每个假设之间的不一致;以及存储该不一致作为各自的无标注区域。
26.根据权利要求1的方法,其中样本多核苷酸序列包括基因组(G),该方法还包括通过迭代地最大化该基因组的后验概率P (GlMatedReads)重组装样本多核苷酸序列,后验概率说明所有映射配合读出。
27.根据权利要求1的方法,其中分析每个序列假设的集合还包括对其中得到不一致的序列假设的集合中的假设的任一部分,生成无标注。
28.一种系统,包括数据存储库,其存储参考多核苷酸序列和从样本多核苷酸序列获得的映射配合读出, 映射配合读出映射到参考多核苷酸序列中的位置;计算机集群,包括通过网络耦接到数据存储库的多个计算机;以及在多个计算机上并行执行的变异标注器,该变异标注器配置为 基于其中一个或多个碱基可能已经从参考多核苷酸序列中的对应碱基变化的映射配合读出,定位样本多核苷酸序列中的局部区域;对每个局部区域优化序列假设,以对每个局部区域得到高概率的序列假设的集合;以及分析每个序列假设的集合以鉴定样本多核苷酸序列中的一系列变异标注。
29.—种系统,包括数据存储库,其存储参考多核苷酸序列和从样本多核苷酸序列获得的映射配合读出, 映射配合读出映射到参考多核苷酸序列中的位置;计算机集群,包括通过网络耦接到数据存储库的多个计算机;以及在多个计算机上并行执行的变异标注器,该变异标注器配置为 部分基于由基于贝叶斯公式和德布鲁因图的算法执行的证据推理的组合,对参考多核苷酸序列和映射配合读出执行统计概率分析;使用统计概率分析关于参考多核苷酸序列鉴定和标注映射配合读出中检测的变异;以及输出变异列表,每个描述以在或靠近具体位置观察到映射配合读出不同于参考多核苷酸序列的方式。
30.根据权利要求四的系统,其中变异文件还包括关于由于计算的不确定性不能标注变异的无标注区域的列表。
31.根据权利要求四的系统,其中变异包括一个或多个碱基的删除、插入、突变、多态性和重复或重排的序列。
32.根据权利要求31的系统,其中变异标注器还配置为使用概率分析从映射配合读出组装样本多核苷酸序列,其中组装的多核苷酸序列实质上基于参考多核苷酸序列,但是包括鉴定的变异。
33.根据权利要求四的系统,其中配置计算机集群,使得在多个计算机的不同计算机上执行的变异标注器的实例对参考多核苷酸序列和映射配合读出的不同部分并行地进行操作。
34.用于与参考多核苷酸序列相比标注从样本多核苷酸序列获得的映射配合读出中的变异的计算机执行方法,该方法包括部分基于由基于贝叶斯公式和德布鲁因图的算法执行的证据推理的组合,对参考多核苷酸序列和映射配合读出执行统计概率分析;使用统计概率分析关于参考多核苷酸序列鉴定和标注映射配合读出中检测的变异;以及输出变异列表,每个描述以在或靠近具体位置观察到映射配合读出不同于参考多核苷酸序列的方式,并且在存储器中存储变异标注。
35.根据权利要求34的方法,还包括在数据存储库中存储参考多核苷酸序列和从映射到参考多核苷酸序列中位置的、从样本多核苷酸序列获得的映射配合读出,并且通过在计算机集群中的多个计算机上并行执行的变异标注器来执行该方法,多个计算机的每个通过网络耦接到数据存储库。
全文摘要
提供用于与参考多核苷酸序列相比标注样本多核苷酸序列中的变异的实施方案。实施方案的各方面包括在至少一个计算机上执行定位参考多核苷酸序列中的局部区域的应用,该局部区域中存在样本多核苷酸序列的一个或多个碱基从参考多核苷酸序列中的对应碱基变化的似然,其中至少部分基于样本多核苷酸序列的映射配对读出确定该似然;对每个局部区域生成至少一个序列假设,并且对局部区域的至少部分优化该至少一个序列假设,以对局部区域得到高概率的一个或多个优化序列假设;以及分析优化序列假设以鉴定样本多核苷酸序列中的一系列变异标注。
文档编号G01N33/48GK102460155SQ201080029207
公开日2012年5月16日 申请日期2010年4月28日 优先权日2009年4月29日
发明者A.L.哈尔珀恩, B.马丁, G.尼尔森, I.纳扎伦科, J.M.巴卡施, P.卡内瓦利, R.德马纳克 申请人:考利达基因组股份有限公司