本文描述的方法通常涉及基因组分析，并且更具体地涉及下一代测序应用。
背景技术：
下一代测序
下一代测序(ngs)或大规模并行测序(mps)技术在过去的十年中显著降低了dna测序的成本。ngs在生物学中有广泛的应用，并极大地改变了研究或诊断方法学的方式。例如，rna表达特性谱分析或dna测序只能用少量基因和传统方法进行，如定量pcr或sanger测序。即使使用微阵列，在全基因组水平上对基因表达进行特性谱分析或鉴定突变只能在基因组尺寸相对较小的生物体上实现。利用ngs技术，rna特性谱分析或全基因组测序已成为生物学研究中的常规实践。另一方面，由于ngs的高通量，已经开发出多重方法，不仅要对更多区域进行测序，还要对更多样品进行测序。与传统的sanger测序技术相比，ngs能够并行检测不同基因中更多样品的突变。由于其优于传统测序方法，ngs测序仪正在取代sanger进行常规诊断。特别是，现在可以对个体的基因组变异进行常规分析，用于在精确医学实践中从遗传疾病诊断到药物基因组学精细调整的许多医学应用。ngs包括处理多个片段化的dna序列读取，通常是短序列(少于300个核苷酸碱基对)。然后可通过多种生物信息学方法将得到的读取结果与参考基因组进行比较，以识别小变体，例如对应于单核苷酸取代的单核苷酸多态性(snp)，以及与其参考相比，dna序列中核苷酸的短的插入和缺失(indel)。
靶向富集
在一些病理学中，特定的基因变体与疾病有关，例如某些形式的遗传性乳腺癌和卵巢癌中的brca1和brca2基因或囊性纤维化中的cftr基因。不是对来自个别样品的全基因组(wgs)进行测序，基因组分析可以集中于与疾病相关联的基因组区域，通过用区域特异性dna引物或探针组进行靶向以及例如使用pcr(聚合酶链式反应)的富集或扩增，生物dna样品专门用于对应于dna链上基因的亚区域。现在已经开发了许多下一代测序测定法及其原理作为即用型生物试剂盒，例如multiplicommastr或illumina
可以通过基于探针的杂交(在阵列上或在溶液中)或高度多重的基于pcr的靶向外显子富集从少量dna样品实现靶富集，以便基因覆盖/读取深度和扩增特异性(扩增正确的区域，如通过对期望的靶区域的进一步比对所测量的)都最大化。市售靶富集系统的实例包括安捷伦sureselecttm靶富集系统，罗氏nimblegenseqcapez，illuminanexterarapidcapture，安捷伦haloplextm和multiplicommastrtm。
为了最大限度地使用大规模并行处理ngs测序仪，大量样品在靶向的ngs实验中复用-例如，48个或更多靶富集样品的池可以因此同时输入到illuminamiseq测序仪。然后可以分析从ngs测序仪取出的原始测序数据以鉴定特定的子序列，例如通过与参考基因组比对。结果，对于患者样品中的给定扩增子，扩增可产生超过一千个读取结果。
cnv检测
在实践中，除了snp和indel之外，许多病理性遗传变异体是由dna序列中更显著的变化引起的。拷贝数变异(“拷贝数值”，“拷贝数异常”或cnv)量化样品dna序列中特定区域的拷贝的数量，与参考基因组相比这些特定区域可以进行长时间复制(高于正常值的拷贝的数量)或者可能超过数百个核苷酸的缺失(低于正常值的拷贝的数量)。虽然下一代测序方法在检测snp和indel方面显示比传统sanger测序更有效，但是在被靶向的ngs中检测cnv引发对于与参考基因组比对或与某些特定子序列匹配的许多特定挑战，因为读取长度通常低于300bp，即比整体cnv区域更短的序列。诸如mlpa(多重连接依赖性探针扩增)之类的现有技术cnv检测方法仍然需要单独的实验和基因组分析工作流。这限制了ngs在实际基因组分析应用中的优势，因为需要进行处理不同患者样品的不同工作流以检测病理重要性的cnv。而且，现有技术的cnv检测方法是低通量的，并且不能同时并行地检查对于大量样品和区域的cnv。因此文献中最近提出了许多解决方案，以更好地解决使用ngs工作流的cnv检测。一种方法，如例如wo2014151511中所述，在于比较靶扩增子的水平与对照扩增子的水平，以确定cnv的存在。然而，这种方法对照扩增子的选择非常敏感，这可能不容易得到。另一种方法在于进一步优化靶富集步骤，以便从靶富集样品池本身获得对于cnv检测的更好参考。例如，wo2015112619公开了使用假引物来将独特的参考核苷酸序列组分配给预先分选的扩增子的每个库，代价是额外的pcr扩增步骤和由cnv检测模块进行的迭代彻底检索。wo2014083147提出优化pcr引物，所述pcr引物具有在3'末端的与待分析序列互补的区域和在5'末端的非互补区域。后一种方法需要使用特定的测定试剂盒，这对于许多当前的应用来说太有限制了。
因此，需要一种更好的解决方案，以有效地检测cnv，可能同时针对大量样品和区域，在单个靶向的下一代测序实验中，无论其中采用哪种靶富集技术，以及以可能的自动化工作流从而优于现有技术方法促进研究和临床实验室实践。
概述
上述优点可以通过从使用靶富集技术富集的dna样品池检测拷贝数值(cnv)的方法来实现，每个富集的dna样品与来自扩增子/区域组的汇集片段的文库相关联，每个扩增子/区域用高通量测序仪进行测序以产生对于每个样品和对于每个扩增子/区域的覆盖计数，包括：使用数据处理单元标准化与每个样品相关联的覆盖计数；使用数据处理单元对于每个样品选择dna样品池内的参考样品组作为与所述样品标准化的覆盖计数具有最接近的标准化的覆盖计数的样品；并且对于每个样品，估计所述样品中的拷贝数值作为所述样品中的至少覆盖计数和对于所述样品的所选择参考样品组中的至少覆盖计数的函数。
参考样品的数量可以是样品的总数量的函数。它可以小于样品的总数量。
选择参考样品组可以包括计算在每个样品/丛内和在每个区域内都被标准化的覆盖计数之间的距离。距离可以是欧几里得距离。
标准化覆盖计数和/或选择参考样品组可以取决于在前一次迭代中对于每个样品和对于每个扩增子/区域的拷贝数值的估计。对拷贝数值的先前估计可以是预先定义的。对拷贝数值的先前估计可以迭代地计算，从对拷贝数值的预先定义的先前估计开始并且使用在每次迭代的cnv检测的结果作为后续迭代中对拷贝数值的先前估计，直到对拷贝数值的估计收敛，达到一个周期，或者迭代次数达到预定的限制。
对于每个样品和对于每个扩增子/区域，可以估计每个可能的拷贝数值的似然度。还可以使用隐马尔可夫模型来估计对于每个扩增子/区域的拷贝数值和它们的置信度水平。可以从结果中过滤掉置信度水平低于最低阈值的可能拷贝数值。
可以使用关于snp分数和覆盖计数的信息来计算对拷贝数值的估计。
可以将主成分过滤器应用于覆盖计数。
参考样品的数量可以取决于迭代指数。在一次迭代中，参考样品的数量nr可以等于样品的总数量n，而在另一次迭代中，参考样品的数量nr可以与样品的总数量n不同。
附图的简要说明
图1代表了靶向ngs基因组分析功能工作流。
图2示意性地显示了根据本公开的可能实施方案，作为对基因组数据分析仪cnv检测模块的输入的覆盖信息表的结构。
图3图示了示例性的a)基于扩增子的和b)基于探针的样品覆盖信息，所述样品覆盖信息用来自不同现有技术的靶向ngs平台的igv软件进行检索和可视化。
图4显示了根据本公开的可能实施方案的基因组数据分析仪cnv检测模块的流程图。
详细说明
基因组分析系统
图1显示了根据本公开的可能实施方案的示例性基因组分析系统，其包括靶向富集测定法100、下一代测序仪110和基因组数据分析仪120。
通过靶向富集测定法100处理dna样品池以产生通过基于扩增子的富集或基于探针的富集而制备的dna片段的文库作为对下一代测序仪110的输入，每组片段对应于不同的dna样品。片段的数量取决于应用。例如，在一些基于扩增子的实验中，靶富集可以需要150个引物来富集待从样品基因组中被靶向出来的75个不同的区域。在其他基于探针的实验中，探针富集可以选择例如来自413个所选择区域的dna片段。样品的数量也可以适应于下一代测序序列仪110的并行处理能力，例如可以通过illuminamiseq测序仪并行测序48个样品。可以使用其他ngs测序仪技术，例如罗氏454tmgsjunior或gsflx、illumina
下一代测序仪110分析输入样品并产生代表原始ngs测序数据的以计算机可读文件格式的序列读取结果。取决于ngs技术，ngs测序仪110可输出一个或多个文件。在一些实施方案中，fastq文件格式可以以用于正向和反向读取结构的两个不同文件或者作为单个联合文件使用。其他实施方案也是可能的。原始ngs测序数据被进一步输入到基因组数据分析器120。
基因组数据分析器120计算机系统(本文中也称为“系统”)120被编程或以其他方式配置为实现不同的基因组数据分析方法，诸如接收和/或组合测序数据和/或注释测序数据。
基因组数据分析器120可以是计算机系统或计算机系统的部分，包括中央处理单元(cpu，本文中的“处理器”或“计算机处理器”)，诸如ram的存储器和诸如硬盘的存储单元，以及通过通信网络例如互联网或本地网络与其他计算机系统通信的通信界面。在一些实施方案中，计算机系统可以包括一个或多个计算机服务器，其使得能够分布式计算，例如云计算，例如在基因组数据场中。在一些实施方案中，基因组数据分析器120可以被集成到大规模并行系统中。在一些实施方案中，基因组数据分析器120可以直接集成到下一代测序系统中。
如图1所示，基因组数据分析器120可以包括比对模块121，其将原始ngs测序数据与参考基因组进行比较。可以进一步分析比对结果(其可以用bam，sam或其他类似格式中的一个或几个文件表示，如生物信息学领域的技术人员已知的)，通过snp/indel检测模块124在对snp和indel多态性的搜索中进行。比对信息可被进一步过滤和分析以检索覆盖信息(或覆盖计数)。在本公开的一个实施方案中，覆盖提取模块122可以根据在整个基因组分析工作流中的上游应用的靶向富集100和ngs测序仪110技术来处理比对数据以提取覆盖信息。根据本公开的cnv检测模块123然后可以分析覆盖信息以识别和限定在原始dna样品中的拷贝数变异(cnv)。在一些实施方案中，可以通过在下一次迭代中使用来自前一步骤的cnv检测信息来迭代地操作cnv检测模块123。从图1中可以看出，具有单个靶富集步骤100和单个测序步骤110的相同ngs实验因此可以用于同时分析不同的snp/indel和cnv基因组变体，而不是如现有技术基因组数据分析工作流中运行单独的ngs/snp-indels检测和mlpa/cnv检测实验。
cnv检测-整体方案
图2示意性地显示覆盖之外的覆盖信息表，其中行表示扩增子(或者在基于探针的技术的情况下适当定义的区域)并且列表示样品。表中的符号'*'表示对于每个样品中每个扩增子/区域的覆盖信息(覆盖计数)。该覆盖计数被定义为与给定样品中的给定扩增子或区域对应的正向读取、反向读取和读取对(在双末端测序的情况下)的数量的合适函数。与扩增子的“对应性”的定义可以基于读取的开始和/或结束与扩增子的开始和/或结束之间的匹配，并且与区域的“对应性”可以基于读取与区域的重叠。然而，其他合适的定义也可以用于不同的实施方案中。也可以使用基于读取参数(例如读取长度或做图质量)的其他过滤器。在一个实施方案中，覆盖被定义为正向读取次数和反向读取次数的总和，而在另一个实施方案中，覆盖被定义为读取对的数量。在其他实施方案中，可以使用其他函数(例如，仅使用正向读取计数或仅使用反向读取计数，或者使用两个计数中的最大值)。此外，取决于靶富集测定法的性质，不同的函数可以用于不同的扩增子/区域。图3图示了使用比对的读取在igv软件中观察到的基于扩增子和基于探针的设计中，覆盖定义的可能选择。图3a)将扩增子覆盖显示为与对应区域比对的读取对的数量。图3b)将区域覆盖显示为与适当选择的靶区域重叠的读取总数(正向和反向)。
在cnv对应于在复制中的整个扩增子或区域的情况下，对应于该区域的ngs覆盖将异常地高。相反地，当cnv对应于在缺失中的整个扩增子或区域时，对应于该区域的ngs覆盖将异常地低。因此可以通过分析跨扩增子/区域和跨样品的覆盖信息表中的覆盖信息分布来检测cnv。覆盖计数ia,s可近似地分解成样品/丛依赖性和扩增子/区域依赖性贡献：
(eq.1)ia,s＝ra.sfapx,s+δia,s
其中：
·ia,s是样品s中对于扩增子(或区域)a的覆盖计数。
·fa是对于扩增子(或区域)a是特异性的扩增因子。
·px,s是表示对于丛管x中样品s的加工的dna材料量的因子(如果适用，例如在具有多个丛管的基于扩增子的技术的情况下，其中x指定对于扩增子a的丛)。
·δia,s是覆盖噪声，其被假设与总覆盖ia,s相比是小的。这种总体覆盖噪声可以来自实验室程序的各个阶段(dna提取和靶向富集)以及来自测序技术本身。在一个优选的简单模型中，覆盖噪声可以建模为δia,s＝ε(1)a,sia,s+ε(2)a,s√ia,s，其中第一项代表噪声的密集部分(与覆盖计数成比例)，第二项代表从有限数量的测序读取的随机波动引起的泊松噪声(该贡献与覆盖计数的平方根成比例)。
·ra.s是通过提出的cnv检测方法推导的拷贝数值(给定样品中给定扩增子/区域的多重性)。在大多数情况下(性染色体和同源区域除外)，ra.s＝2是正常的cnv值，并且偏离该正常值可能表明存在cnv(例如，ra.s＝1和ra.s＝3分别对应于杂合缺失和杂合复制)。请注意，在种系样品中，预计所有细胞都携带原始的个别基因组dna，因此拷贝数为整数。
因为即使在最好的实验室条件下，不同的样品可以具有略微不同的扩增因子fa，通过将每个样品与选择的“参考样品”组进行比较，cnv检测可以更可靠，所述“参考样品”被选为与正在进行测试的样品具有最佳相关联性的样品，而不是将所有样品一起进行比较。在多个丛的情况下，对于不同的丛，参考样品可以不同。
该方法假设一个批次中的样品数量足够大(在一个可能的实施方案中，样品数量至少为8，但其他选择也是可能的)，并且cnv足够稀少，从而对于任何扩增子/区域，大部分样品都有常规的拷贝数。
如图4所示，所提出的方法因此包括以下步骤：
·考虑到在前一次迭代(如果有)(400)中检测到的拷贝数值的估计值，通过样品s(以及丛x，如果适用)对覆盖信息进行标准化。在第一次迭代中，常规(通常)的拷贝数值可以用作假设，例如2。
·对于每个样品s(以下称为“当前样品”)：
ο选择参考样品组作为与当前样品s的标准化覆盖具有最接近的覆盖样式的标准化覆盖的样品(401)；
ο使用当前样品的标准化覆盖和来自其选择的参考样品的标准化覆盖，通过扩增子/区域a进行标准化(402)；
ο对于当前样品中的每个扩增子/区域a，
■估计参考水平和噪音(410)；
■对于每个扩增子/区域a，估计每个可能的拷贝数值的似然度(cnv状态)(420)；
■基于估计的拷贝数值的似然度，确定实际的拷贝数值(cnv)及其置信度水平(421)。
·通过考虑在先前迭代中已经检测到的cnv，从步骤(400)开始迭代，直到检测到的cnv稳定或达到一个周期，或达到最大迭代次数(430)；
·通过残留噪音和通过发现的cnv过滤样品。可选地，如有必要，从一开始重复整个过程，排除一些样品(440)。
现在将详细描述个别步骤如下。
样品方式的和丛方式的标准化
由于靶富集实验过程的具体情况，ngs工作流以外的原始覆盖信息未标准化。例如，在每个样品/丛中可能存在不同量的dna，导致沿着样品/丛轴的原始结果中的不同覆盖信息值。
在一个可能的实施方案中，为了从原始数据组中去除样品/丛偏差，可以采取样品/丛方式的标准化步骤400如下：将关于样品/丛的平均值确定为使用在先前迭代过程中计算的拷贝数值ra,s，标准化为单个拷贝的所有样品覆盖计数ia,s的平均值。对于所有样品/丛的所有扩增子的初始拷贝数值可以设定为正常值(通常ra,s＝2，性染色体和同源区域的情况除外)。然后可以分别对于每个样品/丛，将覆盖数据除以计算的平均值，将得到的样品/丛标准化的覆盖作为一级：
(eq.2)i(norm)a,s＝ia,s/平均(ia',s/ra',s|丛x中的所有扩增子/区域a')。
在一个可能的实施方案中，假定的纯合缺失ra,s＝0的区域可以从平均值计算中排除。在一个可能的实施方案中，平均值被计算为算术平均值。在另一个可能的实施方案中，平均值被计算为几何平均值。其他实施方案也是可能的。对于本领域技术人员来说显而易见的是，标准化方法的不同实施方案也可以用于cnv检测方法的不同迭代中。
参考样品的自动选择
所提出的方法允许从来自靶ngs实验数据的标准化覆盖率信息中参考样品的自动选择，而不需要用户输入来提供或手动选择专用对照样品，如现在将更详细描述的。
通常，对于每个样品/丛，可以通过cnv检测模块123自动选择合适的参考样品(步骤401)，作为与当前样品s0具有最接近的覆盖样式的那些。在一个可能的实施方案中，可以通过对于每个样品计算与当前样品s0的距离，然后按照距离增加的顺序对样品进行排序，并且从列表的顶部选择一定数量的样品(具有最小距离)，来选择最接近的覆盖样式。
对于统计学领域的技术人员来说显而易见的是，有许多可能的方式来定义和计算样品之间的距离。在一个可能的实施方案中，我们首先使用样品/丛标准化的计数i(norm)a,s来计算样品/丛/扩增子/拷贝数标准化的计数作为向量va,s：
(eq.3)va,s＝(i(norm)a,s/ra,s)/中值(i(norm)a,s'/ra,s'|全部样品s')。
在一个可能的实施方案中，具有假定的纯合缺失ra,s＝0的区域可以从中值计算中排除。
在一个可能的实施方案中，任何样品s与当前样品s0之间的距离可以被定义为矢量va,s和va,s0之间的欧几里得距离。在一个可能的实施方案中，可以计算矢量va,s和va,s0之间的相关联性。
在具有多个丛管的基于扩增子的技术的情况下，优选地对于每个丛x分别计算距离，可能导致对于不同丛的不同参考样品组。例如，在两个丛(例如，对于cftrmastrmultiplicom靶向富集试剂盒)的情况下，对于每个样品可以有两个参考组：每个丛一个。
定义和计算样品之间距离的其他方法也是可能的。例如，可以使用算术或几何平均值代替中值或反之，或者可以使用其他类型的度量来代替欧几里得度量。根据不同的标准，算法还可以排除某些区域或将不同的权重归于不同的区域。在一些实施方案中，聚类算法可以用于分配距离。其他类型的算法也是可能的。
在计算每个样品与当前样品s0之间的距离之后，可以选择参考样品作为具有最小距离的特定数量的样品。如本领域技术人员将理解的那样，应仔细选择参考样品的数量。这个数量应该足够大，以获得良好的统计相相关性；特别是在每个扩增子位置，参考样品组应具有大部分正常拷贝数(无突变)。另一方面，这个数量应该足够小，以便只比较运行中的相似样品并过滤掉异常值。
在一个可能的实施方案中，可以作为样品的总数量n的函数来选择参考样品数量nr。在一个可能的实施方案中，该函数可以被定义为：
(eq.4)nr＝[αn+β√n+γ]
具有适当选择的系数α、β和γ，并且[..]表示整数部分。在一个可能的实施方案中，参数α、β和γ可以选择为使得nr＝[0.25*n]+2，以选择大约25％的样品作为参考样品。系数以及更一般地函数的其他选择也是可能的。在另一个实施方案中，nr可以不仅取决于样品的总数量n，而且取决于数据的其他属性，例如取决于覆盖计数的波动水平。此外，在一些实施方式中，对于不同的当前样品s0(或对于不同的丛)，参考样品的数量可以是不同的，例如，如果参考样品被选择为那些低于与当前样品s0的特定截止距离的距离。基于计算出的距离选择参考样品的算法的其他选择也是可能的。
虽然通常保持参考样品的数量nr小于样品的总数量n可能是有益的，为了从参考中排除噪声样品，在可能的实施方案中，参考样品的数量nr可以等于样品的总数量n。
在一些实施方式中，计算距离和/或选择参考样品的具体实施方案的选择可以在迭代之间变化。例如，在一个可能的实施方案中，在cnv检测的迭代之一中(在一个可能的实施方案中，具体地在第二次迭代中)，参考样品的数量nr被认为等于n：结果，对于这个一次迭代所有样品都用作参考样品。对于生物信息学领域的技术人员来说显而易见的是，这种灵活的方法可以在频繁cnv的情况下提高检测性能。
扩增子方式的标准化
除了丛/样品偏差之外，不同扩增子/区域之间还可能存在覆盖信息离散，因为扩增效率倾向于取决于区域，因此也导致沿着扩增子/区域轴的不同覆盖信息值。因此，例如通过将覆盖信息除以每个扩增子a的中值，再次使用标准化为来自先前迭代的拷贝数的假设正常值(或第一次迭代中的正常拷贝数)的覆盖水平，通过cnv检测模块123在步骤402中对丛标准化的数据进行另外的标准化的扩增子方式。在该步骤中，标准化可以特别地在包括步骤401中选择的参考样品的缩减样品组内执行：
(eq.5)ca,s＝i(norm)a,s/中值(i(norm)a,s'/ra,s'|参考采样s')
在一个实施方案中，除了所有参考样品之外，当前样品s0可以被包括在中值计算中。在另一个可能的实施方案中，当前样品s0可以从中值计算中排除。在一个可能的实施方案中，具有假定的纯合缺失ra,s＝0的区域可以从中值计算中排除。
估算对于每个扩增子或区域的参考水平和噪音
在所提出的方法的这个阶段，已经为每个样品s计算了以下数据：
(1)参考样品组(如果适用，对于每个丛x)；
(2)样品/丛/扩增子标准化的覆盖水平ca,s。
估计参考覆盖水平和噪声(不确定性)的另一个步骤410可以进一步应用如下。
首先，通过使用在先前迭代中计算的假设的拷贝数值ra,s，可以将覆盖水平ca,s转换成每个拷贝的标准化的覆盖水平：
(eq.6)c(0)a,s＝ca,s/ra,s
在第一次迭代中，可以假设对于所有样品和对于所有扩增子/区域的拷贝数是正常的，例如ra,s＝2。
其次，可以估计对于每个扩增子/区域a的参考标准化的覆盖水平ca。在一个可能的实施方案中，可以假设标准化的覆盖水平ca等于1。在另一个可能的实施方案中，可以将覆盖水平ca计算为标准化的值c(0)a,s在参考样品上的平均值，其中异常值(偏离平均值的值大于某个阈值，例如三个标准差)被除去。正如统计学领域的技术人员将会理解的那样，估计参考覆盖水平的其他选择也是可能的。
第三，可以估计对于每个扩增子/区域的噪音水平。噪声水平可以被定义为覆盖的预期相对(除以平均值)均方根不确定度。在一个可能的实施方案中，噪声水平σa,s可以被估计为
eq.(7)σa,s＝max(σs,σa,1/√ia,s)，
其中：
·σs和σa是对于给定样品和对于扩增子(在包含给定样品和参考样品的组内)的c(0)a,s的标准偏差。
·1/√ia,s是对应于原始覆盖值ia,s的泊松噪声的相对均方根偏差。
通过使用跨样品和跨扩增子的标准偏差，eq.7考虑了噪音样品和噪音扩增子两者的似然度。为了计算σa，排除异常值可以是有益的。在一个可能的实施方式中，3σ区间之外的数据点可以被排除。对于统计学领域的技术人员来说显而易见的是，也可以使用估算噪声水平σa,s的其他方法来代替eq.7。
在一些实施方案中，估计对于每个扩增子/区域a的参考标准化的覆盖水平ca的步骤可以在第一次迭代之一中(例如，在一个可能的实施方案中，具体地在第一次迭代中)被修改，以便检测跨所有扩增子/区域的大部分的cnv。在这种情况下，通过样品的简单标准化可能无法估计正确的参考水平。这个问题可以通过额外的算法来解决，该算法根据“对照”扩增子/区域的特定组(假设它们大部分是无cnv的)或基于标准化的水平的最佳匹配为整数来确定参考水平。
计算不同拷贝数的似然度
在此步骤中，可以将对于每个扩增子/区域a和对于每个样品s的标准化的覆盖水平和噪声水平进一步转换为对数似然度la(定义为在假设的给定拷贝数并假设给定噪声水平，特定覆盖水平的似然度的负对数)(步骤420)。在一个可能的实施方案中，高斯噪声模型可以用于计算对数似然度la：
eq.(8)la(r)＝min((ca,s/(rca)-1)2/(2σa,s2),lmax)，
其中ca,s、ca和σa,s分别是当前样品s和扩增子/区域a的覆盖、参考标准化的覆盖水平和噪声水平。因此可以针对范围从0到某个最大值的r的所有整数值计算对数似然度la(r)(在一个实施例中，我们选择r的最大值等于6，但是其他实施方案也是可能的)。在r＝0的情况下，eq.8可被la(0)＝min((ca,s/ca)2/(2σ02),lmax)代替，其中σ0是完全缺失的假设噪声水平(在一个实施方案中，我们选择σ0＝0.01，但其他选择也是可能的)。为了考虑到大波动不服从正态分布，对数似然度可能被限制在特定的值lmax。
在一些实施方案中，lmax可以取决于扩增子/区域a。在其他实施方案中，可以使用其他噪声模型来代替eq.8。
寻找cnv及其置信度
使用在前一步中计算的对数似然度，可以找到最可能的cnv状态，并计算它们的置信度。在一个可能的实施方式中，隐马尔可夫模型(hmm)方法可以用于该目的(步骤421)。例如可以在例如s.ivakhno等人的“cnaseg-anovelframeworkforidentificationofcopynumberchangesincancerfromsecondgenerationsequencingdata”，bioinformatics(2010)26(24))：3051-3058中找到在cnv检测中使用hmm的实例。其他实施方案也是可能的，例如la(r)与适当选择的阈值的简单比较，类似于mlpa推荐的程序，如本领域技术人员已知的。
在一个可能的实施方案中，hmm方法可以通过将hmm得分如下定义而实现：
eq.(9)shmm({ra})＝∑a(la(ra)+pnb(ra)+psw(ra,ra+1)，
其中hmm得分shmm({ra})是当前样品中对于每个扩增子/区域的假设拷贝数ra的组的函数，la(ra)是在前一步骤计算的对数似然度，并且pnb(ra)和psw(ra,ra+1)是与非正常拷贝数相关联的附加罚分并且具有在相邻扩增子/区域(表示为a和a+1)之间转换不同拷贝数。可以选择参数pnb(ra)和psw(ra,ra+1)以提供良好的性能并反映样品中具有cnv的贝叶斯先前预期。在一个可能的实施方案中，函数pnb(ra)和psw(ra,ra+1)可以被选择为独立于区域/扩增子a，但是在其他实施方案中它们本身可以取决于区域/扩增子a。例如，在其他可能的实施方案中，函数pnb(ra)和psw(ra,ra+1)可以是区域的长度，区域之间的间隙，扩增子之间可能的重叠或给定区域中cnv的已知收集统计数据的函数。其他实施方案也是可能的。
一旦hmm得分由eq.9所定义，正如本领域技术人员已知的那样，可以使用前向-后向算法来找到使hmm得分和“置信度”值的组最小化的cnv状态的组{ra}。可以将位置a处的置信度值定义为在状态ra不同于其最佳值时hmm得分的最小可能增加。这个置信度的统计意义是在确定ra中的错误概率的负对数。
在实践中，可以引入两种置信度：
·确定拷贝的确切数量的“数字”置信度；
·将状态分类为普通/插入/缺失而不指定拷贝数的确切值(在插入或缺失的情况下)的“变体”置信度。因此，“变体”置信度总是高于或等于“数字”置信度。
所有参数可以进一步优化以获得更好的性能和/或可以以各种方式修改hmm模型，这对于本领域技术人员是显而易见的。其他hmm模型也可以应用。例如，根据所有可能的cnv的统计数据，对插入和缺失使用不同的罚分，并且引入切换罚分对cnv状态ra和ra+1之间差异的依赖性是值得的。
主迭代
从上述步骤可以看出，噪声和参考值计算取决于假定的拷贝数值(所谓的“对拷贝数值的先前估计”)，其中正常值例如2是起点。为了自身一致地确定cnv状态，可以使用在每次迭代的cnv检测结果作为后续迭代中对拷贝数值的先前估计来对算法进行数次迭代(步骤430)。实际的实验表明，在高质量数据的情况下，少数迭代足以有效检测真实的cnv值。有时(在数据噪音的情况下)，迭代可以进入周期性循环。在那种情况下，一旦检测到循环模式，该算法就可以停止。
一些靶向富集技术可以包括所谓的“对照扩增子/区域”：在目的区域外部的扩增子或区域，通常跨基因组而广泛分布。这种对照扩增子可以用于标准化覆盖信息。在一个实施方案中，具体在第一次迭代中，在样品/丛标准化(步骤400)中仅使用来自对照扩增子的覆盖信息。这可以允许更稳健地检测可能的大cnv(例如，整个基因的缺失)，如在这种情况下，在第一次迭代时，对于大cnv区域的拷贝数值ra,s将被设置为正确的值，并且将在随后的迭代中保持这个值。否则，在计算噪声和参考值时，对照扩增子可以与测试扩增子处于同一基础上，但它们不需要包含在算法的hmm部分中。
最终过滤
在最后一次迭代之后，所提出的方法提供cnv水平的结果值以及它们的置信度水平。在一个可能的实施方案中，可以为置信度值设置最小阈值，在该阈值以下，可以假设个别扩增子的结果是“不可靠的”并且可以从结果中过滤掉(步骤440)。在其他实施方案中，某些样品也可以基于剩余样品噪声σs或者其中一个丛中的残余噪声或者检测到的不切实际的大量cnv被排除为“不可靠”。将样品标记为“不可靠的”确切条件可以取决于靶富集和测序技术的细节。
在步骤440中从最终结果中丢弃“不可靠”样品的cnv结果。最后，如果过多样品(在一个实施方案中超过所有样品的一半以上)被过滤掉为“不可靠的”，则整个过程需要从开始进行重复并且排除“不可靠的”样品。该选项可为运行提供更好的性能，其中大部分样品在总体分析工作流中的初始靶富集或测序步骤中存在技术问题。
优化(改进)
在一些实施方案中，可以应用算法的进一步改进。一种可能的进一步实施方案可以应用主成分过滤，类似于m.fromer等人，“discoveryandstatisticalgenotypingofcopy-numbervariationfromwhole-exomesequencingdepth”，am.j.hum.genet(2012)91：597-607。例如在主cnv检测算法步骤400至430之前，或在初步标准化的数据集上，例如在通过样品/丛标准化(步骤401)之后或在通过扩增子/区域的标准化(步骤402)之后，可以将主成分过滤器应用于原始数据集。在一个实施方案中，过滤器在特别选择的训练数据集上被训练一次，随后在没有进一步更新的情况下使用。在其他实施方案中，过滤器可以包括从新数据集中学习。
在又一个实施方案中，如图1所示，来自并行snp/indel检测模块124的结果也可以用作对cnv检测模块123的输入，以进一步加强cnv检测。在这种情况下，对于杂合snp的覆盖分数的信息可以用来偏倚对cnv值的决定。例如，33％的snp分数可以是支持复制的有力论据(拷贝数等于3)。在一个可能的实施方案中，可以通过对发现一个或多个杂合snp的区域添加对pnb(ra)的适当选择的贡献，在hmm步骤421引入该偏差。
在另一个可能的实施方案中，所提出的cnv检测方法可适用于同源(相同或几乎相同)区域或假基因的情况。在这种情况下，正常拷贝数值可以不同于2(例如，在一对同源区域的情况下，正常拷贝数值等于4)。通过调整正常拷贝数的假设值(例如从2到4)并且通过使用在主cnv算法步骤400-430中与所考虑的区域同源所有区域中的读取总数，cnv检测算法可以被推广应用于这种情况。另外，同源区域之间的覆盖差异可以以与前述实施方案描述中的杂合snp类似的方式使用，以偏置参数pnb(ra)。
另一种正常拷贝数可以不同于2的情况是性染色体(x和y染色体)的情况。在一个可能的实施方案中，根据患者的性别调整x染色体和y染色体中的区域的正常拷贝数。在另一个可能的实施方案中，可以通过比较x染色体、y染色体和常染色体之间的覆盖信息，根据它们在靶扩增技术中的存在，自动地确定患者的性别。
实验结果
如图4流程图所描绘的用于检测cnv变体的所提出的方法的功效已经在实验中与mlpa方法进行了比较，该实验包括源自一个实验室的11批次中的474个样品，并使用brcatrusseq技术和miseq下一代测序流水线。brcatruseq测定法的一个已知特征是它具有许多短的扩增子，但是覆盖噪音相对高，并且没有对照扩增子。根据所提出的迭代方法，仍然可以通过比较基因之间的覆盖水平来检测整个基因(brca1或brca2)的缺失。
这些样品包含由mlpa确认的16个cnv。使用我们的cnv模块独立地分析相同的样品。然后将结果与mlpa确认的变体进行比较。作为通过算法与mlpa方法相比捕获的cnv百分比而测量的检测灵敏度测量为100％，即成功检测到全部16个cnv。此外，采用相同或替代方法(越低越好)建议对于重新测试cnv的拒绝或假阳性样品的百分比测量为4.2％，而最大值为10％可能在最好的实验室实践中可接受的。因此，所提出的基因组数据分析方法使得能够达到与现有技术mlpa方法相似的cnv检测灵敏度和准确度，同时能够使用单个靶ngs实验流水线，其在研究或临床实践中带来显著的实际优势。
其他实施例和应用
尽管上面的详细描述包含许多具体细节，但是这些不应该被解释为限制实施方案的范围，而是仅仅提供几个实施方案中的一些实施方案的图示。
虽然以上已经描述了各种实施方案，但应该理解的是，它们已经以示例而非限制的方式呈现。对于相关联领域的技术人员显而易见的是，在不脱离本发明的精神和范围的情况下，可以在其中进行形式和细节上的各种改变。实际上，在阅读了以上描述之后，相关联领域的技术人员将明白如何实现替代实施方案。
另外，应该理解的是，任何突出显示功能性和优点的图仅出于示例目的而被呈现。所公开的方法足够灵活且是可配置的，使得它们可以以不同于所示的方式来使用。
虽然在说明书，权利要求和附图中经常使用术语“至少一个”，但是术语“一个”，“一种”，“该”，“所述”等也在说明书、权利要求和附图中表示“至少一个”或“至少一种”。
最后，申请人的意图是，只有包含明确表述“用于......的装置”或“用于......的步骤”的权利要求才能被解释为属于35u.s.c.112，第6段。没有明确包括短语“用于......的装置”或“用于......的步骤”的权利要求不应根据35u.s.c.112，第6段进行解释。