位置相关变体识别计算流水线的制作方法

位置相关变体识别计算流水线的制作方法
【专利摘要】本公开提供了一种利用计算机辅助实现的分析多个核酸序列片段中变体的方法。该方法使用的计算流水线包括使用位置相关参数的模块。该方法包括在处理器上执行下列步骤：接收多个核酸序列片段，这些核酸序列片段至少包括一个第一核酸序列片段和一个第二核酸序列片段；分别将第一核酸序列片段和第二核酸序列片段比对到基因组中的第一位置和第二位置；基于第一位置和第二位置，分别向该位置相关参数赋予第一数值和第二数值；将第一核酸序列片段和第二核酸序列片段输送通过使用位置相关参数的模块，分别使用第一数值和使用第二数值；并且产生变体识别。
【专利说明】位置相关变体识别计算流水线
[0001] 相关申请的交叉引用
[0002] 本申请要求2015年3月27日提交的美国临时申请62/139,148和2015年4月3日提交 的美国临时申请62/143,013的优先权。其全部内容在此参考并入。 发明领域
[0003] 本发明主要设及基因测序数据的分析。
【背景技术】
[0004] 二代测序技术(NGS)为大量生产生物数据提供了强有力的工具，并为取得个性化 医疗提供了帮助。虽然仅就取得序列数据来说，高通量基因测序的成本有所降低，但是分析 与解读运些大规模测序数据依然存在巨大的挑战。为了识别NGS数据中的变体，大量序列比 对器和变体识别器被研发出来并且被整合成各式各样的计算流水线。一个典型的计算流水 线包含一个序列比对器和一个变体识别器:前者可W将序列片段与参考基因组进行比对， 后者确定变异点并且向对象分配一个基因型。在计算过程中，用户为了正确分析序列数据， 常常需要设定许多参数。更重要的是，一些数据需要基于细胞种类或者用于制备样品的族 群来进行优化，从而准确识别变体。具体来说，一个参数的最佳取值可能取决于其在基因组 中的位置。例如，与变体出现概率相关的参数就可能取决于其在基因组中的位置。然而由于 计算流水线每次运行都需要巨大的计算量，通过运行整个变体识别计算流水线来测试每一 个参数设定在实际操作中是很难实现的。因此，有持续的需求研发新的方法和系统来针对 基因组位置进行分析NGS数据的参数优化。
[00化]发明概述
[0006] 一方面，本发明提供了一种通过用电脑执行计算流水线来分析多个核酸序列片段 中的变体的方法。该方法包括基于位置相关参数的模块。该方法包括用处理器执行一系列 步骤:获取多个核酸序列片段，其至少包含第一核酸序列片段和一个第二核酸序列片段;将 第一核酸序列片段比对到基因组的第一位置;将第二核酸序列片段比对到基因组的第二位 置，其中第一位置与第二位置不相同。根据基因组第一位置为位置相关参数设定第一数值； 根据基因组第二位置为位置相关参数设定第二数值，其中第一数值与第二数值不相同；将 第一核酸序列片段输送通过所述模块，其中使用第一数值;将第二核酸序列片段输送通过 所述模块，其中使用第二数值;并且生成变体识别。
[0007] 在某些实施方式中，所述模块是变体识别模块。在某些实施方式中，所述位置相关 参数是变体识别的先验概率或者阔值。在某些实施方式中，所述先验概率是全基因组单核 巧酸多态性概率，或者插入缺失概率。
[000引在某些实施方式中，所述模块是局域比对模块。在某些实施方式中，所述位置相关 参数是碱基错配惩罚，间隙开口，间隙延伸，或者一个比对候选者的阔值。
[0009]在某些实施方式中，所述第一和/或第二位置是外显子或内含子。在某些实施方式 中，所述第一和/或第二位置是细胞核周边或细胞核中屯、。
[0010] 在某些实施方式中，所述基因组从一个族群或者区域群体中获得。在某些实施方 式中，所述基因组从一个健康主体或者一个患有疾病的主体中获得。在某些实施方式中，所 述疾病是癌症。
[0011] 另一方面，本发明提供了一种非暂时性计算机可读介质W及通过使用其存储的计 算流水线来识别来自多个核酸序列片段的变体的方法。该计算流水线包括至少一个基于位 置相关参数的模块。处理器执行的指令行实现一系列步骤，包括:获取多个核酸序列片段， 其至少包括一个第一核酸序列片段和一个第二核酸序列片段;将第一核酸序列片段比对到 基因组中的第一位置;将第二核酸序列片段比对到基因组中的第二位置，其中第一位置不 同于第二位置;基于基因组中的第一位置为所述位置相关参数赋予第一数值;基于基因组 中的第二位置，为所述位置相关参数赋予第二数值，其中第一数值不同于第二数值;将第一 核酸序列片段输送通过所述模块，其中使用第一数值;将第二核酸序列片段输送通过所述 模块，其中使用第二数值;并且生成变体识别。
[0012] 另一方面，本发明提供了一种通过用电脑执行计算流水线来分析大量核酸序列片 段中的变体的方法。该方法包括定位模块和变体识别模块，其中该变体识别模块取决于先 验概率参数。所述方法包括在处理器上执行一系列步骤:获取多个核酸序列片段，其至少包 含第一核酸序列片段和第二核酸序列片段的;将所述第一核酸序列片段比对到基因组中的 第一位置;将所述第二核酸序列片段比对到基因组中的第二位置，其中第一位置不同于第 二位置;基于基因组中的第一位置为所述先验概率参数赋予第一数值;基于基因组中的第 二位置，为先验概率参数赋予第二数值，其中第一数值不同于第二数值;将第一核酸序列片 段输送通过所述变体识别模块，其中使用第一数值;将第二核酸序列片段输送通过所述变 体识别模块，其中使用第二数值;并且生成变体识别。
[0013] 另一方面，本发明提供了一种非暂时性计算机可读介质W及通过使用其存储的计 算流水线来识别来自多个核酸序列片段的变体的方法。该计算流水线包括比对模块和变体 识别模块，其中变体识别模块依赖于先验概率参数。处理器执行的指令行实现一系列步骤， 包括:获取多个核酸序列片段，其至少包括第一核酸序列片段和第二核酸序列片段;将第一 核酸序列片段比对到基因组中的第一位置进行比对;将第二核酸序列片段比对到基因组中 的第二位置，其中第一位置不同于第二位置;基于基因组中的第一位置为先验概率参数赋 予第一数值;基于基因组中的第二位置，为先验概率参数赋予第二数值，其中第一数值不同 于第二数值;将第一核酸序列片段输送通过所述变体识别模块，其中使用第一数值;将第二 核酸序列片段输送通过所述变体识别模块，其中使用第二数值;并且生成变体识别。
[0014] 该发明的上述特征优势可W通过W下的描述、权利要求和附图，得到更好的理解。 [001引附图简要说明
[0016] 图1.展示了一个示例性计算流水线。
[0017] 图2.展示了一个包含了一个使用位置相关参数的模块的示例性计算流水线。
[001引发明详述
[0019]在关于W上发明的概述、具体描述、如下的权利要求W及附图中，引用了本发明中 的特定特征(包括步骤）。应当理解的是本发明的说明书中包含了对运些特定特征的所有可 能的组合。例如，当发明的实施例或者一个特定方面或者一个权利要求中展示了一个特定 的特征，该特征也可W在可能的程度上在本发明的其他方面或实施例中使用。
[0020] 应当理解，除非根据上下文不允许，本说明书和权利要求中使用的单数形式"一 个"包括复数形式。比如，一个"模块"包括一个或多个模块，W及本领域技术人员知道的等 价形式。
[0021] -个具有两个或者更多特定步骤的方法可W被W任意顺序或者同时执行(除非在 文本中排除了运种可能性）。该方法可W包含一个或多个其他步骤，运些步骤可W在任意特 定步骤前面，或者在两个特定步骤之间，或者在所有的特定步骤之后被执行(除非文本中排 除了运种可能性）。
[0022] 在提供的值的范围内，可W理解，每个居中值，到下限的单位的十分之一，除非上 下文清楚地另有规定，该范围的上限和下限和任何在该所述范围中的所述或者居中值，都 被本公开内容涵盖，同时符合在所述范围内明确排除的极限。当所述范围包括一个或者两 个极限，排除运两个极限或其中一个极限的范围也被包含在运个公开内容中。
[0023] 为了简单清楚的阐述，当合适的时候，标号在不同的附图中重复使用，W指示相应 的或类似的元件。此外，大量的具体细节被提供，W便透彻理解运里所描述的实施例的阐 述。然而，本文描述的实施例可W在不存在具体细节的情况下实施。在其他实例中，方法、程 序和组件没有详细描述，但没有模糊所描述的相关功能。此外，描述不应被认为是限制本文 所述的实施方式的范围。应该理解的是，除非另有说明，在本公开中阐述的实施例的描述和 表征并非相互排斥。
[0024] 定义
[0025] 本公开内容使用了如下定义：
[0026] 术语"包括"W及它们在本文中的同义词代指其它组分，成分，步骤等是任选存在。 例如，物品"包括"（或"其包含"）组分A，B和C可W由（即只包含)组分A，B和C，或可W不仅包 含组分A，B，和C，还有一种或多种其它组分。
[0027] 后跟数字的术语"至少"在本文中用于表示W该数字开始的一个范围的开始(可W 是具有上限或没有上限的范围，运取决于所限定的变量）。例如，"至少1"表示1或大于1。后 跟数字的术语"至多"在本文中用于表示一个W该数字结尾的范围（它可W是具有1或0作为 其下限的范围，也可W是不具有下限的范围，运取决于被定义的变量）。例如，"至多4"是指4 或者小于4, W及"至多40%"是指40%或低于40%。在本公开中，当一个范围被给定为"(第 一数字)至(第二个数字r或"(第一数字)-(第二数字)"，运表示范围的下限是第一数字，上 限为第二数字。例如，25至IOOmm表示下限为25毫米，上限为100毫米的范围。
[0028] 如本文所用，术语"核酸序列片段"指的是由测序方法确定的核酸序列。该核酸序 列可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。
[0029] 术语"变体"当在核酸序列的语境中使用时是指一种与参照不同的核酸序列片段。 典型的核酸序列变体包括但不限于单核巧酸多态性(SNP)，短缺失插入多态性（Indel),拷 贝数变异(CNV)，微卫星标记或短串联重复序列和结构的变化。
[0030] 如本文所使用的，术语"计算机实现的方法"是指相关的方法是在计算机中执行， 例如，一个由CPU执行的计算机程序。一种计算机，如本文所使用的，指的是可被编程W自动 执行的一组算术或逻辑运算的设备(为了一般或特定目的）。计算机，如在此使用的，包括但 不限于个人计算机，工作站，服务器，大型机和超级计算机。计算机可W是独立的系统，网络 系统或位于计算云中的虚拟机。本公开中的方法可W通过多线程或其他并行计算方式实 现。
[0031] 如本文所使用的，"计算流水线"或"流水线"是指一组串联连接的数据处理元件， 其中一个元件的输出是下一个元件的输入。在某些实施方案中，一个操作的输出被自动馈 送到下一个操作。如本文所使用的，计算流水线中的元素可W被称为"模块"。在某些实施方 案中，流水线是线性的，单向的。在某些实施方案中，主要是单向的流水线可W在其他方向 上的有一些交流。在某些实施方案中，流水线可W是完全双向的。
[0032] 如本文所使用的，"模块"指的是一个在计算流水线内的数据处理元件。一组模块 W串联形式连接，从而形成一个计算流水线。通常，一个模块接收一个输入数据，基于所输 入的数据执行特定的功能，并产生输出数据，然后将其用于随后模块的输入数据。在某些实 施方案中，一个模块可W被进一步分成几个子模块，例如，子模块可W W串联形式连接。
[0033] 术语"参数"，如本文中所使用的，指的是用户需要在计算流水线或模块中设定的 基准，特征或数值。当数据集通过计算流水线或模块时，基准，功能或数值被传递给函数，程 序，子例程，指令，或程序。相关术语"位置相关参数"被使用在将核酸序列片段输送通过一 个模块的语境中时，指的是参数的数值随着核酸序列片段在参考(例如，参考基因组）中的 位置而变化。
[0034] 术语"数值"，如本文中所使用的，是指参数设定的具体数或特征。相应的，当使用 计算流水线分析核酸序列片段时，该具体的数或特征被传递给计算流水线的函数，程序，子 例程，指令，或程序。为某个参数设定的具体数或特征可W根据所讨论的模块和参数确定。 例如，在变体识别模块中使用的先验概率参数的数值可W是0.0005,0.0008,0.OOl或 0.002。又例如，在比对模块中使用的错配惩罚参数的数值可W是+1到巧。
[0035] 在"使数据(例如，核酸序列片段)通过流水线或一个模块"中所使用的术语"通过" 指的是用流水线或者模块分析所述数据。典型的，在通过一个流水线或模块时，所述数据被 作为输入数据输送到流水线或模块。该流水线或模块使用该数据运行流水线或模块中的函 数，程序，子例程，指令，或程序并产生输出数据。在某些实施方式中，所产生的输出数据作 为输入数据输送到另一流水线或模块中。在某些实施方案中，使数据通过一个模块还包括 用所述数据作为该模块的间接输入。例如，在一个包括至少两个W串联连接的模块的计算 流水线中，所述数据作为第一模块的输入被传递，其生成的输出又作为所述第二模块的输 入被传递。在运样的情况下，所述数据被视为通过第二模块，尽管他们并没有被用作所述第 二模块的直接输入。
[0036] 术语"非暂时性计算机可读介质"指的是任何计算机可读介质，唯一的例外是一个 短暂的传播信号。非暂时性计算机可读介质包括，但不限于，易失性存储器，非易失性存储 器，软盘，硬盘，存储棒，寄存器存储器，处理器高速缓存和随机存取存储器。
[0037] 本文所用术语"定位"或"定位到参照"是指将核酸序列片段与参照，例如，序列已 知的参考基因组，比对。各种程序和算法已经被开发来将核酸序列片段定位到一个参照(参 见，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)二代测序哦按段对准算法的对比分 析。Bioinformatics27:2790-96;Patna化S等（2012)使用组合方法的外显子组数据分析流 水线的定制化OS 0肥7:e30080)。在各种程序和算法中，基于己路±惠勒改造的己路±惠 勒比对法（BWA)，化i H，Durbin R(2009年）快速，准确的己路±惠勒短弧度排列。 Bioinformatics 25:1754-60)体现了运行时间，存储器的使用和准确性之间的良好平衡， 并常常被使用在不同的计算流水线中。 陶]分析核酸序列片段的计算流水线
[0039] 二代测序技术的迅速发展在过去几年变革了生物和生物医学研究。根据所使用的 测序方法和系统，产生的核酸片段的数量通常在数百万W上。例如，Illumina的MiniSeq系 统每次运行产生高达2500万的片段，而Illumina公司的化Seq系列每次运行可W产生高达 50亿片段。随着大量测序数据的生成，对可W用来分析和解释运些大规模测序数据的强大 的计算工具的需求迫在眉睫。
[0040] 已经开发了许多具有多个序列比对器和多个变体识别器的计算流水线，其中包 括，但不限于，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)。
[0041 ]用于分析核酸序列片段的示例性计算流水线如图1所示。图1中，用于分析核酸序 列片段的计算流水线100包括串联连接的模块阵列。
[0042]首先，原始读取的数据被馈送到定位模块110来使短测序片段与参照比对对。定位 模块110将核酸序列片段与参照比对，例如，一个序列已知的参考基因组。各种程序和算法 已经被开发出来将核酸序列片段定位到一个参照（参见，Flicek P，Birney E. (2009)从序 列片段中检测：用于对齐和装配的方法，化tMethods 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)，展现了运行时间， 存储器的使用情况和准确性之间的良好平衡，并被频繁使用在不同的计算流水线中。
[00创比对模块的输出（例如，SAM(序列比/映射)文件或BAM(SAM的二进制版本)文件)被 馈送到重复标记模块120, W去除PCT重复。在制备DNA测序样品的过程中，PCR经常被用来扩 增所述片段，从而产生复制品。理想的情况下，制备的样品常常产生百分之几(例如，约4%) 的彼此完全一样的片段拷贝，即，复制品。有时，30%至70%的片段是重复的。Wysoker A等 (PicardTools)和Li H等化i H等（2009)序列比对/映射格式和SAMtools,Bioinformatics 25:2078-79)。
[0044] 复制被标记或者移除的片段随后被送入到局部重比对模块130, W改善片段的比 对。通常情况下，相对于基准来说，局部重新比对发生在片段插入和缺失（Indel)的区域周 围，并且是将片段与Indel-侧的一端相对齐，然后再对另一侧的其余部分进行定位。当片 段最初定位到参照时，并不能获得任何关于插入缺失的信息。因此，定位到运样的区域的片 段，只具有一小段代表插入缺失一侧的区域，并通常不会被正确定位整个indel，而是会有 一端没有对齐，或者整个indel区域有很多不匹配的定位。局部重比对模块130使用其余定 位到含有indel区域的片段信息，包括位于插入缺失区域更居中位置的片段，并且因此已经 与所述插入缺失的两侧端部对齐。其结果是，生成了另一个和之前一样好甚至更好的定位。
[0045] 局部重比对算法已由化mer等人描述。巧omer N(2010)用srma通过局域重新调二 代测序数据短片段来改进变体发现。Genome Bioll 1(10):R99)。在第一步骤中，所有输入 片段的对齐信息被整理在一个高效的基于图的数据结构中，其基本上类似于de-BruUn的 图表。运种调整图表示了片段是如何对齐于参考序列W及片段是如何彼此重叠。在第二步 骤中，元数据是从可W反映出可W潜在地改善片段定位的对准位置的图表结构中获取，该 图表结构还提供了如何重新比对片段W得到最简洁多个对准的假设。在第=步骤中，重新 比对图和其元数据被用于实际执行各个片段的局域比对。
[0046] DePristo等描述了用于局域重新比对的替代算法(DePristo MA等（2011)。一种使 用二代DNA测序数据的变体发现和基因型的框架。化t Genet 43:491-98)。该算法首先识别 用于重新对准的区域，其中（i)至少一个片段包含插入缺失区域，（ii)存在错配碱基或 (iii) 一种已知的插入缺失的聚集(例如，从化SNP数据库(单核巧酸多态性数据库），运是一 个公共档案，其涵盖由国家生物技术信息中屯、(NCBI)开发和托管的各种不同物种的普通变 异，其中包含了一系列分子的变异，包括（1)单核巧酸多态性；（2)短缺失和插入多态性，（3) 微卫星或短串联重复序列（STR)，（4)复核多态性(MNPs)，（5)杂序列，W及(6)命名的变体。 在每个区域中，单倍型是由向参考序列定点渗入任何已知的插入缺失而生成，片段中的插 入缺失来自整个位置点或者来自于所有不能完美定位参照序列的片段的史密斯-沃特曼比 对(Durbin等（1998)生物序列分析：蛋白质和核酸的概率模型（剑桥大学出版社，剑桥， UK))。对于每一个单倍型出，片段无间距地和化对齐，W如下标准进行打分：
[0047]
[004引其中Rj是第j个片段，k是R神日出无缝对准的偏移，e化是依据片段Rj的第k个碱基的 声明的质量分数而决定的错误率。可W最大化U出)的单倍型出被选为最好的可替代的单倍 型。接着，所有片段根据最好的单倍型和参照化0)重新对准，并且每个片段R非皮分配给出或 化，取决于哪一个可W最大化L(RjlH)。如果指数几率比值或者两单倍型模型比单个基准单 倍型要好至少五个指数单位，那么片段就需要被重新调整：
[0049]
[0050] 此离散化反映了精度和完整统计量的有效计算之间的折衷。在某些例子中，算法 对所有个体中的所有片段同时运行，从而确保了所有个体之中的推断单倍型的一致性，运 对于可靠的插入缺失识别和对比分析，例如躯体SNP和插入缺失识别，是十分关键的。通常 情况下，重新调整的片段被写入到SAM/BAM中W进行进一步分析。
[0051] 局域调整模块的输出随后被馈送到碱基质量重校准模块140,运为每个片段中的 每个碱基提供了经验性准确的碱基质量分数。在某些实施例中，碱基质量重校准模块140还 校正了误差协变量，例如机器周期和二核巧酸，W及测序平台特定的误差协变量例如SOLiD 的颜色空间不匹配和454的流循环。碱基质量重校准模块140的示例性算法已由DePristo MA等人所述。（DePristo MA等(2011)，一种使用二代DNA测序数据进行变异发现和基因型识 别的框架。化t Genet 43:491-98)。通常情况下，该算法首先列出每一个道与在所有已知数 量上不改变（dbSNP build 129)的位点的参照的经验不匹配度，根据所报道的质量分数 (R)、片段的机器周期(C)和二核巧酸化)分类。每个类别的经验质量分数可W根据如下标准 估计：
[0化2]
[0化3]
[0054] Qempiricai (R, C, D) = (mismatch(R,C,D)+l)/(bases (R,C,D)+1)
[0055] 该协变量随后被分解成线性可分的误差估计，重新校正的质量分数Qreeal可W用下 式计算：
[0化6]
[0化7]
[0化引 [0化9]
[0060]
[0061] 其中，每个A中和A A A是经验不匹配率和报导的所有仅包含了Qr或者包含了协 变量和化的观察的质量分数之间的剩余差异。其中Qr是碱基的所报导的质量分数和Er为它 预期的错误率;br，E，d是有特定协变量的碱基，且r，c，d和R，C，D分别是一系列的报导的质量 得分，机器周期和二核巧酸。
[0062] 碱基质量重校准模块140的输出随后被送入变体识别模块150, W在包括单核巧酸 多态性，短插入缺失和拷贝数变异的片段中发现具有替代等位基因的统计学证据的所有站 点。通常，变体识别模块使用由设及的一个或多个随机变量和可能的其他非随机变量的数 学方程规定的算法或统计模型。例如，根据片段深度和变体计数，表明在特定位置的特定变 体的置信水平为真阳性的概率可W用基于统计模型的方法W及用基准样品的局域化方法 计算出来。
[0063] 已经开发了各种用于变体识别的算法。例如，定位和有质量的组装(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)。
[0064] 在某些实施方案中，计算流水线可如上所述具有更少的模块。例如，Liu X等人描 述的跳过局域重组模块、重复标记模块和碱基重新校正模块的流水线化iu X等人(2013年） 二代测序数据的变体识别：比较研究化OS 0肥8(9): 675619)。
[0065] 具有位置相关设定的计算流水线
[0066] 在片段通过计算流水线时，用户需要为许多参数设值。例如，在变体识别模块中， 用户需要指定先验概率。先验概率"指的是反映之前变体分布信息的概率分布，可W被用来 进行变体识别。例如，在采用GATK算法的变体识别模块中，用户需要指定全基因组单核巧酸 多态性(SNP)概率，和插入缺失概率。"全基因组SNP概率"是指单个核巧酸在基因组的一些 特定位置发生变异的可能性，其中每个变异有一定明显的程度地存在在族群之中。"插入缺 失概率"，是指在一个生物体DNA中发生插入或者缺失碱基的概率。类似的，在采用SAMTools 算法的变体识别模块中，用户需要指定化red-scaled间隙延伸错误概率和化red-scaled间 隙开口错误概率。
[0067] 用户需要在流水线中设置许多其它参数。例如，一个局域比对模块中，用户需要指 定针对不匹配的惩罚（一个用于对准少量零散基因序列的打分系统，例如，为了更准确地对 准片段，突变被注释为序列中的间隙，该间隙经由各种惩罚计分方法惩罚，从而实现了序列 对准的优化，W获得基于可用的信息的最好的对准），间隙开口（打开任何长度的间隙所需 的成本)和缺口延伸(延长现有间隙的长度所需要的成本）。碱基质量校准时，用户需要指定 校准表(在表格中，为每组协变量建立跟踪匹配/不匹配的数据，并且给在一个识别组中的 每个变体识别分配一个已经校准过的概率）。用户可能还需要指定变种识别模式（即在一组 表示变种频率的估计值和自信指数的数据中最常出现的数值)和变体识别阔值(该程序应 该发出变异位点的最小置信度阔值;如果该点相关的基因型有着比识别阔值更低的信度指 数，程序会发出过滤的站点；运个阔值将高信度指数识别和低信度指数识别分离开来），成 对隐含马尔可夫模型(成对-HMM，基本HMM的一个变体对查找序列比对并评估对准符号的显 著程度尤其有用）W及其中使用的转换概率等。
[0068] 对于某些参数来说，他们的最佳数值的选取取决于在基因组中的位置。运些参数 包括，但并不限制于:先验概率;错配惩罚；间隙开口；局域对准中的间隙延伸和成对HMM;变 体识别的阔值;等等。运些参数和变体出现的概率相关，并且取决于在基因组中的位置。
[0069] 例如，已知变体在外显子区域中出现的概率比在内含子区域中要低(Zhao et al. (2003)Investigating single nucleotide poIymorphism(SNP)density in the human genome and its implications for molecular evolution .Gene 312:207-13)。变体出现 频率会因为在不同染色体和染色体区域而有所不同。因此，当给位置相关参数赋值时应该 考虑到基因组中的位置运一因素。例如，对于单核巧酸多态性和缺失插入多态性中的任何 一个或者两者而言，先验概率的一个更为准确的描述是一个基于参考基因组的地图，例如， 一个先验概率的地图描述了整个基因组中，变体在不同位置的出现频率。在某些实施例中， 运样一个可W描述参数在整个基因组中最佳取值的地图被称为参数地图。
[0070] 在某些实施例中，参数地图可W根据一个族群来具体设置，例如，白种人、亚洲人 等等。在某些实施例中，参数地图可W根据一个国家或者一个区域的人口来具体设置，例 如，日本、英国等等。
[0071] 在某些实施例中，参数地图可W根据患有一种特定疾病或者素乱的人口来具体设 置。例如，参数地图可W根据一种特定癌症的患者，例如，乳腺癌、肝癌患者。
[0072] 在某些实施例中，参数地图可W根据特定细胞种类来具体设置，例如皮肤细胞、血 细胞等等。
[0073] 在某些实施例中，先验概率的参数地图可W通过分析大量基因组W获得变体出现 频率来获得，例如在整个基因组中所有位置的单核巧酸多态性或者插入缺失概率。在某些 实施例中，一个族群例如白种人、亚洲人等的参数地图可W通过分析该族群的基因组来获 得。同样地，对于一个特定细胞种类或者患有特定疾病/素乱的群体的参数地图可W通过分 析运个患病/素乱群体的细胞中的基因组来获得。
[0074] 因此，一方面，本发明提供了一个使用计算流水线来分析多个核酸序列片段中的 变体的方法和系统，该系统包括使用位置相关参数的模块。该方法包括在处理器上实施一 系列步骤:获取多个核酸序列片段，其至少包含第一核酸序列片段和第二核酸序列片段;将 第一核酸序列片段比对到基因组的第一位置;将第二核酸序列片段比对到基因组的第二位 置，其中第一位置与第二位置不相同；根据基因组第一位置为位置相关参数设定第一数值； 根据基因组第二位置为位置相关参数设定第二数值，其中第一数值与第二数值不相同；将 第一核酸序列片段输送通过所述模块，其中使用第一数值;将第二核酸序列片段输送通过 所述模块，其中使用第二数值;并且生成变体识别。
[0075] 图2展示了一个示例性的包括使用位置相关参数的模块的计算流水线。如图2所 示，一个用于分析核酸序列片段的计算流水线包括一列串联的模块。
[0076] 包括片段1和片段2的原始数据被馈送到比对模块210,用W使短序列片段向参考 基因组对准。片段1被定位到参考基因组的位置1(例如，外显子），片段2被定位到参考基因 组的位置2(例如，内含子）。比对模块210的输出被馈送到重复标记模块220,其输出被馈送 到局域重新比对模块230。局域比对模块230的输出被亏送到碱基重新校正模块240。碱基重 新校正模块240的输出被馈送到变体识别模块250。当输入数据经过变体识别模块250时，参 数，例如先验概率，基于片段的位置被赋予一个数值。如图所示，当片段1通过变体识别模块 250时，先验概率被赋予第一数值(例如外显子区域中单核巧酸多态概率为0.0005(见Zhao et al.(2003)Investigating single nucleotide poIymorphism(SNP)density in the human genome and its implications for molecular evolution.Gene 312:207-13)))。 当片段2通过变体识别模块250时，先验概率被赋予第二数值(例如内含子区域中单核巧酸 多态概率为0.0008)。并且，产生变体识别。
[0077] 在一个典型的目前正在使用的变体识别流水线中，单一一组先验概率(例如，变体 频率的整体平均值)被用于整个基因组。因而，在先验概率低的区域(例如，外显子内）一些 额外的非必要变体被识别出来，在高先验概率区域(例如，在内含子内），或遗漏一些变体识 另IJ。尽管可W调整先验概率参数，但只要是整个基因组用单一一组数值，就需要牺牲部分区 域来使其他区域获益，而不可能在整个基因组内同时达到优化识别。
[0078] 如上述例子，通过根据片段在基因组中的位置，使用不同的先验概率，对每个变体 识别都可W使用基于其在基因组位置，例如外显子或内含子，的最佳先验概率来实现。运样 一个位置相关的变体识别大大提高了变体识别的准确性。
[0079] 更进一步地，如果已知一个新基因组的种类，例如，种族，W及/或者它是否是癌症 样品和他的细胞种类，计算流水线可W选用一个不同的参数地图。运个选择过程可W被自 动实现:流水线可W计算并且输出使用常规参数数值的第一次通过产生的结果。根据第一 次通过产生的结果，特定的额外信息，例如族群，可W被确定下来。根据第一次通过所获得 的信息确定一个参数地图，使用该参数地图再次运行流水线W提高变体识别的整体准确 度。在必要情况下，该过程可W重复多次。
[0080] 通过使用上述位置相关的变体识别方法，一个多通迭代过程可W被用来改善参数 地图。对于一个特定类型的基因组(族群，癌症细胞种类等），一组该类型基因组数据被输送 通过流水线，并且该流水线在第一次通过中使用一组常规参数。初始先验概率的参数地图 可W通过对在基因组各个位置的数据组中的变体出现频率取平均值来获得。整个数据组随 后被第二次输送通过一个使用一通结果产生的参数地图的流水线，从而提高变体识别的准 确率。然后根据二通的结果可W生成第二个先验概率的参数地图。运样生成的第二个先验 概率的参数地图又可W被应用到第=次通过。运个迭代过程可W-直继续下去直到参数地 图收敛，即第n次和(n-1)次通过的参数地图的差别小于一个特定阔值。运样收敛的参数地 图对于该基因组数据集是最优的。
[0081] 除了先验概率，其余参数的最优值可能也会取决于基因组的位置，例如，错配惩 罚、间隙开口 W及用于局部比对的间隙延伸。因为运些惩罚数值与单核巧酸多态性和插入 缺失多态性的出现概率相关，所W可W据此调整W产生他们对应的参数地图。其余最优值 取决于基因组位置的参数还包括变体识别的阔值。与先验概率参数相同，通过使用运些参 数的参数地图，根据运些参数的位置相关的最优值，可W实现变体识别。
[0082] 在某些实施例中，参数地图可W使用在变体识别W外的其他流水线步骤，例如，比 对步骤。在比对过程中，用户可W选择部分参数，例如错配惩罚、间隙开口、不精确比对中的 间隙延伸W及比对候选者的阔值等。运些参数的最优值也可能是位置相关的，因为使用运 些参数的参数地图也可W改进比对的结果。
[0083] 在某些实施例中，位置相关的参数地图也可W被用在RNA测序流水线，例如，使用 基因组位置相关的参数地图来提高RNA测序流水线的准确性。
[0084] 单通多变量地图的计算流水线
[0085] 本发明的另一部分提供了一个在单通流水线中使用多个参数地图实现多个变体 识别的方法和系统。
[0086] 在某些实施例中，用户不清楚具体的样品特征，例如族群、或者它是否是癌症样品 等。在该情况下，用户可W使用多个参数地图进行样品识别。在某些实施例中，即使知道现 有基因组样品的种类，为了要学习变体在不同假设下如何变化，用户可能依然想要探索不 同的参数地图。一个探索多个参数地图的方法是对于每个参数地图，将测序数据多次输送 通过整个流水线，然而该方法会产生巨大的计算量。为了节省计算量并且节省时间，用户可 W只将测序数据通过流水线一次，同时多次运行使用了多个参数地图的变体识别模块。
[0087] 在运个单通多设定流水线中，不同变体识别共用的许多计算可W只运行一次，再 重复使用其结果，运样可W大大节约整体计算时间和探索大量可选参数地图的成本。例如， 如果用户想探索先验概率的多个参数地图，他们可W选择一列地图。在一个示例性实施例 中，流水线前期的一些不取决于先验概率参数的子步骤可W只运行一次，同时多次运行统 计模型一一对于每个参数地图运行一次。更进一步地，统计模型中的部分计算也不取决于 先验概率，所W也可W只运行一次，例如，片段过滤只取决于片段碱基质量和比对，和对于 特定碱基质量和基因型的碱基条件概率。总的来说，所有不受不同参数地图影响的共用计 算只需要被运行一次，然后被重复使用在多个变体识别中。通过移除冗余计算，在单通流水 线中，我们可W输出被识别的变体在使用各种参数地图的情况下的特性，与对每个参数地 图都运行整个流水线所需要的计算量相比，该计算量仅为其很小一部分。在单通多地图运 行时，被识别变体的特性(例如，它的概率)可W W图像形式显示为参数地图的选择的函数， 从而使得用户可W更好地从视觉上理解变体。运个单通多地图流水线方法可W被应用于多 种参数地图的探索，W通过重复使用不同参数地图中的共有中间结果来减少计算量。同样 地，运不仅仅适用于DNA组装流水线，也同样适用于RNA测序流水线，使得可W仅通过少量计 算实现单通多地图的探索。
[0088] 运个单通多地图流水线方法也可W被应用在向参考比对的计算中。当用户探索比 对参数地图时，对于可W完美向参考比对的片段，不需要对每个不同的比对参数地图都进 行重新比对，而可W仅运行一次但是将其结果应用到所有参数地图。只有不能比对的片段 W及和参考错配的片段需要对不同参数地图进行多次比对。
[00例逐步勘探计算流水线
[0090] 在某些实施例中，所述位置相关流水线可W被用在一个逐步勘探模式中。不同于 计算并且在一次通过中输出多个参数地图的结果，该流水线可W在第一次通过中的使用一 个参数地图的情况下计算输出结果，同时保存第一次通过产生的许多中间结果。如果用户 决定探索额外的参数地图，该流水线可W自动判别出那些不需要被再次计算的中间结果， 并将它们直接应用于一个新的通过。运样一来，使用新参数地图的额外通过只需要计算那 些会被新旧参数地图影响的流水线部分，从而节省了参数地图探索所需的计算量，哪怕并 不是一个单次通过。在某些实施例中，运个模式的优势在于用户可W在审阅了之前的参数 地图产生的结果后，对将进行的参数地图探索进行判断，从而可能可W减少一些对于不必 要参数地图的探索。
[0091] 在某些实施例中，运两个使用模式，即"单通多参数地图"和"逐步勘探"模式可W 更进一步地被结合在一起。在运种结合模式中，测序数据被输送通过才第一个单次通过中 使用多个参数地图的流水线，并且保存产生的中间结果。在审阅了第一个通过的结果之后， 用户可能决定在第二次流水线通过中探索额外的一个或多个参数地图。在第二次通过中， 流水线可W重复使用第一次通过的中间结果，从而节省新通过中的计算量。
[0092]虽然本发明已经被具体地展示，并依据具体实施方案给予了详细说明（其中部分 是优选的实施方案），但需要被本领域技术人员明白的是，在不脱离本文公开的发明的精神 和范围的情况下，依然可W对其做出各种形式和细节上的调整。
【主权项】
1. 一种使用计算流水线分析多个核酸序列片段中的变体的计算机辅助实现方法，其中 所述计算流水线包括取决于位置相关参数的模块，其特征在于，所述方法包括在处理器上 执行如下步骤： 接收多个核酸序列片段，其至少包含第一核酸序列片段和第二核酸序列片段； 将所述第一核酸序列片段和所述第二核酸序列片段分别比对到基因组中的第一位置 与第二位置，其中所述第一位置不同于所述第二位置； 基于基因组中的所述第一位置和所述第二位置，分别给所述位置相关参数赋予第一数 值和第二数值； 将所述第一核酸序列片段和所述第二核酸序列片段通过所述模块，分别使用所述第一 数值和所述第二数值的； 生成变体识别。2. 如权利要求1所述的方法，其特征在于，所述模块是变体识别模块。3. 如权利要求2所述的方法，其特征在于，所述位置相关参数是先验概率或变体识别阈 值。4. 如权利要求3所述的方法，其特征在于，所述先验概率是全基因组单核苷酸多态性概 率或插入缺失多态性概率。5. 如权利要求1所述的方法，其特征在于，所述模块是局域比对模块。6. 如权利要求5所述的方法，其特征在于，所述位置相关参数是碱基错配惩罚、间隙开 口、间隙延伸或比对候选者的阈值。7. 如权利要求1所述的方法，其特征在于，所述第一位置和/或第二位置是在外显子或 内含子内。8. 如权利要求1所述的方法，其特征在于，所述基因组来自于一个族群或区域群体。9. 如权利要求1所述的方法，其特征在于，所述基因组来自于健康对象或者患有一种疾 病的对象。10. 如权利要求9所述的方法，其特征在于，所述疾病是癌症。11. 一种非暂态计算机可读介质，其附带指令可以使其利用其所存储的计算流水线分 析来自多个核酸序列片段的变体，所述计算流水线包括使用位置相关参数的模块，所述指 令当被处理器执行时，会进行如下操作： 接收多个核酸序列片段，其至少包含第一核酸序列片段和第二核酸序列片段； 分别将所述第一核酸序列片段和所述第二核酸序列片段比对到基因组中的第一位置 和第二位置，其中所述第一位置不同于所述第二位置； 基于所述基因组中的所述第一位置和所述第二位置，分别给所述位置相关参数赋予第 一数值和第二数值； 将所述第一核酸序列片段和所述第二核酸序列片段通过所述模块，分别使用所述第一 数值和所述第二数值;并且 产生变体识别。12. 如权利要求11所述的方法，其特征在于，所述模块是变体识别模块。13. 如权利要求12所述的方法，其特征在于，所述位置相关参数是先验概率或变体识别 阈值。14. 如权利要求13所述的方法，其特征在于，所述先验概率是全基因组单核苷酸多态性 概率或插入缺失概率。15. 如权利要求11所述的方法，其特征在于，所述模块是局域比对模块。16. 如权利要求15所述的方法，其特征在于，所述位置相关参数是碱基错配惩罚、间隙 开口、间隙延伸或比对候选者阈值。17. 如权利要求11所述的方法，其特征在于，所述第一位置和/或所述第二位置是在一 个外显子或内含子内。18. 如权利要求11所述的方法，其特征在于，所述基因组来自一个族群或区域群体。19. 如权利要求11所述的方法，其特征在于，所述基因组来自健康对象或患有疾病的对 象。20. 如权利要求19所述的方法，其特征在于，所述疾病是癌症。
【文档编号】G06F19/22GK106021992SQ201610172078
【公开日】2016年10月12日
【申请日】2016年3月24日
【发明人】叶军, 周巍, 陈洛祁, 冯汉鹰, 陈洪, 刘晓峰
【申请人】知源生信公司（美国硅谷）