CN106021998A - 单通多变体识别计算流水线 - Google Patents

Abstract

本公开提供了一种利用计算机辅助实现的分析多个核酸片段中变体的方法。该方法使用的计算流水线包括至少一个恒定模块和至少一个可变模块，其中可变模块基于可变参数。该方法包括在处理器上执行下列步骤：接收多个核酸序列片段；为可变参数设定多个数值；将多个核酸片段通过不变模块以产生中间输出结果；将中间输出结果多次通过可变模块，每次使用可变参数的多个数值中的一个；并生成多个变体识别。

Description

单通多变体识别计算流水线

相关申请的交叉引用

本申请要求2015年3月27日提交的美国临时申请62/139,148的优先权。其全部内容在此参考并入。

发明领域

本发明主要涉及基因测序数据的分析。

背景技术

二代测序技术(NGS)为大量生产生物数据提供了强有力的工具，并为取得个性化医疗提供了帮助。虽然仅就取得序列数据来说，高通量基因测序的成本有所降低，但是分析与解读这些大规模测序数据依然存在巨大的挑战。为了识别NGS数据中的变体，大量序列比对器和变体识别器被研发出来并且被整合成各式各样的计算流水线。一个典型的计算流水线包含一个序列比对器和一个变体识别器：前者可以将序列片段与参考基因组进行比对，后者确定变异点并且向对象分配一个基因型。在计算过程中，用户为了正确分析序列数据，常常需要设定许多参数。更重要的是，一些数据需要基于细胞种类或者用于制备样品的族群来进行优化，从而准确识别变体。然而由于计算流水线每次运行都需要巨大的计算量，通过运行整个变体识别计算流水线来测试每一个参数设定在实际操作中是很难实现的。因此，有持续的需求研发新的方法和系统来优化用于分析NGS数据的参数设定。

发明概述

一方面，本发明提供了一种通过用电脑执行计算流水线来分析多个核酸序列片段中的变体的方法。该流水线包括恒定模块和可变模块，所述可变模块基于可变参数。在一些实施例中，这个方法包括用处理器执行一系列步骤：获取多个核酸序列片段；为可变参数设置多个数值；将多个核酸序列片段输送通过恒定模块以产生中间输出结果；将所述中间输出结果多次输送通过所述可变模块，每次使用所述可变参数设定的多个数值中的一个；产生多个变体识别。

在某些实施方式中，所述可变模块为变体识别模块。

在某些实施方式中，所述可变参数是先验概率。在某些实施方式中，所述先验概率是全基因组单核苷酸多态性概率，插入缺失概率，Phred-scaled间隙延伸错误概率或Phred-scaled间隙开口错误概率。

在某些实施方式中，所述恒定模块是定位模块，重复片段标记模块，局域重新比对模块，碱基质量重新校正模块或其组合。

在某些实施方式中，所述恒定模块基于恒定参数，并且所述方法包括为所述恒定参数设置数值的步骤。

在某些实施方式中，所述恒定参数选自下组：定位模块中的比对种子，定位模块中的比对种子长度，定位模块中的可允许的最大错配数，重复片段标记模块中的选择非重复片段的评分策略，重复片段标记模块中的两组重复片段的最大补偿，局域重新比对模块中的错配惩罚，局域重新比对模块中的间隙开口，局域重新比对模块中的间隙延伸，碱基质量重新校正模块中的重新校准表格，变体识别模块中的成对隐含马可夫模型，变体识别模块中的变体识别模式和变体识别模块中的变体识别阈值

在某些实施方式中，前述方法进一步包括如下步骤：基于所述多个变体识别，为所述可变参数设置矫正数值；将所述中间输出结果通过所述可变模块，其中使用所述可变参数设定的所述矫正数值；并且生成一个矫正的变体识别。

另一方面，本发明提供了一种通过用电脑执行计算流水线来分析多个核酸序列片段中的变体的方法。该流水线包括恒定模块和可变模块，所述可变模块基于可变参数。在一些实施例中，这个方法包括用处理器执行一系列步骤：获取多个核酸序列片段；为可变参数设置多个数值；将多个核酸序列片段输送通过可变模块多次以产生多个中间输出结果，每次使用所述可变参数设定的多个数值中的一个；将所述中间输出结果的每一个输送通过所述恒定模块；产生多个变体识别。

另一方面，本发明提供了一种通过用电脑执行计算流水线来分析大量核酸序列片段中的变体的方法。该流水线包括恒定模块和可变模块，所述可变模块基于可变参数。该方法包括用处理器执行一系列步骤：获取多个核酸序列片段；为可变参数设置第一数值；将多个核酸序列片段通过所述恒定模块以产生中间输出结果；将所述中间输出结果第一次通过所述可变模块，其中使用所述可变参数的第一数值；产生第一变体识别；根据第一变体识别，为可变参数设置第二数值；将所述中间输出结果第二次通过所述可变模块，使用所述可变参数的第二数值；产生第二变体识别。

另一方面，本发明提供了一种非暂时性计算机可读介质以及通过使用其存储的如前文所述的计算流水线来识别来自多个核酸序列片段的变体的方法。该计算流水线包括恒定模块和可变模块，所述可变模块基于可变参数。

在某些实施方式中，处理器执行的指令行使一系列步骤，包括：获取多个核酸序列片段；为可变参数设置多个数值；将所述多个核酸序列片段输送通过恒定模块以产生中间输出结果；将所述中间输出结果多次输送通过可变模块，其中每次使用可变参数设定的多个数值中的一个；产生变体识别。

在某些实施方式中，处理器执行的指令行使一系列步骤，包括：获取多个核酸序列片段；为可变参数设置第一数值；将多个核酸序列片段通过所述恒定模块以产生中间输出结果；将所述中间输出结果第一次通过所述可变模块，其中使用所述可变参数的第一数值；产生第一变体识别；根据第一变体识别，为可变参数设置第二数值；将所述中间输出结果第二次通过所述可变模块，使用所述可变参数的第二数值；产生第二变体识别。

该发明的上述特征优势可以通过以下的描述、权利要求和附图，得到更好的理解。

附图简要说明

图1.展示了一个示例性计算流水线。

图2.展示了一个包含了至少一个恒定模块和至少一个可变模块的示例性计算流水线。

图3.展示了一个包含了至少一个恒定模块和至少一个可变模块的示例性的计算流水线。

发明详述

在关于以上发明的概述、具体描述、如下的权利要求以及附图中，引用了本发明中的特定特征(包括步骤)。应当理解的是本发明的说明书中包含了对这些特定特征的所有可能的组合。例如，当发明的实施例或者一个特定方面或者一个权利要求中展示了一个特定的特征，该特征也可以在可能的程度上在本发明的其他方面或实施例中使用。

应当理解，除非根据上下文不允许，本说明书和权利要求中使用的单数形式“一个”包括复数形式。比如，一个“可变模块”包括一个或多个可变模块，以及本领域技术人员知道的等价形式。

一个具有两个或者更多特定步骤的方法可以被以任意顺序或者同时执行(除非在文本中排除了这种可能性)。该方法可以包含一个或多个其他步骤，这些步骤可以在任意特定步骤前面，或者在两个特定步骤之间，或者在所有的特定步骤之后被执行(除非文本中排除了这种可能性)。

在提供的值的范围内，可以理解，每个居中值，到下限的单位的十分之一，除非上下文清楚地另有规定，该范围的上限和下限和任何在该所述范围中的所述或者居中值，都被本公开内容涵盖，同时符合在所述范围内明确排除的极限。当所述范围包括一个或者两个极限，排除这两个极限或其中一个极限的范围也被包含在这个公开内容中。

为了简单清楚的阐述，当合适的时候，标号在不同的附图中重复使用，以指示相应的或类似的元件。此外，大量的具体细节被提供，以便透彻理解这里所描述的实施例的阐述。然而，本文描述的实施例可以在不存在具体细节的情况下实施。在其他实例中，方法、程序和组件没有详细描述，但没有模糊所描述的相关功能。此外，描述不应被认为是限制本文所述的实施方式的范围。应该理解的是，除非另有说明，在本公开中阐述的实施例的描述和表征并非相互排斥。

定义

本公开内容使用了如下定义：

术语“包括”以及它们在本文中的同义词代指其它组分，成分，步骤等是任选存在。例如，物品“包括”(或“其包含”)组分A，B和C可以由(即只包含)组分A，B和C，或可以不仅包含组分A，B，和C，还有一种或多种其它组分。

后跟数字的术语“至少”在本文中用于表示以该数字开始的一个范围的开始(可以是具有上限或没有上限的范围，这取决于所限定的变量)。例如，“至少1”表示1或大于1。后跟数字的术语“至多”在本文中用于表示一个以该数字结尾的范围(它可以是具有1或0作为其下限的范围，也可以是不具有下限的范围，这取决于被定义的变量)。例如，“至多4”是指4或者小于4，以及“至多40％”是指40％或低于40％。在本公开中，当一个范围被给定为“(第一数字)至(第二个数字)”或“(第一数字)-(第二数字)”，这表示范围的下限是第一数字，上限为第二数字。例如，25至100mm表示下限为25毫米，上限为100毫米的范围。

如本文所用，术语“核酸序列片段”指的是由测序方法确定的核酸序列。该核酸序列可以是DNA或RNA序列。在某些实施方案中，该核酸片段是基因组DNA测序数据。在某些实施方案中，该核酸片段是外显子测序数据。经典的DNA测序方法包括链终止法(桑格测序)。在某些实施方案中，“核酸片段”指的是由二代测序(高通量测序)的方法确定的核酸序列，其中并行测序过程中，同时产生数千或数百万的序列。二代测序方法包括，例如，通过合成技术测序(Illumina公司)，焦磷酸测序(454)，离子半导体技术(离子洪流测序)，单分子实时测序(太平洋生物科学)和通过连接测序(SOLiD测序)。根据测序方法的不同，每个核酸片段的长度可以从约30bp到大于10,000bp。例如，例如，使用SOLiD测序仪的Illumina测序产生约50bp的核酸序列片段。又例如，离子洪流测序产生高达400bp的核酸序列片段，而454焦磷酸测序产生约为700bp的核酸序列片段。再例如，单分子实时测序方法可以产生长度为10,000bp至15,000bp的核算序列片段。因此，在某些实施方式中，核酸序列片段的长度为30-100bp,50-200bp,或50-400bp。

术语“变体”当在核酸序列的语境中使用时是指一种与参照不同的核酸序列片段。典型的核酸序列变体包括但不限于单核苷酸多态性(SNP)，短缺失和插入多态性(Indel)，拷贝数变异(CNV)，微卫星标记或短串联重复序列和结构的变化。

如本文所使用的，术语“计算机实现的方法”是指相关方法是在计算机中执行，例如，一个由CPU执行的计算机程序。一种计算机，如本文所使用的，指的是可被编程以自动执行的一组算术或逻辑运算的设备(为了一般或特定目的)。计算机，如在此使用的，包括但不限于个人计算机，工作站，服务器，大型机和超级计算机。计算机可以是独立的系统，网络系统或位于计算云中的虚拟机。本公开中的方法可以通过多线程或其他并行计算方式实现。

如本文所使用的，“计算流水线”或“流水线”是指一组串联连接的数据处理元件，其中一个元件的输出是下一个元件的输入。在某些实施方案中，一个操作的输出被自动馈送到下一个操作。如本文所使用的，计算流水线中的元素可以被称为“模块”。在某些实施方案中，流水线是线性的，单向的。在某些实施方案中，主要是单向的流水线可以在其他方向上的有一些交流。在某些实施方案中，流水线可以是完全双向的。

如本文所使用的，“模块”指的是一个在计算流水线内的数据处理元件。一组模块以串联形式连接，从而形成一个计算流水线。通常，一个模块接收一个输入数据，基于所输入的数据执行特定的功能，并产生输出数据，然后将其用于随后模块的输入数据。在某些实施方案中，一个模块可以被进一步分成几个子模块，例如，子模块可以以串联形式连接。

相关术语“恒定模块”是指不具有可变参数的模块，即，当一个数据集通过一个计算流水线时，只有一组的值被设定为模块的(多个)参数。然而，应当注意的是，当向计算流水线传递不同的数据集时，不同组的值可被设定为模块的(多个)参数。

术语“可变模块”是指具有至少一个可变参数的模块。应当理解，除了包含至少一个可变参数，可变模块还可能包含至少一个恒定参数。

术语“参数”，如本文中所使用的，指的是用户需要在计算流水线或模块中设定的基准，特征或数值。当数据集通过计算流水线或模块时，基准，功能或数值被传递给函数，程序，子例程，指令，或程序。相关术语“恒定参数”指的是当运行计算流水线或模块时，该参数被设置为一个值。与此相反，“可变参数”指的是当运行计算流水线或模块时，该参数被设置为多于一个值。

术语“数值”，如本文中所使用的，是指参数设定的具体数或特征。相应的，当使用计算流水线分析核酸序列片段时，该具体的数或特征被传递给计算流水线的函数，程序，子例程，指令，或程序。为某个参数设定的具体数或特征可以根据所讨论的模块和参数确定。例如，在变体识别模块中使用的先验概率参数的数值可以是0.0005，0.0008，0.001或0.002。又例如，在比对模块中使用的错配惩罚参数的数值可以是+1到+5。

在“使数据(例如，核酸序列片段)通过流水线或一个模块”中所使用的术语“通过”指的是用流水线或者模块分析所述数据。典型的，在通过一个流水线或模块时，所述数据被作为输入数据输送到流水线或模块。该流水线或模块使用该数据运行流水线或模块中的函数，程序，子例程，指令，或程序并产生输出数据。在某些实施方式中，所产生的输出数据作为输入数据输送到另一流水线或模块中。在某些实施方案中，使数据通过一个模块还包括用所述数据作为该模块的间接输入。例如，在一个包括至少两个以串联连接的模块的计算流水线中，所述数据作为第一模块的输入被传递，其生成的输出又作为所述第二模块的输入被传递。在这样的情况下，所述数据被视为通过第二模块，尽管他们并没有被用作所述第二模块的直接输入。

术语“非暂时性计算机可读介质”指的是任何计算机可读介质，唯一的例外是一个短暂的传播信号。非暂时性计算机可读介质包括，但不限于，易失性存储器，非易失性存储器，软盘，硬盘，存储棒，寄存器存储器，处理器高速缓存和随机存取存储器。

本文所用术语“定位”或“定位到参照”是指将核酸序列片段与参照，例如，序列已知的参考基因组，比对。各种程序和算法已经被开发来将核酸序列片段定位到一个参照(参见，Flicek P，Birney E.(2009)从序列检测片段中感应：用于对其和装配的方法，Nat Methods 6(11增刊)：S6-S12；Neilsen R，Paul JS等(2011)基因型和二代测序数据中的SNP识别。Nat Rev Genet 12：443-52；Ruffalo M等(2011)二代测序哦按段对准算法的对比分析。Bioinformatics27:2790-96；Patnaik S等(2012)使用组合方法的外显子组数据分析流水线的定制PLoS ONE 7：e30080)。在各种程序和算法中，基于巴路士惠勒改造的巴路士惠勒比对法(BWA)，(Li H，Durbin R(2009年)快速，准确的巴路士惠勒短弧度排列。Bioinformatics 25：1754-60)体现了运行时间，存储器的使用和准确性之间的良好平衡，并常常被使用在不同的计算流水线中。

分析核酸序列片段的计算流水线

二代测序技术的迅速发展在过去几年变革了生物和生物医学研究。根据所使用的测序方法和系统，产生的核酸片段的数量通常在数百万以上。例如，Illumina的MiniSeq系统每次运行产生高达2500万的片段，而Illumina公司的HiSeq系列每次运行可以产生高达50亿片段。随着大量测序数据的生成，对可以用来分析和解释这些大规模测序数据的强大的计算工具的需求迫在眉睫。

已经开发了许多具有多个序列比对器和多个变体识别器的计算流水线，其中包括，但不限于，SAMtools(Li H等人(2009)序列比对/映射格式和SAMtools.Bioinformatics 25：2078-79)，glftools(Abecasis lab(2010)Abecasis lab GLF工具)，GATK(DePristo MA等(2011)用二代DNA测序数据进行基因分型和变体发现的框架，Nat Genet 43:491-98；McKenna A等。(2010)基因组分析工具包：用于分析的下一代DNA测序数据的MapReduce框架。Genome Res 20：1297-1303)和Atlas(Challis D等(2012)全外显子组的二代测序数据的综合变体分析套件，BMC Bioinformatics 13：8)。

用于分析核酸序列片段的示例性计算流水线如图1所示。图1中，用于分析核酸序列片段的计算流水线100包括串联连接的模块阵列。

首先，原始读取的数据被馈送到定位模块110来使短测序片段与参照比对对。定位模块110将核酸序列片段与参照比对，例如，一个序列已知的参考基因组。各种程序和算法已经被开发出来将核酸序列片段定位到一个参照(参见，Flicek P，Birney E.(2009)从序列片段中检测：用于对齐和装配的方法，NatMethods 6(11增刊):S6-S12；Neilsen R，Paul JS等人(2011)二代测序数据的基因型和SNP识别。Nat Rev Genet 12：443-52；Ruffalo M等(2011)二代测序片段对准算法的对比分析。Bioinformatics 27:2790-96；Patnaik S等(2012)使用组合方法的外显子组数据分析流水线的客制。PLoS ONE 7:e30080)在众多的程序和算法中，基于巴路士惠勒改造(Li H，Durbin R(2009)使用巴路士惠勒改造的快速准确短弧度排列。Bioinformatics 25:1754-60)的巴路士惠勒比对法(BWA)，展现了运行时间，存储器的使用情况和准确性之间的良好平衡，并被频繁使用在不同的计算流水线中。

比对模块的输出(例如，SAM(序列比/映射)文件或BAM(SAM的二进制版本)文件)被馈送到重复标记模块120，以去除PCT重复。在制备DNA测序样品的过程中，PCR经常被用来扩增所述片段，从而产生复制品。理想的情况下，制备的样品常常产生百分之几(例如，约4％)的彼此完全一样的片段拷贝，即，复制品。有时，30％至70％的片段是重复的。Wysoker A等(PicardTools)和Li H等(Li H等(2009)序列比对/映射格式和SAMtools,Bioinformatics 25:2078-79)。

复制被标记或者移除的片段随后被送入到局部重比对模块130，以改善片段的比对。通常情况下，相对于基准来说，局部重新比对发生在片段插入和缺失(Indel)的区域周围，并且是将片段与Indel一侧的一端相对齐，然后再对另一侧的其余部分进行定位。当片段最初定位到参照时，并不能获得任何关于插入缺失的信息。因此，定位到这样的区域的片段，只具有一小段代表插入缺失一侧的区域，并通常不会被正确定位整个indel，而是会有一端没有对齐，或者整个indel区域有很多不匹配的定位。局部重比对模块130使用其余定位到含有indel区域的片段信息，包括位于插入缺失区域更居中位置的片段，并且因此已经与所述插入缺失的两侧端部对齐。其结果是，生成了另一个和之前一样好甚至更好的定位。

局部重比对算法已由Homer等人描述。(Homer N(2010)用srma通过局域重新调二代测序数据短片段来改进变体发现。Genome Biol1 1(10):R99)。在第一步骤中，所有输入片段的对齐信息被整理在一个高效的基于图的数据结构中，其基本上类似于de-Bruijn的图表。这种调整图表示了片段是如何对齐于参考序列以及片段是如何彼此重叠。在第二步骤中，元数据是从可以反映出可以潜在地改善片段定位的对准位置的图表结构中获取，该图表结构还提供了如何重新比对片段以得到最简洁多个对准的假设。在第三步骤中，重新比对图和其元数据被用于实际执行各个片段的局域比对。

DePristo等描述了用于局域重新比对的替代算法(DePristo MA等(2011)。一种使用二代DNA测序数据的变体发现和基因型的框架。Nat Genet 43:491-98)。[君合：这个算法是和Homer的一样吗？我们使用哪一种？]该算法首先识别用于重新对准的区域，其中(i)至少一个片段包含插入缺失区域，(ii)存在错配碱基或(iii)一种已知的插入缺失的聚集(例如，从dbSNP数据库(单核苷酸多态性数据库)，这是一个公共档案，其涵盖由国家生物技术信息中心(NCBI)开发和托管的各种不同物种的普通变异，其中包含了一系列分子的变异，包括(1)单核苷酸多态性；(2)短缺失和插入多态性，(3)微卫星或短串联重复序列(STR)，(4)复核多态性(MNPs)，(5)杂序列，以及(6)命名的变体。在每个区域中，单倍型是由向参考序列定点掺入任何已知的插入缺失而生成，片段中的插入缺失来自整个位置点或者来自于所有不能完美定位参照序列的片段的史密斯-沃特曼比对(Durbin等(1998)生物序列分析：蛋白质和核酸的概率模型(剑桥大学出版社，剑桥，UK))。对于每一个单倍型H_i，片段无间距地和Hi对齐，以如下标准进行打分：

$L\left(R_j\right|H_i)=\underset{k}{\Π}L\left(R_{j,k}\right|H_{i,k})$

$L\left(R_{j,k}\right|H_{i,k})=\left\{\begin{array}{cc}1-{\mathrm{\ϵ}}_{j,k}\≈1& R_{j,k}=H_{i,k}\\ {\mathrm{\ϵ}}_{j,k}& R_{j,k}\≠H_{i,k}\end{array}\right.$

$L\left(H_i\right)=\underset{j}{\Π}L\left(R_j\right|H_i)$

其中R_j是第j个片段，k是R_j和H_i无缝对准的偏移，ε_jk是依据片段R_j的第k个碱基的声明的质量分数而决定的错误率。可以最大化L(H_i)的单倍型H_i被选为最好的可替代的单倍型。接着，所有片段根据最好的单倍型和参照(H₀)重新对准，并且每个片段R_j被分配给H_i或H₀，取决于哪一个可以最大化L(Rj|H)。如果指数几率比值或者两单倍型模型比单个基准单倍型要好至少五个指数单位，那么片段就需要被重新调整：

$\frac{L(H_0,H_i)}{L\left(H_0\right)}=\frac{{\Π}_j\max \[L\left(R_j\right|H_i),L\left(R_j\right|H_0)\]}{{\Π}_{j}L\left(R_j\right|H_0)}$

此离散化反映了精度和完整统计量的有效计算之间的折衷。在某些例子中，算法对所有个体中的所有片段同时运行，从而确保了所有个体之中的推断单倍型的一致性，这对于可靠的插入缺失识别和对比分析，例如躯体SNP和插入缺失识别，是十分关键的。通常情况下，重新调整的片段被写入到SAM/BAM中以进行进一步分析。

局域调整模块的输出随后被馈送到碱基重新校准模块140，这为每个片段中的每个碱基提供了经验性准确的碱基质量分数。在某些实施例中，碱基重新校准模块140还校正了误差协变量，例如机器周期和二核苷酸，以及测序平台特定的误差协变量例如SOLiD的颜色空间不匹配和454的流循环。碱基重新校准模块140的示例性算法已由DePristo MA等人所述。(DePristo MA等(2011)，一种使用二代DNA测序数据进行变异发现和基因型识别的框架。Nat Genet 43：491-98)。通常情况下，该算法首先列出每一个道与在所有已知数量上不改变(dbSNP build 129)的位点的参照的经验不匹配度,根据所报道的质量分数(R)、片段的机器周期(C)和二核苷酸(D)分类。每个类别的经验质量分数可以根据如下标准估计：

$mismatch(R,C,D)=\underset{r\∈R}{\Σ}\underset{c\∈C}{\Σ}\underset{d\∈D}{\Σ}\underset{br,c,d}{\Σ}br,c,d\≠bref$

$bases(R,C,D)=\underset{r\∈R}{\Σ}\underset{c\∈C}{\Σ}\underset{d\∈D}{\Σ}\left|\{b,r,c,d\}\right|$

Qempirical(R，C,D)＝(mismatch(R,C,D)+1)/(bases(R,C，D)+1)

该协变量随后被分解成线性可分的误差估计，重新校正的质量分数Q_recal可以用下式计算：

recal(r,c,d)＝Qr+ΔQ(r)+ΔΔQ(r,c)+ΔΔQ(r,d)

$\mathrm{\Δ}Q=Q_{empirical}(R,C,D)-(\underset{\mathrm{\ρ}}{\Σ}{\mathrm{\ϵ}}_r\×N_r)/bases(R,C,D)$

ΔQ(r)＝Q_empirical(r,C，D)-Q_r-ΔQ

ΔQ(r,c)＝Q_empirical(r，c，D)-(ΔQ_r+ΔQ(r))

ΔQ(r，d)＝Q_empirical(r，C，d)-(ΔQ_r+ΔQ(r))

其中，每个Δ中和ΔΔΔ是经验不匹配率和报导的所有仅包含了Q_r或者包含了协变量和Qr的观察的质量分数之间的剩余差异。其中Q_r是碱基的所报导的质量分数和ε_r为它预期的错误率；b_r,c,d是有特定协变量的碱基，且r，c，d和R，C，D分别是一系列的报导的质量得分，机器周期和二核苷酸。

碱基重新校准模块140的输出随后被送入变体识别模块150，以在包括单核苷酸多态性，短插入缺失和拷贝数变异的片段中发现具有替代等位基因的统计学证据的所有站点。通常，变体识别模块使用由涉及的一个或多个随机变量和可能的其他非随机变量的数学方程规定的算法或统计模型。例如，根据片段深度和变体计数，表明在特定位置的特定变体的置信水平为真阳性的概率可以用基于统计模型的方法以及用基准样品的局域化方法计算出来。

已经开发了各种用于变体识别的算法。例如，定位和有质量的组装(MAQ)(Li H等人(2008)使用定位质量分数以定位短DNA序列片段和识别变体。Genome Res 18：1851-58)和SOAPsnp(Li R等人(2009)大规模并行的全基因组重测序的SNP检测。Genome Res 19:1124-32)使用固定的杂合子和核苷酸-片段错误的先验概率值。SeqEM(Martin ER等人(2010)SeqEM：二代测序研究的自适应的基因型识别方法Bioinformatics 26：2803-10)引入了通过自适应方法调用期望最大化(EM)算法来估计模型参数的多样品基因型识别。SAMtools使用修订后的MAQ模型来估计测序错误。该glftools家族(glfSingle，glfMultiples和polymutt)从预先生成的基因型可能性文件(GLF)中调用SNP。GATK采用MapReduce的理念，以并行编程简单贝叶斯模型(Dean J，Ghemawat S(2008)MapReduce：大型集群的简化数据处理Commu ACM 51：107-13)。Atlas2采用已验证的全外显子组捕获测序数据，而不是常规的可能性计算的物流回归模型，并已被证明具有高灵敏度(Ji HP(2012)识别外显子变体的改进的生物信息学流水线Genome Med 4：7)。

在某些实施方案中，计算流水线可如上所述具有更少的模块。例如，Liu X等人描述的跳过局域重组模块、重复标记模块和碱基重新校正模块的流水线(Liu X等人(2013年)二代测序数据的变体识别：比较研究PLoS ONE 8(9)：e75619)。

单通多设定的计算流水线

在片段通过计算流水线时，用户需要为许多参数设值。例如，在变体识别模块中，用户需要指定先验概率。先验概率”指的是反映之前变体分布信息的概率分布，可以被用来进行变体识别。例如，在采用GATK算法的变体识别模块中，用户需要指定全基因组单核苷酸多态性(SNP)概率，和插入缺失概率。“全基因组SNP概率”是指单个核苷酸在基因组的一些特定位置发生变异的可能性，其中每个变异有一定明显的程度地存在在族群之中。“插入缺失概率”，是指在一个生物体DNA中发生插入或者缺失碱基的概率。类似的，在采用SAMTools算法的变体识别模块中，用户需要指定Phred-scaled间隙延伸错误概率和Phred-scaled间隙开口错误概率。

用户需要在流水线中设置许多其它参数。例如，一个局域对准模块中，用户需要指定针对不匹配的惩罚(一个用于对准少量零散基因序列的打分系统，例如，为了更准确地对准片段，突变被注释为序列中的间隙，该间隙经由各种惩罚计分方法惩罚，从而实现了序列对准的优化，以获得基于可用的信息的最好的对准)，间隙开口(打开任何长度的间隙所需的成本)和缺口延伸(延长现有间隙的长度所需要的成本)。碱基质量校准时，用户需要指定校准表(在表格中，为每组协变量建立跟踪匹配/不匹配的数据，并且给在一个识别组中的每个变体识别分配一个已经校准过的概率)。用户可能还需要指定变种识别模式(即在一组表示变种频率的估计值和自信指数的数据中最常出现的数值)和变体识别阈值(该程序应该发出变异位点的最小置信度阈值；如果该点相关的基因型有着比识别阈值更低的信度指数，程序会发出过滤的站点；这个阈值将高信度指数识别和低信度指数识别分离开来)，成对隐含马尔可夫模型(成对-HMM，基本HMM的一个变体对查找序列比对并评估对准符号的显著程度尤其有用)以及其中使用的转换概率等。

为了更好地评估变种，用户需要更改参数或根据一系列参数设置来研究变种。例如，用户可能希望看到当SNP/Indel的先验概率被改变时，所识别的变体的概率如何改变，所以要探索样品的各种情况。例如，如果所述样品来自癌细胞，可能需要使用较高的先验概率。根据样品的族群，需要使用不同的先验概率，因为某些族群可能有比其他族群更高的变体概率。如果用户之前不知道最佳先验概率，他们可能希望在一定范围的先验概率下探索变体特性，所以要覆盖所有可能的情况。

参数探索的一种方法是再次经历整个计算流水线来设置每个参数。例如，在用户用最可能的参数设置运行流水线后，并已经审查了被识别的变体，用户可以决定返回改变参数并重新运行整个变体识别流水线，看看这些识别的变体会有什么改变。整个管道的这样的迭代重新运行是非常耗时和极度昂贵的，从而导致只能探索有限参数的设置，因此，可能会导致某些生物学显著信息被错过。

因此，在一个方面，本公开内容提供了用计算流水线分析来自多个核酸序列片段的变体的方法和系统，该计算流水线包括至少一个恒定模块和至少一个可变模块，其中可变模块基于可变参数。在某些实施方案中，该方法包括在处理器上执行一系列步骤：接收所述多个核酸序列片段；为可变参数设定多个数值；将多个核酸序列片段通过至少一个不变模块一次；将多个核酸序列片段通过至少一个可变模块多次，每次使用可变参数设定的多个数值中的一个；并且产生多个变体识别。

图2展示了一个示例性的计算流水线，其包括至少一个恒定模块和至少一个可变模块。如图2所示，用于分析核酸序列片段的计算流水线200包括串联连接的模块阵列。

原始片段数据被馈送到定位模块210将短测序片段定位到参照。定位模块210的参数，如对准种子、种子长度和允许错配的最大数量分别被设置了一个值。映射模块210的输出被馈送到重复标记模块220。重复标记模块220的参数，例如选择非重复片段的评分策略，和两组重复片段的最大补偿分别被设置了单一数值。重复标记模块220的输出被馈送到局域重比对模块230。局域调整模块230的参数，包括不匹配惩罚，间隙开口和缺口延伸分别被设置了单一数值。局域重比对模块230的输出被馈送到碱基重新校准模块240。碱基重新校准模块240的参数包括重新校准表(上文详细描述过)中的值被分别设置为单一数值。碱基重新校准模块240的输出被馈送到变体识别模块250。多组数值251，252，253被赋予给各个参数，例如，先验概率。因而，输入的数据通过变体识别模块250多次，通过将每组数值(251，252和253)赋予给各个参数。在一个实例中，对一个3×3矩阵：SNP(0.002，0.001，0.0005)×Indel(0.0002，0.0001,0.00005)进行了探讨，其中，SNP 0.002指SNP先验概率被设定为0.002，如此等等。因此，该变体识别模块用3x3不同的参数组合运行9次，生产9个变体识别输出。

在某些实施方案中，一些变体识别模块250中的部分计算，如片段过滤，仅取决于片段碱基的品质和对准，以及碱的某些条件概率，不依赖于现有的概率的，所以只需被运行一次。

通过去除冗余运算，变体识别模块250中使用多个参数设置的多个输出可以靠通过计算流水线200一次来获得。一次通过需要的总计算量仅仅是多次通过完整流水线的计算量的一小部分。

在这个单次通过多参数设定的运行之后，变体识别模块250的输出可以在视觉上显示为参数集合的一个函数，使用户对变体有更好的视觉上的认识。

逐步勘探计算流水线

相对于计算并且输出多个参数在一次通过的运行结果，该流水线可以计算并输出在一个参数设置下的第一次通过结果。第一次通过时，每个模块的输出，例如，输出中间结果，被保存下来。用户可以根据该结果，为至少一个流水线的模块，例如，第一次通过的变体识别，调整参数设定，并使片段第二次通过流水线。在第二次通过中，流水线可以自动识别不需要重新计算的中间结果，并且直接在新参数设定中使用它们。这样，对于附加的参数设置仅需要计算那些被新旧参数设置差异影响的流水线的部分，因此，节省了参数勘探中的计算量。该方法的优势也在于它可以避免对一些不必要的参数设定所进行的探索。

因此，在另一个方面，本公开内容提供了使用包括至少一个恒定模块和至少一个可变模块的计算流水线，来分析多个核算片段中的变体的方法与系统，其中可变模块基于可变参数。在某些实施方案中，该方法包括在处理器上执行一系列步骤：接收多个核酸序列片段，其中至少一个核酸片段包括一个变体；为可变参数设置第一数值；让多个核酸片段通过至少一个恒定模块一次；通过至少一个可变模块，其中使用可变参数的第一数值；产生第一变体识别；根据第一变体识别为可变参数设置第二数值；让多个核酸片段第二次通过至少一个可变模块，其中使用可变参数的第二数值；并且在不通过至少一个恒定模块的情况下，产生第二个变体识别。

图3展示了用于分析核酸序列片段的示例性增量勘探计算流水线300。计算流水线300包括串联连接的模块阵列。如上文所详细描述地，当原始片段数据被馈送到计算流水线300时，它会通过定位模块310，重复标记模块320，局域重比对模块330，碱基重新校准模块340和变体识别模块350。第一个数值集合被用于变体识别模块350中的先验概率，并且获得第一个变体识别的结果。在第一个变体识别的结果的基础上，第二组数值被设置为变体识别模块350中的概率。碱基重新校准340的输出被馈送到变体识别模块350，以产生一个第二个变体识别结果。

在某些实施方案中，“单通多设置”模式和“增量勘探”模式可以结合起来。例如，流水线在第一个单次通过中运行参数设置，产生多个变体识别，同时还保存了中间值的输出。在审查完在第一单次通过中产生的多个变体识别后，用户可以调整一个或多个模块的参数。然后流水线可以将第一通中的中间结果，重新用在新的通过中，从而节省了在新的通过中需要的计算。

定位模块

在另一个方面，本公开提供了改进对比参照计算的方法和系统。在定位过程中，用户也可以选择一些参数，例如，对准种子，种子长度，允许的最大错配数量，等等。当用户探索对准参数时，已经与基准完美匹配的片段不需要在每次参数设定中重复被对准，而是可以只运行一次，其结果可以被用在所有情况下。只有没有定位的片段以及与参照有错配的片段才需要在不同的参数制定下，反复多次对准。

组装流水线

在另一方面，本公开提供了从头组装核酸片段的方法和系统，例如不需借助参照(即，参照基因组序列)来进行组装。从头组件通常被用于研究非模式生物体，因为他们的参照基因组序列常常不存在。从头组件还被用于转录序列数据。一些从头组装的程序已经被开发出来，包括SOAPdenovo-Trans，Velvet/Oases，Trans-ABySS，和Trinity。两种一般被用于从头装配的基本算法是：重叠图形和德布鲁因图。虽然重叠图大多用于桑格测序片段，因为它比德布鲁因图具有更大的计算强度，并且以最高的效率装配更少的高度重叠片段。德布鲁因图基于k-1的序列守恒对准K-mers(通常为25-50bp)以创建重叠群(连续，重复序列片段)。对于比片段长度更短的K-mers的使用，使得德布鲁因图表法降低了计算的强度，并且适用于二代测序数据。

在某些实施方案中，从头装配的计算流水线包括至少一个恒定模块，例如，片段清洁和过滤模块，和至少一个可变模块，例如，重叠群模块。在通过组装流水线时，多个值被赋予给可变模块的可变参数，例如，重叠群模块的K-mers，同时单值被赋给不变模块的恒定参数。因此，不同参数设置之间的共同计算可以被节省下来。

本发明可以用如下实施例进行说明。这些实施例不旨在从任何方面限制本发明。

实施例1

以下是一个用于分析核酸序列片段的单通多设置的计算流水线的例子。

该计算流水线由以下模块串联而成：定位模块，重复标记模块，局域重新比对模块，碱基质量重新校准模块和变体识别模块。

计算流水线的实施使用SuperServer 6016T-NTF服务器(Intel Dual XeonE5620 2.4 Hz CPU，48 GB DDR3 1333MHz Ecc Memory,和2 TB SATA SeagateST32000644NS HDD)。

用于分析的核酸序列片段来自人类基因组外显子组测序数据。

具有FastQ格式的核酸片段(由Illumina公司系统生产)被馈送到定位模块，该定位模块用BWA算法将片段对应到参照基因组(human glk v37)。SAM格式的定位文件被转换为BAM格式并且由GATK进行排序。排序后的BAM文件被馈送到重复标记模块，采用Picard算法删除重复的PCR。去掉重复部分的BAM文件被输送到使用GATK算法的局域重新比对模块以及作为基准的变量识别格式(vcf)的dbSNP。重新对准的BAM文件被馈送到使用GATK算法的碱基质量重新校准模块。

上述模块中所用参数的具体数值设置如下：重复标记模块中的两组重复片段的最大补偿为100，局部重比对模块的LOD阈值为5.0。

重新校准的BAM文件随后被馈送到使用GATK算法的变体识别模块。这个模块被运行了9次以探索3×3矩阵：SNP(0.002,0.001,0.0005)×插入缺失(Indel)(0.0002,0.0001,0.00005)。

结果：大约需要2300分钟来运行上述单通多设置计算流水线(定位模块大约需要40分钟，重复片段标记模块大约需要40分钟，局部重比对模块大约需要40分钟，碱基质量重新校准模块需要大约180分钟，变体识别模块大约需要2000分钟)。作为比较，如果使用单设置的计算流水线运行9次以探索上述3×3矩阵则需要大约4800分钟。通过使用变体识别模块的矩阵中重新校准的BAM文件，上述流水线中所使用的总计算量与通过全流水线9次以探讨矩阵所需计算量相比，仅仅是很小的一部分。

虽然本发明已经被具体地展示，并依据具体实施方案给予了详细说明(其中部分是优选的实施方案)，但需要被本领域技术人员明白的是，在不脱离本文公开的发明的精神和范围的情况下，依然可以对其做出各种形式和细节上的调整。

Claims (20)

1.一种使用计算流水线分析多个核酸序列片段中的变体的计算机辅助实现方法，其中所述计算流水线包括恒定模块和可变模块，所述可变模块基于可变参数，其特征在于，所述方法包括在处理器上执行如下步骤：

接收多个核酸序列片段；

为所述可变参数设定多个数值；

将所述多个核酸序列片段通过所述不变模块以产生中间输出结果；

将所述中间输出结果通过所述可变模块多次，其中每次使用所述可变参数设定的所述多个数值中的一个；并且

生成多个变体识别。

2.如权利要求1所述的方法，其特征在于，所述可变模块为变体识别模块。

3.如权利要求2所述的方法，其特征在于，所述可变参数是先验概率。

4.如权利要求4所述的方法，其特征在于，所述先验概率是全基因组单核苷酸多态性概率，插入缺失概率，Phred-scaled间隙延伸错误概率或Phred-scaled间隙开口错误概率。

5.如权利要求1所述的方法，其特征在于，所述恒定模块是定位模块，重复片段标记模块，局域重新比对模块，碱基质量重新校正模块或其组合。

6.如权利要求1所述的方法，其特征在于，所述恒定模块基于恒定参数，并且所述方法包括为所述恒定参数设置数值的步骤。

7.如权利要求1所述的方法，进一步包括如下步骤：

基于所述多个变体识别，为所述可变参数设置矫正数值；

将所述中间输出结果通过所述可变模块，其中使用所述矫正数值；并且

生成一个矫正的变体识别。

8.一种非暂态计算机可读介质，其附带指令可以使其利用其所存诸的计算流水线分析来自多个核酸序列片段的变体，所述计算流水线包括恒定模块和可变模块，所述至少一个可变模块基于可变参数，所述指令当被处理器执行时，会进行如下操作：

接收多个核酸序列片段；

为所述可变参数设定多个数值；

将所述多个核酸序列片段通过所述恒定模块以产生中间输出结果；

将所述中间输出结果通过所述可变模块多次，其中每次使用所述可变参数设定的所述多个数值中的一个；并且

生成多个变体识别。

9.如权利要求8所述的非暂态计算机可读介质，其特征在于，所述可变模块为变体识别模块。

10.如权利要求9所述的非暂态计算机可读介质，其特征在于，所述可变参数是先验概率。

11.如权利要求10所述的非暂态计算机可读介质，其特征在于，所述先验概率是全基因组单核苷酸多态性概率，插入缺失概率，Phred-scaled间隙延伸错误概率或Phred-scaled间隙开口错误概率。

12.如权利要求8所述的非暂态计算机可读介质，其特征在于，所述恒定模块是定位模块，重复片段标记模块，局域重新比对模块，碱基质量重新校正模块或其组合。

13.如权利要求8所述的非暂态计算机可读介质，其特征在于，所述恒定模块基于恒定参数，并且所述方法包括为所述恒定参数设置数值的步骤。

14.如权利要求8所述的非暂态计算机可读介质，其特征在于，所述指令当被处理器执行时，会进一步进行如下操作：

基于所述多个变体识别，为所述可变参数设置一个矫正数值；

将所述中间输出结果通过所述可变模块，其中使用所述可变参数设定的所述矫正数值；并且

生成一个矫正的变体识别。

15.一种使用计算流水线分析多个核酸序列片段中的变体的计算机辅助实现方法，其中所述计算流水线包括恒定模块和可变模块，所述可变模块基于可变参数，其特征在于，所述方法包括在处理器上执行如下步骤：

接收多个核酸序列片段；

为所述可变参数设定第一数值；

将所述多个核酸序列片段通过所述不变模块以产生中间输出结果；

将所述中间输出结果第一次通过所述可变模块，并使用所述可变参数的所述第一数值；

生成第一变体识别；

根据所述第一变体识别，为所述可变参数设置第二数值；

将所述中间输出结果第二次通过所述可变模块，并使用所述可变参数的所述第二数值；并且

生成第二变体识别。

16.如权利要求15所述的方法，其特征在于，所述可变模块为变体识别模块。

17.如权利要求16所述的方法，其特征在于，所述可变参数是先验概率。

18.如权利要求17所述的方法，其特征在于，所述先验概率是全基因组单核苷酸多态性概率，插入缺失概率，Phred-scaled间隙延伸错误概率或Phred-scaled间隙开口错误概率。

19.如权利要求15所述的方法，其特征在于，所述恒定模块是定位模块，重复片段标记模块，局域重新比对模块，碱基质量重新校正模块或其组合。

20.如权利要求15所述的方法，其特征在于，所述恒定模块基于恒定参数，并且所述方法包括为所述恒定参数设置数值的步骤。