单通多变体识别计算流水线的制作方法

单通多变体识别计算流水线的制作方法
【专利摘要】本公开提供了一种利用计算机辅助实现的分析多个核酸片段中变体的方法。该方法使用的计算流水线包括至少一个恒定模块和至少一个可变模块，其中可变模块基于可变参数。该方法包括在处理器上执行下列步骤：接收多个核酸序列片段；为可变参数设定多个数值；将多个核酸片段通过不变模块以产生中间输出结果；将中间输出结果多次通过可变模块，每次使用可变参数的多个数值中的一个；并生成多个变体识别。
【专利说明】单通多变体识别计算流水线
[0001] 相关申请的交叉引用
[0002] 本申请要求2015年3月27日提交的美国临时申请62/139,148的优先权。其全部内 容在此参考并入。 发明领域
[0003] 本发明主要设及基因测序数据的分析。
【背景技术】
[0004] 二代测序技术(NGS)为大量生产生物数据提供了强有力的工具，并为取得个性化 医疗提供了帮助。虽然仅就取得序列数据来说，高通量基因测序的成本有所降低，但是分析 与解读运些大规模测序数据依然存在巨大的挑战。为了识别NGS数据中的变体，大量序列比 对器和变体识别器被研发出来并且被整合成各式各样的计算流水线。一个典型的计算流水 线包含一个序列比对器和一个变体识别器:前者可W将序列片段与参考基因组进行比对， 后者确定变异点并且向对象分配一个基因型。在计算过程中，用户为了正确分析序列数据， 常常需要设定许多参数。更重要的是，一些数据需要基于细胞种类或者用于制备样品的族 群来进行优化，从而准确识别变体。然而由于计算流水线每次运行都需要巨大的计算量，通 过运行整个变体识别计算流水线来测试每一个参数设定在实际操作中是很难实现的。因 此，有持续的需求研发新的方法和系统来优化用于分析NGS数据的参数设定。
[00化]发明概述
[0006] 一方面，本发明提供了一种通过用电脑执行计算流水线来分析多个核酸序列片段 中的变体的方法。该流水线包括恒定模块和可变模块，所述可变模炔基于可变参数。在一些 实施例中，运个方法包括用处理器执行一系列步骤:获取多个核酸序列片段;为可变参数设 置多个数值;将多个核酸序列片段输送通过恒定模块W产生中间输出结果;将所述中间输 出结果多次输送通过所述可变模块，每次使用所述可变参数设定的多个数值中的一个;产 生多个变体识别。
[0007] 在某些实施方式中，所述可变模块为变体识别模块。
[000引在某些实施方式中，所述可变参数是先验概率。在某些实施方式中，所述先验概率 是全基因组单核巧酸多态性概率，插入缺失概率，Phred-scaled间隙延伸错误概率或 I%red-scaled间隙开口错误概率。
[0009] 在某些实施方式中，所述恒定模块是定位模块，重复片段标记模块，局域重新比对 模块，碱基质量重新校正模块或其组合。
[0010] 在某些实施方式中，所述恒定模炔基于恒定参数，并且所述方法包括为所述恒定 参数设置数值的步骤。
[0011] 在某些实施方式中，所述恒定参数选自下组:定位模块中的比对种子，定位模块中 的比对种子长度，定位模块中的可允许的最大错配数，重复片段标记模块中的选择非重复 片段的评分策略，重复片段标记模块中的两组重复片段的最大补偿，局域重新比对模块中 的错配惩罚，局域重新比对模块中的间隙开口，局域重新比对模块中的间隙延伸，碱基质量 重新校正模块中的重新校准表格，变体识别模块中的成对隐含马可夫模型，变体识别模块 中的变体识别模式和变体识别模块中的变体识别阔值
[0012] 在某些实施方式中，前述方法进一步包括如下步骤:基于所述多个变体识别，为所 述可变参数设置矫正数值;将所述中间输出结果通过所述可变模块，其中使用所述可变参 数设定的所述矫正数值;并且生成一个矫正的变体识别。
[0013] 另一方面，本发明提供了一种通过用电脑执行计算流水线来分析多个核酸序列片 段中的变体的方法。该流水线包括恒定模块和可变模块，所述可变模炔基于可变参数。在一 些实施例中，运个方法包括用处理器执行一系列步骤:获取多个核酸序列片段;为可变参数 设置多个数值;将多个核酸序列片段输送通过可变模块多次W产生多个中间输出结果，每 次使用所述可变参数设定的多个数值中的一个;将所述中间输出结果的每一个输送通过所 述恒定模块;产生多个变体识别。
[0014] 另一方面，本发明提供了一种通过用电脑执行计算流水线来分析大量核酸序列片 段中的变体的方法。该流水线包括恒定模块和可变模块，所述可变模炔基于可变参数。该方 法包括用处理器执行一系列步骤:获取多个核酸序列片段;为可变参数设置第一数值;将多 个核酸序列片段通过所述恒定模块W产生中间输出结果;将所述中间输出结果第一次通过 所述可变模块，其中使用所述可变参数的第一数值;产生第一变体识别；根据第一变体识 另IJ，为可变参数设置第二数值;将所述中间输出结果第二次通过所述可变模块，使用所述可 变参数的第二数值;产生第二变体识别。
[0015] 另一方面，本发明提供了一种非暂时性计算机可读介质W及通过使用其存储的如 前文所述的计算流水线来识别来自多个核酸序列片段的变体的方法。该计算流水线包括恒 定模块和可变模块，所述可变模炔基于可变参数。
[0016] 在某些实施方式中，处理器执行的指令行使一系列步骤，包括:获取多个核酸序列 片段;为可变参数设置多个数值;将所述多个核酸序列片段输送通过恒定模块W产生中间 输出结果;将所述中间输出结果多次输送通过可变模块，其中每次使用可变参数设定的多 个数值中的一个;产生变体识别。
[0017] 在某些实施方式中，处理器执行的指令行使一系列步骤，包括:获取多个核酸序列 片段;为可变参数设置第一数值;将多个核酸序列片段通过所述恒定模块W产生中间输出 结果;将所述中间输出结果第一次通过所述可变模块，其中使用所述可变参数的第一数值； 产生第一变体识别;根据第一变体识别，为可变参数设置第二数值;将所述中间输出结果第 二次通过所述可变模块，使用所述可变参数的第二数值;产生第二变体识别。
[0018] 该发明的上述特征优势可W通过W下的描述、权利要求和附图，得到更好的理解。
[0019] 附图简要说明
[0020] 图1.展示了一个示例性计算流水线。
[0021] 图2.展示了一个包含了至少一个恒定模块和至少一个可变模块的示例性计算流 水线。
[0022] 图3.展示了一个包含了至少一个恒定模块和至少一个可变模块的示例性的计算 流水线。
[00剖发明详述
[0024] 在关于W上发明的概述、具体描述、如下的权利要求W及附图中，引用了本发明中 的特定特征(包括步骤）。应当理解的是本发明的说明书中包含了对运些特定特征的所有可 能的组合。例如，当发明的实施例或者一个特定方面或者一个权利要求中展示了一个特定 的特征，该特征也可W在可能的程度上在本发明的其他方面或实施例中使用。
[0025] 应当理解，除非根据上下文不允许，本说明书和权利要求中使用的单数形式"一 个"包括复数形式。比如，一个"可变模块"包括一个或多个可变模块，W及本领域技术人员 知道的等价形式。
[0026] -个具有两个或者更多特定步骤的方法可W被W任意顺序或者同时执行(除非在 文本中排除了运种可能性）。该方法可W包含一个或多个其他步骤，运些步骤可W在任意特 定步骤前面，或者在两个特定步骤之间，或者在所有的特定步骤之后被执行(除非文本中排 除了运种可能性）。
[0027] 在提供的值的范围内，可W理解，每个居中值，到下限的单位的十分之一，除非上 下文清楚地另有规定，该范围的上限和下限和任何在该所述范围中的所述或者居中值，都 被本公开内容涵盖，同时符合在所述范围内明确排除的极限。当所述范围包括一个或者两 个极限，排除运两个极限或其中一个极限的范围也被包含在运个公开内容中。
[0028] 为了简单清楚的阐述，当合适的时候，标号在不同的附图中重复使用，W指示相应 的或类似的元件。此外，大量的具体细节被提供，W便透彻理解运里所描述的实施例的阐 述。然而，本文描述的实施例可W在不存在具体细节的情况下实施。在其他实例中，方法、程 序和组件没有详细描述，但没有模糊所描述的相关功能。此外，描述不应被认为是限制本文 所述的实施方式的范围。应该理解的是，除非另有说明，在本公开中阐述的实施例的描述和 表征并非相互排斥。
[0029]
[0030] 本公开内容使用了如下定义：
[0031] 术语"包括"W及它们在本文中的同义词代指其它组分，成分，步骤等是任选存在。 例如，物品"包括"（或"其包含"）组分A，B和C可W由（即只包含)组分A，B和C，或可W不仅包 含组分A，B，和C，还有一种或多种其它组分。
[0032] 后跟数字的术语"至少"在本文中用于表示W该数字开始的一个范围的开始(可W 是具有上限或没有上限的范围，运取决于所限定的变量）。例如，"至少1"表示1或大于1。后 跟数字的术语"至多"在本文中用于表示一个W该数字结尾的范围（它可W是具有1或0作为 其下限的范围，也可W是不具有下限的范围，运取决于被定义的变量）。例如，"至多4"是指4 或者小于4, W及"至多40%"是指40%或低于40%。在本公开中，当一个范围被给定为"(第 一数字)至(第二个数字r或"(第一数字)-(第二数字)"，运表示范围的下限是第一数字，上 限为第二数字。例如，25至IOOmm表示下限为25毫米，上限为100毫米的范围。
[0033] 如本文所用，术语"核酸序列片段"指的是由测序方法确定的核酸序列。该核酸序 列可W是DNA或RNA序列。在某些实施方案中，该核酸片段是基因组DNA测序数据。在某些实 施方案中，该核酸片段是外显子测序数据。经典的DNA测序方法包括链终止法(桑格测序）。 在某些实施方案中，"核酸片段"指的是由二代测序(高通量测序）的方法确定的核酸序列， 其中并行测序过程中，同时产生数千或数百万的序列。二代测序方法包括，例如，通过合成 技术测序（Illumina公司），焦憐酸测序(454)，离子半导体技术(离子洪流测序），单分子实 时测序(太平洋生物科学)和通过连接测序(s化iD测序）。根据测序方法的不同，每个核酸片 段的长度可W从约30bp到大于10，000bp。例如，例如，使用SOLiD测序仪的Illumina测序产 生约50bp的核酸序列片段。又例如，离子洪流测序产生高达4(K)bp的核酸序列片段，而454焦 憐酸测序产生约为700bp的核酸序列片段。再例如，单分子实时测序方法可W产生长度为 10 ,OOObp至15 ,OOObp的核算序列片段。因此，在某些实施方式中，核酸序列片段的长度为 30-100bp，50-200bp，或50-400bp。
[0034] 术语"变体"当在核酸序列的语境中使用时是指一种与参照不同的核酸序列片段。 典型的核酸序列变体包括但不限于单核巧酸多态性(SNP)，短缺失和插入多态性（Indel), 拷贝数变异(CNV)，微卫星标记或短串联重复序列和结构的变化。
[0035] 如本文所使用的，术语"计算机实现的方法"是指相关方法是在计算机中执行，例 如，一个由CPU执行的计算机程序。一种计算机，如本文所使用的，指的是可被编程W自动执 行的一组算术或逻辑运算的设备(为了一般或特定目的）。计算机，如在此使用的，包括但不 限于个人计算机，工作站，服务器，大型机和超级计算机。计算机可W是独立的系统，网络系 统或位于计算云中的虚拟机。本公开中的方法可W通过多线程或其他并行计算方式实现。
[0036] 如本文所使用的，"计算流水线"或"流水线"是指一组串联连接的数据处理元件， 其中一个元件的输出是下一个元件的输入。在某些实施方案中，一个操作的输出被自动馈 送到下一个操作。如本文所使用的，计算流水线中的元素可W被称为"模块"。在某些实施方 案中，流水线是线性的，单向的。在某些实施方案中，主要是单向的流水线可W在其他方向 上的有一些交流。在某些实施方案中，流水线可W是完全双向的。
[0037] 如本文所使用的，"模块"指的是一个在计算流水线内的数据处理元件。一组模块 W串联形式连接，从而形成一个计算流水线。通常，一个模块接收一个输入数据，基于所输 入的数据执行特定的功能，并产生输出数据，然后将其用于随后模块的输入数据。在某些实 施方案中，一个模块可W被进一步分成几个子模块，例如，子模块可W W串联形式连接。
[0038] 相关术语"恒定模块"是指不具有可变参数的模块，即，当一个数据集通过一个计 算流水线时，只有一组的值被设定为模块的（多个)参数。然而，应当注意的是，当向计算流 水线传递不同的数据集时，不同组的值可被设定为模块的（多个)参数。
[0039] 术语"可变模块"是指具有至少一个可变参数的模块。应当理解，除了包含至少一 个可变参数，可变模块还可能包含至少一个恒定参数。
[0040] 术语"参数"，如本文中所使用的，指的是用户需要在计算流水线或模块中设定的 基准，特征或数值。当数据集通过计算流水线或模块时，基准，功能或数值被传递给函数，程 序，子例程，指令，或程序。相关术语"恒定参数"指的是当运行计算流水线或模块时，该参数 被设置为一个值。与此相反，"可变参数"指的是当运行计算流水线或模块时，该参数被设置 为多于一个值。
[0041] 术语"数值"，如本文中所使用的，是指参数设定的具体数或特征。相应的，当使用 计算流水线分析核酸序列片段时，该具体的数或特征被传递给计算流水线的函数，程序，子 例程，指令，或程序。为某个参数设定的具体数或特征可W根据所讨论的模块和参数确定。 例如，在变体识别模块中使用的先验概率参数的数值可W是0.0005,0.0008,0.OOl或 0.002。又例如，在比对模块中使用的错配惩罚参数的数值可W是+1到巧。
[0042] 在"使数据(例如，核酸序列片段)通过流水线或一个模块"中所使用的术语"通过" 指的是用流水线或者模块分析所述数据。典型的，在通过一个流水线或模块时，所述数据被 作为输入数据输送到流水线或模块。该流水线或模块使用该数据运行流水线或模块中的函 数，程序，子例程，指令，或程序并产生输出数据。在某些实施方式中，所产生的输出数据作 为输入数据输送到另一流水线或模块中。在某些实施方案中，使数据通过一个模块还包括 用所述数据作为该模块的间接输入。例如，在一个包括至少两个W串联连接的模块的计算 流水线中，所述数据作为第一模块的输入被传递，其生成的输出又作为所述第二模块的输 入被传递。在运样的情况下，所述数据被视为通过第二模块，尽管他们并没有被用作所述第 二模块的直接输入。
[0043] 术语"非暂时性计算机可读介质"指的是任何计算机可读介质，唯一的例外是一个 短暂的传播信号。非暂时性计算机可读介质包括，但不限于，易失性存储器，非易失性存储 器，软盘，硬盘，存储棒，寄存器存储器，处理器高速缓存和随机存取存储器。
[0044] 本文所用术语"定位"或"定位到参照"是指将核酸序列片段与参照，例如，序列已 知的参考基因组，比对。各种程序和算法已经被开发来将核酸序列片段定位到一个参照(参 见，Flicek P，Birney E. (2009)从序列检测片段中感应：用于对其和装配的方法，Nat Methods 6(11 增刊）：S6-S12;Neilsen R，Paul JS等（2011)基因型和二代测序数据中的SNP 识别。化t Rev Genet 12:443-52;Ruffalo M等（2011)二代测序哦按段对准算法的对比分 析。Bioinformatics 27: 2790-96;化tnaik S等（2012)使用组合方法的外显子组数据分析 流水线的定制化OS 0肥7:630080)。在各种程序和算法中，基于己路±惠勒改造的己路± 惠勒比对法（BWA)，化i H，Durbin R(2009年）快速，准确的己路±惠勒短弧度排列。 Bioinformatics 25:1754-60)体现了运行时间，存储器的使用和准确性之间的良好平衡， 并常常被使用在不同的计算流水线中。
[0045] 分析核酸序列片段的计算流水线
[0046] 二代测序技术的迅速发展在过去几年变革了生物和生物医学研究。根据所使用的 测序方法和系统，产生的核酸片段的数量通常在数百万W上。例如，Illumina的MiniSeq系 统每次运行产生高达2500万的片段，而Illumina公司的化Seq系列每次运行可W产生高达 50亿片段。随着大量测序数据的生成，对可W用来分析和解释运些大规模测序数据的强大 的计算工具的需求迫在眉睫。
[0047] 已经开发了许多具有多个序列比对器和多个变体识别器的计算流水线，其中包 括，但不限于，SMtoolS(Li H等人(2009)序列比对/映射格式和SAMtoolS.Bioinformatics 25:2078-79)，glftools(Abecasis lab(2010)Abecasis Iab GLF工具），GATK化ePristo MA 等（2011)用二代DNA测序数据进行基因分型和变体发现的框架，Nat Genet 43:491-98; McKenna A等。（2010)基因组分析工具包：用于分析的下一代DNA测序数据的MapReduce框 架。Genome Res 20:1297-1303)和Atlas(化allis D等(2012)全外显子组的二代测序数据 的综合变体分析套件，BMC Bioinformatics 13:8)。
[0048] 用于分析核酸序列片段的示例性计算流水线如图1所示。图1中，用于分析核酸序 列片段的计算流水线100包括串联连接的模块阵列。
[0049] 首先，原始读取的数据被馈送到定位模块110来使短测序片段与参照比对对。定位 模块110将核酸序列片段与参照比对，例如，一个序列已知的参考基因组。各种程序和算法 已经被开发出来将核酸序列片段定位到一个参照（参见，Flicek P，Birney E. (2009)从序 列片段中检测：用于对齐和装配的方法，化t Methods 6(11增刊）：S6-S12;Neilsen R，Paul JS等人(2011)二代测序数据的基因型和SNP识别。化t Rev Genet 12:443-52;Ruffalo M等 (2011) 二代测序片段对准算法的对比分析。Bioinformatics 27 : 2790-96 ;化tnaik S等 (2012) 使用组合方法的外显子组数据分析流水线的客制。PLoS 0肥7:630080)在众多的程 序和算法中，基于己路±惠勒改造化i H，Durbin R(2009)使用己路±惠勒改造的快速准确 短弧度排列。Bioinformatics 25:1754-60)的己路±惠勒比对法(BWA)，展现了运行时间， 存储器的使用情况和准确性之间的良好平衡，并被频繁使用在不同的计算流水线中。
[0050] 比对模块的输出（例如，SAM(序列比/映射)文件或BAM(SAM的二进制版本)文件)被 馈送到重复标记模块120, W去除PCT重复。在制备DNA测序样品的过程中，PCR经常被用来扩 增所述片段，从而产生复制品。理想的情况下，制备的样品常常产生百分之几(例如，约4%) 的彼此完全一样的片段拷贝，即，复制品。有时，30%至70%的片段是重复的。Wysoker A等 (PicardTools)和Li H等化i H等（2009)序列比对/映射格式和SAMtools,Bioinformatics 25:2078-79)。
[0051] 复制被标记或者移除的片段随后被送入到局部重比对模块130, W改善片段的比 对。通常情况下，相对于基准来说，局部重新比对发生在片段插入和缺失（Indel)的区域周 围，并且是将片段与Indel-侧的一端相对齐，然后再对另一侧的其余部分进行定位。当片 段最初定位到参照时，并不能获得任何关于插入缺失的信息。因此，定位到运样的区域的片 段，只具有一小段代表插入缺失一侧的区域，并通常不会被正确定位整个indel，而是会有 一端没有对齐，或者整个indel区域有很多不匹配的定位。局部重比对模块130使用其余定 位到含有indel区域的片段信息，包括位于插入缺失区域更居中位置的片段，并且因此已经 与所述插入缺失的两侧端部对齐。其结果是，生成了另一个和之前一样好甚至更好的定位。 [0化2] 局部重比对算法已由化mer等人描述。巧omer N(2010)用srma通过局域重新调二 代测序数据短片段来改进变体发现。Genome Bioll 1(10):R99)。在第一步骤中，所有输入 片段的对齐信息被整理在一个高效的基于图的数据结构中，其基本上类似于de-BruUn的 图表。运种调整图表示了片段是如何对齐于参考序列W及片段是如何彼此重叠。在第二步 骤中，元数据是从可W反映出可W潜在地改善片段定位的对准位置的图表结构中获取，该 图表结构还提供了如何重新比对片段W得到最简洁多个对准的假设。在第=步骤中，重新 比对图和其元数据被用于实际执行各个片段的局域比对。
[0053] DePristo等描述了用于局域重新比对的替代算法(DePristo MA等（2011)。一种使 用二代DNA测序数据的变体发现和基因型的框架。化t Genet 43:491-98)。[君合:运个算法 是和化mer的一样吗?我们使用哪一种？]该算法首先识别用于重新对准的区域，其中（i)至 少一个片段包含插入缺失区域，（ii)存在错配碱基或（iii)一种已知的插入缺失的聚集(例 如，从化SNP数据库(单核巧酸多态性数据库），运是一个公共档案，其涵盖由国家生物技术 信息中屯、(NCBI)开发和托管的各种不同物种的普通变异，其中包含了一系列分子的变异， 包括（1)单核巧酸多态性；（2)短缺失和插入多态性，（3)微卫星或短串联重复序列（STR)， (4)复核多态性(MNPs)，（5)杂序列，W及(6)命名的变体。在每个区域中，单倍型是由向参考 序列定点渗入任何已知的插入缺失而生成，片段中的插入缺失来自整个位置点或者来自于 所有不能完美定位参照序列的片段的史密斯-沃特曼比对(Durbin等（1998)生物序列分析： 蛋白质和核酸的概率模型（剑桥大学出版社，剑桥，UK))。对于每一个单倍型出，片段无间距
[0化4] 地和化对齐，W如下标准进行打分：
[0化5]
[0化6]
[0057]其中Rj是第j个片段，k是R神日出无缝对准的偏移，e化是依据片段Rj的第k个碱基的 声明的质量分数而决定的错误率。可W最大化U出)的单倍型出被选为最好的可替代的单倍 型。接着，所有片段根据最好的单倍型和参照化0)重新对准，并且每个片段R非皮分配给出或 化，取决于哪一个可W最大化L(RjlH)。如果指数几率比值或者两单倍型模型比单个基准单 倍型要好至少五个指数单位，那么片段就需要被重新调整：
[0化引
[0059] 此离散化反映了精度和完整统计量的有效计算之间的折衷。在某些例子中，算法 对所有个体中的所有片段同时运行，从而确保了所有个体之中的推断单倍型的一致性，运 对于可靠的插入缺失识别和对比分析，例如躯体SNP和插入缺失识别，是十分关键的。通常 情况下，重新调整的片段被写入到SAM/BAM中W进行进一步分析。
[0060] 局域调整模块的输出随后被馈送到碱基重新校准模块140,运为每个片段中的每 个碱基提供了经验性准确的碱基质量分数。在某些实施例中，碱基重新校准模块140还校正 了误差协变量，例如机器周期和二核巧酸，W及测序平台特定的误差协变量例如SOLiD的颜 色空间不匹配和454的流循环。碱基重新校准模块140的示例性算法已由DePristo MA等人 所述。（DePristo M等(2011)，一种使用二代DNA测序数据进行变异发现和基因型识别的框 架。化t Genet 43:491-98)。通常情况下，该算法首先列出每一个道与在所有已知数量上不 改变(化SNP build 129)的位点的参照的经验不匹配度，根据所报道的质量分数(R)、片段 的机器周甜化)巧^巧巧酸化苯-倍个*别的经验质量分数可W根据如下标准估计：
[0061]
[0062]
[0063] Qempirical(R,C,D)=(mismatch(R,C,D)+l)/(bases(R,C,D)+1)
[0064] 该协变量随后被分解成线性可分的误差估计，重新校正的质量分数Qreeal可W用下 式计算：
[00化]recal(;r,c,d) =Qr+A Q(T)+A A Q(；r,c)+A A Q(；r,d)
[0066]
[0067]
[006引
[0069] A Q(r,d) =Qempiricai(r,C,d)-( A Qr+A Q(r))
[0070] 其中，每个A中和A A A是经验不匹配率和报导的所有仅包含了Qr或者包含了协 变量和化的观察的质量分数之间的剩余差异。其中Qr是碱基的所报导的质量分数和Er为它 预期的错误率;br,e,d是有特定协变量的碱基，且r，c，d和R，C，D分别是一系列的报导的质量 得分，机器周期和二核巧酸。
[0071] 碱基重新校准模块140的输出随后被送入变体识别模块150, W在包括单核巧酸多 态性，短插入缺失和拷贝数变异的片段中发现具有替代等位基因的统计学证据的所有站 点。通常，变体识别模块使用由设及的一个或多个随机变量和可能的其他非随机变量的数 学方程规定的算法或统计模型。例如，根据片段深度和变体计数，表明在特定位置的特定变 体的置信水平为真阳性的概率可W用基于统计模型的方法W及用基准样品的局域化方法 计算出来。
[0072] 已经开发了各种用于变体识别的算法。例如，定位和有质量的组装(MAQKLi H等 人(2008)使用定位质量分数W定位短DNA序列片段和识别变体。Genome Res 18:1851-58) 和SOAPs叩（Li R等人（2009)大规模并行的全基因组重测序的SNP检测。Genome Res 19: 1124-32)使用固定的杂合子和核巧酸-片段错误的先验概率值。SeqEM(Martin ER等人 (2010)SeqEM:二代测序研究的自适应的基因型识别方法Bioinformatics 26:2803-10)引 入了通过自适应方法调用期望最大化化M)算法来估计模型参数的多样品基因型识别。 SAMtoo Is使用修订后的MAQ模型来估计测序错误。该gif too Is家族（gif Single， glfMultipies和polymutt)从预先生成的基因型可能性文件（GLF)中调用SNPeGATK采用 MapReduce的理念，W并行编程简单贝叶斯模型（Dean J，Ghemawat S(2008)MapReduce:大 型集群的简化数据处理Commu ACM 51:107-13) "Atlas2采用已验证的全外显子组捕获测序 数据，而不是常规的可能性计算的物流回归模型，并已被证明具有高灵敏度(Ji HP(2012) 识别外显子变体的改进的生物信息学流水线Genome Med 4:7)。
[0073] 在某些实施方案中，计算流水线可如上所述具有更少的模块。例如，Liu X等人描 述的跳过局域重组模块、重复标记模块和碱基重新校正模块的流水线化iu X等人(2013年） 二代测序数据的变体识别：比较研究化OS 0肥8(9): 675619)。
[0074] 单通多设定的计算流水线
[0075] 在片段通过计算流水线时，用户需要为许多参数设值。例如，在变体识别模块中， 用户需要指定先验概率。先验概率"指的是反映之前变体分布信息的概率分布，可W被用来 进行变体识别。例如，在采用GATK算法的变体识别模块中，用户需要指定全基因组单核巧酸 多态性(SNP)概率，和插入缺失概率。"全基因组SNP概率"是指单个核巧酸在基因组的一些 特定位置发生变异的可能性，其中每个变异有一定明显的程度地存在在族群之中。"插入缺 失概率"，是指在一个生物体DNA中发生插入或者缺失碱基的概率。类似的，在采用SAMTools 算法的变体识别模块中，用户需要指定化red-scaled间隙延伸错误概率和化red-scaled间 隙开口错误概率。
[0076] 用户需要在流水线中设置许多其它参数。例如，一个局域对准模块中，用户需要指 定针对不匹配的惩罚（一个用于对准少量零散基因序列的打分系统，例如，为了更准确地对 准片段，突变被注释为序列中的间隙，该间隙经由各种惩罚计分方法惩罚，从而实现了序列 对准的优化，W获得基于可用的信息的最好的对准），间隙开口（打开任何长度的间隙所需 的成本)和缺口延伸(延长现有间隙的长度所需要的成本）。碱基质量校准时，用户需要指定 校准表(在表格中，为每组协变量建立跟踪匹配/不匹配的数据，并且给在一个识别组中的 每个变体识别分配一个已经校准过的概率）。用户可能还需要指定变种识别模式（即在一组 表示变种频率的估计值和自信指数的数据中最常出现的数值)和变体识别阔值(该程序应 该发出变异位点的最小置信度阔值;如果该点相关的基因型有着比识别阔值更低的信度指 数，程序会发出过滤的站点；运个阔值将高信度指数识别和低信度指数识别分离开来），成 对隐含马尔可夫模型(成对-HMM，基本HMM的一个变体对查找序列比对并评估对准符号的显 著程度尤其有用）W及其中使用的转换概率等。
[0077] 为了更好地评估变种，用户需要更改参数或根据一系列参数设置来研究变种。例 如，用户可能希望看到当SNPAndel的先验概率被改变时，所识别的变体的概率如何改变， 所W要探索样品的各种情况。例如，如果所述样品来自癌细胞，可能需要使用较高的先验概 率。根据样品的族群，需要使用不同的先验概率，因为某些族群可能有比其他族群更高的变 体概率。如果用户之前不知道最佳先验概率，他们可能希望在一定范围的先验概率下探索 变体特性，所W要覆盖所有可能的情况。
[0078] 参数探索的一种方法是再次经历整个计算流水线来设置每个参数。例如，在用户 用最可能的参数设置运行流水线后，并已经审查了被识别的变体，用户可W决定返回改变 参数并重新运行整个变体识别流水线，看看运些识别的变体会有什么改变。整个管道的运 样的迭代重新运行是非常耗时和极度昂贵的，从而导致只能探索有限参数的设置，因此，可 能会导致某些生物学显著信息被错过。
[0079] 因此，在一个方面，本公开内容提供了用计算流水线分析来自多个核酸序列片段 的变体的方法和系统，该计算流水线包括至少一个恒定模块和至少一个可变模块，其中可 变模炔基于可变参数。在某些实施方案中，该方法包括在处理器上执行一系列步骤:接收所 述多个核酸序列片段;为可变参数设定多个数值;将多个核酸序列片段通过至少一个不变 模块一次;将多个核酸序列片段通过至少一个可变模块多次，每次使用可变参数设定的多 个数值中的一个;并且产生多个变体识别。
[0080] 图2展示了一个示例性的计算流水线，其包括至少一个恒定模块和至少一个可变 模块。如图2所示，用于分析核酸序列片段的计算流水线200包括串联连接的模块阵列。
[0081] 原始片段数据被馈送到定位模块210将短测序片段定位到参照。定位模块210的参 数，如对准种子、种子长度和允许错配的最大数量分别被设置了一个值。映射模块210的输 出被馈送到重复标记模块220。重复标记模块220的参数，例如选择非重复片段的评分策略， 和两组重复片段的最大补偿分别被设置了单一数值。重复标记模块220的输出被馈送到局 域重比对模块230。局域调整模块230的参数，包括不匹配惩罚，间隙开口和缺口延伸分别被 设置了单一数值。局域重比对模块230的输出被馈送到碱基重新校准模块240。碱基重新校 准模块240的参数包括重新校准表(上文详细描述过）中的值被分别设置为单一数值。碱基 重新校准模块240的输出被馈送到变体识别模块250。多组数值251,252,253被赋予给各个 参数，例如，先验概率。因而，输入的数据通过变体识别模块250多次，通过将每组数值(251， 252和253)赋予给各个参数。在一个实例中，对一个3 X3矩阵:SNP(0.002,0.001,0.0005) X Indel(0.0002,0.0001,0.00005)进行了探讨，其中，SNP 0.002指SNP先验概率被设定为 0.002,如此等等。因此，该变体识别模块用3x3不同的参数组合运行9次，生产9个变体识别 输出。
[0082] 在某些实施方案中，一些变体识别模块250中的部分计算，如片段过滤，仅取决于 片段碱基的品质和对准，W及碱的某些条件概率，不依赖于现有的概率的，所W只需被运行 一次。
[0083] 通过去除冗余运算，变体识别模块250中使用多个参数设置的多个输出可W靠通 过计算流水线200-次来获得。一次通过需要的总计算量仅仅是多次通过完整流水线的计 算量的一小部分。
[0084] 在运个单次通过多参数设定的运行之后，变体识别模块250的输出可W在视觉上 显示为参数集合的一个函数，使用户对变体有更好的视觉上的认识。
[0085] 逐步勘探计算流水线
[0086] 相对于计算并且输出多个参数在一次通过的运行结果，该流水线可W计算并输出 在一个参数设置下的第一次通过结果。第一次通过时，每个模块的输出，例如，输出中间结 果，被保存下来。用户可W根据该结果，为至少一个流水线的模块，例如，第一次通过的变体 识别，调整参数设定，并使片段第二次通过流水线。在第二次通过中，流水线可W自动识别 不需要重新计算的中间结果，并且直接在新参数设定中使用它们。运样，对于附加的参数设 置仅需要计算那些被新旧参数设置差异影响的流水线的部分，因此，节省了参数勘探中的 计算量。该方法的优势也在于它可W避免对一些不必要的参数设定所进行的探索。
[0087] 因此，在另一个方面，本公开内容提供了使用包括至少一个恒定模块和至少一个 可变模块的计算流水线，来分析多个核算片段中的变体的方法与系统，其中可变模炔基于 可变参数。在某些实施方案中，该方法包括在处理器上执行一系列步骤:接收多个核酸序列 片段，其中至少一个核酸片段包括一个变体;为可变参数设置第一数值;让多个核酸片段通 过至少一个恒定模块一次;通过至少一个可变模块，其中使用可变参数的第一数值;产生第 一变体识别；根据第一变体识别为可变参数设置第二数值;让多个核酸片段第二次通过至 少一个可变模块，其中使用可变参数的第二数值;并且在不通过至少一个恒定模块的情况 下，产生第二个变体识别。
[0088] 图3展示了用于分析核酸序列片段的示例性增量勘探计算流水线300。计算流水线 300包括串联连接的模块阵列。如上文所详细描述地，当原始片段数据被馈送到计算流水线 300时，它会通过定位模块310，重复标记模块320，局域重比对模块330，碱基重新校准模块 340和变体识别模块350。第一个数值集合被用于变体识别模块350中的先验概率，并且获得 第一个变体识别的结果。在第一个变体识别的结果的基础上，第二组数值被设置为变体识 别模块350中的概率。碱基重新校准340的输出被馈送到变体识别模块350, W产生一个第二 个变体识别结果。
[0089] 在某些实施方案中，"单通多设置"模式和"增量勘探"模式可W结合起来。例如，流 水线在第一个单次通过中运行参数设置，产生多个变体识别，同时还保存了中间值的输出。 在审查完在第一单次通过中产生的多个变体识别后，用户可W调整一个或多个模块的参 数。然后流水线可W将第一通中的中间结果，重新用在新的通过中，从而节省了在新的通过 中需要的计算。
[0090] 定位模块
[0091] 在另一个方面，本公开提供了改进对比参照计算的方法和系统。在定位过程中，用 户也可W选择一些参数，例如，对准种子，种子长度，允许的最大错配数量，等等。当用户探 索对准参数时，已经与基准完美匹配的片段不需要在每次参数设定中重复被对准，而是可 W只运行一次，其结果可W被用在所有情况下。只有没有定位的片段W及与参照有错配的 片段才需要在不同的参数制定下，反复多次对准。
[009。组装流水线
[0093] 在另一方面，本公开提供了从头组装核酸片段的方法和系统，例如不需借助参照 (即，参照基因组序列）来进行组装。从头组件通常被用于研究非模式生物体，因为他们的参 照基因组序列常常不存在。从头组件还被用于转录序列数据。一些从头组装的程序已经被 开发出来，包括 SOAPdeno VO-Trans，Velvet/0ases，^ans-ABySS,和IYinity。两种一般被用 于从头装配的基本算法是:重叠图形和德布鲁因图。虽然重叠图大多用于桑格测序片段，因 为它比德布鲁因图具有更大的计算强度，并且W最高的效率装配更少的高度重叠片段。德 布鲁因图基于k-1的序列守恒对准K-mers(通常为25-50bp)W创建重叠群(连续，重复序列 片段）。对于比片段长度更短的K-mers的使用，使得德布鲁因图表法降低了计算的强度，并 且适用于二代测序数据。
[0094] 在某些实施方案中，从头装配的计算流水线包括至少一个恒定模块，例如，片段清 洁和过滤模块，和至少一个可变模块，例如，重叠群模块。在通过组装流水线时，多个值被赋 予给可变模块的可变参数，例如，重叠群模块的K-mers，同时单值被赋给不变模块的恒定参 数。因此，不同参数设置之间的共同计算可W被节省下来。
[0095] 本发明可W用如下实施例进行说明。运些实施例不旨在从任何方面限制本发明。
[0096] 实施例1
[0097] W下是一个用于分析核酸序列片段的单通多设置的计算流水线的例子。
[0098] 该计算流水线由W下模块串联而成：定位模块，重复标记模块，局域重新比对模 块，碱基质量重新校准模块和变体识别模块。
[0099] 计算流水线的实施使用SuperServer 6016T-NTF服务器（Intel Dual Xeon E5620 2.4 Hz CPU,48 GB DDR3 1333MHz Ecc Memory,和2 TB SATA Seagate ST32000644NS 皿D)。
[0100] 用于分析的核酸序列片段来自人类基因组外显子组测序数据。
[0101]具有化StQ格式的核酸片段（由Illumina公司系统生产)被馈送到定位模块，该定 位模块用BWA算法将片段对应到参照基因组化uman glk v37) dSAM格式的定位文件被转换 为BAM格式并且由GATK进行排序。排序后的BAM文件被馈送到重复标记模块，采用Picard算 法删除重复的PCR。去掉重复部分的BAM文件被输送到使用GATK算法的局域重新比对模块W 及作为基准的变量识别格式(VCf)的化SNP。重新对准的BAM文件被馈送到使用GATK算法的 碱基质量重新校准模块。
[0102] 上述模块中所用参数的具体数值设置如下:重复标记模块中的两组重复片段的最 大补偿为100，局部重比对模块的LOD阔值为5.0。
[0103] 重新校准的BAM文件随后被馈送到使用GATK算法的变体识别模块。运个模块被运 行了9次 W探索3 X 3矩阵:SNP(0.002，0.0 Ol，0.0005) X 插入缺失（Indel) (0.0002，0.0 OOl， 0.00005)O
[0104] 结果:大约需要2300分钟来运行上述单通多设置计算流水线（定位模块大约需要 40分钟，重复片段标记模块大约需要40分钟，局部重比对模块大约需要40分钟，碱基质量重 新校准模块需要大约180分钟，变体识别模块大约需要2000分钟）。作为比较，如果使用单设 置的计算流水线运行9次W探索上述3X3矩阵则需要大约4800分钟。通过使用变体识别模 块的矩阵中重新校准的BAM文件，上述流水线中所使用的总计算量与通过全流水线9次W探 讨矩阵所需计算量相比，仅仅是很小的一部分。
[0105]虽然本发明已经被具体地展示，并依据具体实施方案给予了详细说明（其中部分 是优选的实施方案），但需要被本领域技术人员明白的是，在不脱离本文公开的发明的精神 和范围的情况下，依然可W对其做出各种形式和细节上的调整。
【主权项】
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所述的方法，其特征在于，所述恒定模炔基于恒定参数，并且所述方 法包括为所述恒定参数设置数值的步骤。
【文档编号】G06F19/22GK106021998SQ201610173000
【公开日】2016年10月12日
【申请日】2016年3月24日
【发明人】叶军, 周巍, 陈洛祁, 冯汉鹰, 陈洪, 刘晓峰
【申请人】知源生信公司（美国硅谷）