用于检测异常核型的方法和系统
1.本技术是申请日为2017年2月13日，申请号为201780021833.8(国际申请号pct/us2017/017734)，发明名称为“用于检测异常核型的方法和系统”的发明专利申请的分案申请。
背景技术：2.对人类基因组样本的准确医学解释需要了解潜在核型。用于鉴定异常核型，例如拷贝数变体(cnv)的方法，包括在比较基因组杂交(cgh)中使用dna微阵列，例如使用荧光原位杂交(fish)、克隆和pcr产物测定、寡核苷酸阵列、基因分型阵列(carter np,nature genetics 2007；39s16-21))。然而，阵列技术的缺点是可能难以定义(调用)推定的cnv。
3.用于由下一代测序数据检测染色体异常的方法很少。已经使用了某些下一代测序全基因组拷贝数变体方法，例如基于读对、分解读段、读段深度和组装的方法(pirooznia等，front.genet.2015；6；138)。然而，现有应用集中于分析来自母体血浆样本的非常轻的脱脂全基因组测序(wgs)数据，以检测非整倍体胎儿游离dna的部分以进行非侵入性产前检测(nipt)。在癌症基因组学中已经在一定程度上探索了新一代测序，但考虑到准确测量体细胞染色体异常中克隆镶嵌(mosaicism)程度所必需的覆盖深度，这些分析通常基于snp阵列。
4.已开发的现有方法不能用于由群体规模的全外显子组测序(wes)数据检测异常核型。这些和其他缺点在本公开中得以解决。
技术实现要素：5.应理解，下面的概述和下面的详述都仅为示例性和说明性而非限制性的。公开了用于检测异常核型的方法和系统。示例方法可以包括测定多个样本中每条染色体的读段覆盖度数据、杂合snp的等位基因平衡分布和未观测到杂合性的染色体区段，其中每条染色体包含多个基因组区域；测定所述多个样本中每条染色体的期望读段覆盖度数据；测定所述多个样本中至少一条染色体的读段覆盖度数据与期望读段覆盖度数据之间的偏差；测定对于多个样本中至少一条染色体的多个双等位基因snp而言，等位基因平衡分布与1:1的期望比的偏差；确定所述偏差是出现在整条染色体上还是仅出现在鉴定的染色体的一部分上；使用补体中的读段覆盖度和等位基因平衡数据进一步细化和验证多个样本中至少一条染色体的鉴定偏差，并将所述至少一条染色体鉴定为异常核型。
6.其他优点将在下面的描述中进行部分阐述或者可以通过实践来了解。所述优点将借助于所附权利要求中特别指出的要素和组合来实现和获得。
附图说明
7.并入本说明书中并构成其一部分的附图说明了实施方案，并与说明书一起用于解释所述方法和系统的原理：
8.图1是说明异常核型检测示例方法的流程图；
9.图2是说明示例线性回归模型的图表；
10.图3是说明表现出大残差的异常核型的图表；
11.图4是说明异常核型检测示例方法的另一流程图；
12.图5显示了说明gc含量和覆盖度的关系的图表；
13.图6是说明鉴定的异常核型和异常值的图表；
14.图7a、7b、7c、7d、7e和7f是等位基因平衡图，其显示了样本的9、13和20号染色体上的异常。子图编号为染色体编号。阴影条(701)表示杂合snp等位基因平衡为0.5的期望正常变异范围。实线(702)表示全染色体中值等位基因平衡。虚线(703)表示约20snp滚动窗口中的局部中值等位基因平衡。线(704)表示连续性纯合片段(runs-of-homozygosity)；
15.图8是与图7a-f中相同的样本的读段覆盖度图；
16.图9a、9b、9c、9d、9e和9f是等位基因平衡图，其显示了样本中21号染色体上的异常和涵盖整个x染色体的连续性纯合片段，暗示核型正常的雄性样本仅有一条x染色体。阴影条(901)表示杂合snp等位基因平衡为0.5的期望正常变异范围。实线(902)表示全染色体中值等位基因平衡。虚线(903)表示约20snp滚动窗口中的局部中值等位基因平衡。线(904)表示连续性纯合片段；
17.图10是与图9a-f中相同的样本的读段覆盖度图；
18.图11是说明异常核型检测示例方法的流程图；
19.图12是样本4号染色体上的示例等位基因平衡图，其中检测到大的连续性纯合片段(1202)，其由于异常区域中纯合snp间有少量非零等位基因平衡而具有重叠localhetab事件(1204)；
20.图13a是所有样本的x染色体相对于y染色体覆盖度比率的图和实线1306所示用于测定雄性(1302)和雌性(1304)样本的阈值。另外，可以使用y染色体覆盖度比率阈值(虚线1308)鉴定有y染色体重复的雄性样本；
21.图13b是21号染色体的示例图，证明期望的全染色体中值杂合snp等位基因平衡(chromhetab)相对于在特定读段深度阈值下或高于特定读段深度阈值(例如，50x覆盖度，“pcttargetbases50x”qc度量)的覆盖碱基分率增加；可以基于覆盖度度量，基于观测到的相对于期望的chromhetab的偏差显著性来分配“层级”评级；
22.图14是所有雄性样本x染色体上的chromhetab值(x轴)相对于chromhetab值的计算中所包括的snp数量(推定的杂合snp；y轴)的图。线表示用于基于由大量snp所支持的高的非零chromhetab值来区分x染色体上具有重复的雄性样本的阈值；
23.图15是所有localhetab事件(黑色和灰色点)(面积大于示例阈值即垂直线)相对于事件中包括的杂合snp的数量(y轴)的图，其中对角线表示示例层级评级阈值。灰色点表示具有重叠roh事件的事件；以及
24.图16是说明用于执行所公开的方法的示例性操作环境的方框图。
具体实施方式
25.在公开和描述本方法和系统之前，应理解所述方法和系统不限于特定方法、特定部件或特定实施方式。还应理解，本文使用的术语仅仅是为了描述特定实施方案的目的，而非旨在为限制性。
26.如说明书和所附权利要求书中所用，除非上下文另外明确指出，否则单数形式“一”、“一种(个)”和“所述(该)”包括复数指示物。本文可以将范围表述为“约”一个特定值，和/或至“约”另一个特定值。在表述此类范围时，另一个实施方案包括从所述一个特定值和/或至另一个特定值。类似地，通过使用先行词“约”将数值表述为近似值时，应该理解该特定值形成了另一个实施方案。应该进一步理解的是，每个范围的端点对于另一个端点很重要并且独立于另一个端点。
[0027]“任选的”或“任选地”意指随后描述的事件或情形可能发生或可能不发生，并且该描述包括所述事件或情形发生的情况和不发生的情况。
[0028]
在本说明书的描述和权利要求书全篇中，词语“包含”和该词的变型如“含有”意指“包括但不限于”，并非旨在排除例如其他部件、整数或步骤。“示例性”意指“一个实例”，并非旨在表达优选或理想实施方案。“例如”不是在限制性意义上使用，而是用于解释目的。
[0029]
应理解，所公开的方法和组合物不限于所描述的特定方法、方案和试剂，因为它们可以变化。还应理解，本文所用的术语仅仅是为了描述特定实施方案的目的，而非旨在限制本方法和系统的范围，本方法和系统的范围将仅受所附权利要求限定。
[0030]
除非另有定义，否则本文使用的所有技术和科学术语具有与所公开的方法和组合物所属领域的技术人员通常所理解的相同含义。尽管与本文描述的那些方法和材料类似或等同的任何方法和材料可用于本发明的实践或试验，但特别有用的方法、设备和材料如文所述。本文引用的出版物及其引用的材料特此通过引用明确并入。本文的任何内容均不应解释为承认本发明无权凭借在先发明而先于此类公开。不承认任何参考文献构成现有技术。对参考文献的讨论陈述了作者所声明的内容，并且申请人保留质疑所引用文件的准确性和切合性的权利。将清楚地理解，尽管本文提及了许多出版物，但这样提及并不构成承认任何这些文献形成本领域公知常识的一部分。
[0031]
公开了可用于实现所公开的方法和系统的部件。本文公开了这些和其他部件，并且应当理解，当公开这些部件的组合、子集、交互作用、分组等时，虽然可能未明确公开对这些中每种不同的单个和集体组合排列的特定提及，但对于所有方法和系统，每一种在本文中都有特别考虑和描述。这适用于本技术的所有方面，包括但不限于所公开方法中的步骤。因此，如果存在可以执行的多个附加步骤，则应理解这些附加步骤中的每一个都可以利用公开方法的任何特定实施方案或实施方案的组合来执行。
[0032]
通过参考以下对优选实施方案和其中包括的实例的详细描述以及附图及其之前和之后的描述，可以更容易地理解本方法和系统。
[0033]
如本领域技术人员将理解的，所述方法和系统可以采用全硬件实施方案、全软件实施方案或结合软件和硬件方面的实施方案的形式。此外，所述方法和系统可以采用计算机可读存储介质上的计算机程序产品的形式，其具有包含在存储介质中的计算机可读程序指令(例如，计算机软件)。更具体地，本方法和系统可以采用web实现的计算机软件的形式。可以利用任何合适的计算机可读存储介质，包括硬盘、cd-rom、光存储设备或磁存储设备。
[0034]
下面参考方法、系统、装置和计算机程序产品的方框图和流程图图解来描述所述方法和系统的实施方案。应当理解，方框图和流程图图解的每个方框以及方框图和流程图图解中方框的组合分别可以通过计算机程序指令实现。这些计算机程序指令可以加载到通用计算机、专用计算机或其他可编程数据处理装置上以产生机器，使得在计算机或其他可
编程数据处理装置上执行的指令产生用于实现流程图一个或多个方框中指定的功能的方式。
[0035]
这些计算机程序指令也可以存储在计算机可读存储器中，该计算机可读存储器可以指示计算机或其他可编程数据处理装置以特定方式起作用，使得存储在计算机可读存储器中的指令产生包括用于实现流程图一个或多个方框中指定的功能的计算机可读指令的制品。计算机程序指令也可以加载到计算机或其他可编程数据处理装置上，使一系列操作步骤在计算机或其他可编程装置上执行，以产生计算机实现的过程，使得在计算机或其他可编程装置上执行的指令提供用于实现流程图一个或多个方框中指定的功能的步骤。
[0036]
因此，方框图和流程图图解的方框支持用于执行指定功能的方式组合、用于执行指定功能的步骤组合和用于执行指定功能的程序指令装置。还应理解，方框图和流程图图解的每个方框以及方框图和流程图图解中方框的组合可以由执行指定功能或步骤的基于专用硬件的计算机系统、或专用硬件和计算机指令的组合来实现。
[0037]
一方面，公开了用于从群体规模的全外显子组测序数据(也称为karyoscan)中检测具有异常核型的样本的方法。可以通过染色体上的读段深度分布来检测异常核型，但是多种因素会影响区分真正染色体异常与噪声的能力。pcr扩增受因gc含量和实验条件而有偏倚，从而往往导致基因组上dna片段的不均匀扩增。另外，外显子组捕获技术不能产生均匀的靶标覆盖度。因此，任何特定染色体或染色体区域的期望覆盖度取决于多种因素，其中一些因素可测量，而另一些因素则不可测量。
[0038]
所公开的方法即图1说明的示例方法100可以在102处计算单个样本相对于每条染色体的读段覆盖度曲线。为了减小读段覆盖度的偏倚，可以在104处测定外显子区的代表性gc含量和可映射性度量，因为在gc含量接近50％且可映射性高的区域中变异最小。可以按外显子组区域上gc含量在一定范围(例如45-55％)内且可映射性高于阈值的读段深度的总和，为每条染色体i测定稳健读段覆盖度曲线ri。与中值染色体标签密度不同，这种度量允许分辨亚染色体。
[0039]
然后可以在106处对染色体读段覆盖度曲线进行归一化，以表示每条染色体相对于其他常染色体的全外显子组读段覆盖度比率。染色体i的全外显子组覆盖度比率γi可以表示为：
[0040][0041]
其中(a-i)是除染色体i以外的常染色体集合，并且γi针对所有常染色体和x染色体而测定(可以独立考虑y染色体)。因此，染色体i的覆盖度比率是染色体i与所有其他常染色体相比的读段比率。
[0042]
染色体异常表现为γi与期望值有偏差。然而，即使在正常(二倍体)核型样本间，γi的期望值也不是恒定的，并取决于实验条件。可以在108处使用线性回归模型预测每个个体每条染色体上γi的期望值图2示出了22号染色体进行线性回归拟合后的观测值(γi)和期望值的实例。来自picard的与读段深度变异相关的测序质量控制(qc)度量可以用作该模型中的协变量。qc度量可以包括例如gcdropout、atdropout、meaninsertsize、
onbaitvsselected、pctpfuqreads、pcttargetbases10x、pcttargetbases50x和/或诸如此类中的一个或多个。
[0043]
虽然这些qc度量可以描述在读段覆盖度上观测到的很大一部分变异，但是在使用先前已知的方法获得的结果中可以反映不可测量的附加偏倚。这些偏倚在具有相似外显子组gc含量分布的染色体之间相关，并且包括相似染色体的γi值作为附加协变量可以将方差降低至可接受的水平。一方面，虽然这对于模型特异性而言是有益的，但有一个缺点是这些其他染色体本身可能核型异常，这可能在靶染色体上导致假阳性调用。本发明方法提供的优点是通过限制来自于其他染色体的协变量数量将靶染色体上的假阳性调用减到最少。例如，来自于其他染色体的协变量数量可以限制为两个。
[0044]
因此，可以为n个样本的整个集合的每条染色体回归线性模型，其中：
[0045]
其中染色体j、k定义为相对于染色体i的gc含量分布具有最小d统计值的两条常染色体。在一些方面，性别(由y染色体覆盖度阈值定义)可以用作x染色体的附加协变量。
[0046]
异常核型的检测可以基于110处检测由残差定义的、γi与特定样本期望值的偏差。然而，落在qc度量空间极值上的样本估计值固有地可以产生具有更高方差的平均估计值，使得不能假设在所有样本中对原始残差的解释是一致的。在112处，对于具有协变量x的单个样本而言，公开的方法可以相对于平均估计值的标准误差对残差进行z分数归一化(参见图6)：
[0047][0048]
其中se为残差标准误，n是用于拟合模型的样本数量，以及：
[0049][0050]
可以在114处测定每条染色体基于z分数的p值以鉴定显著较大的残差，其表示染色体i的异常核型。一方面，p《0.05且q《0.05(fdr调整的p)的p值截止值可用于鉴定显著较大的残差。参见图3，其中示出了线性回归拟合后的观测值(γi)和期望值另一方面，可以使用高达0.1的p值。
[0051]
大的残差可以是目标染色体真正异常核型以及异常协变量值的结果(由于qc度量空间中的异常值或者协变量染色体之一的异常核型)。在116处，由于不寻常的协变量，通过在每条染色体的线性模型上标记具有极端杠杆值(常常表示为hi，其中1/n《hi《1)的样本来检测异常值。杠杆值量化了样本的x值(协变量)对模型的影响程度。杠杆值可用于标记不表示目标染色体上的真正异常核型的异常值。杠杆值和标准误差是相关的，因此高杠杆值应该具有高(不显著)的p值。杠杆值可以报告为n和p的函数：
[0052]
[0053]
其中p是模型中协变量的数量。一方面，可以标记hi(n,p)值大于阈值的样本。例如，阈值可以为约3至约5。这通常可用于确保最优拟合。可以使用更保守的阈值来标记最极端的值，例如，对应于第99.5和99.9百分位数的值(～10和～26)。在某些情况下，去除高杠杆值样本并重新拟合模型很有用，从而减小了不具有高杠杆值的样本的标准误差并改善(降低)p值估计值。
[0054]
图4是说明了用于检测异常核型的示例方法400的流程图。在步骤402，可以测定多个样本中每条染色体的读段覆盖度数据。一方面，每条染色体可包含多个基因组区域。测定多个样本中每条染色体的读段覆盖度数据包括测定外显子组区域上gc含量在一定范围内并且可映射性分数高于阈值的读段深度的总和。
[0055]
方法400还可以包括过滤所述读段覆盖度数据。过滤所述读段覆盖度数据可包括基于所述多个基因组区域中一个或多个基因组区域内的鸟嘌呤-胞嘧啶(gc)含量水平过滤所述读段覆盖度数据。基于所述多个基因组区域中一个或多个基因组区域内的鸟嘌呤-胞嘧啶(gc)含量水平过滤所述读段覆盖度数据可包括测定所述多个基因组区域中每一个的gc含量水平并排除所述多个基因组区域中gc含量水平在一定范围之外的一个或多个基因组区域。
[0056]
一方面，本方法可以过滤出具有极端gc含量的一个或多个基因组区域。当偏倚对于任何特定水平的gc含量大致一致时，可以校正gc-扩增偏倚。然而，在极低或极高的gc含量下，覆盖度随机波动可能急剧增加，从而难以有效地归一化。因此，本方法可以过滤出gc分率在可配置(例如，或预定义)范围或阈值之外的一个或多个基因组区域。作为说明，可配置范围可以包括[0.3,0.7]，如图5所示。然而，应理解，可以视情况利用其他范围(例如阈值)。图5显示了说明gc含量和覆盖度的关系的图表。例如，y轴上显示覆盖的变异系数(例如，标准偏差除以平均值)而x轴上显示gc含量。该图表显示了50个样本(例如，为了可见而抖动的点)。高于可配置范围的默认上限(例如，gc＝0.7)，覆盖度方差相对于平均值可以非常高。低于可配置范围的默认下限(例如，gc含量＝0.3)，出现了其他问题。例如，样本间覆盖度本身的方差可以具高变异性。该方差使得难以精确估计特定样本在特定窗口的覆盖度期望方差，因为每个参照组样本的覆盖度值是来自不同分布的观测值。
[0057]
所述方法400中过滤所述读段覆盖度数据可包括基于所述多个基因组区域中一个或多个基因组区域的可映射性分数来过滤所述多个基因组区域中的所述一个或多个基因组区域。基于所述多个基因组区域中一个或多个基因组区域的可映射性分数来过滤所述多个基因组区域中的所述一个或多个基因组区域可包括测定所述多个基因组区域中每个基因组区域的可映射性分数并且如果所述多个基因组区域的一个或多个基因组区域的可映射性分数低于预定阈值，则排除所述一个或多个基因组区域。
[0058]
例如，本方法和系统可以过滤所述多个基因组区域中的所述一个或多个基因组区域，其中从窗口中的每个碱基开始的k聚体的平均可映射性分数小于0.75(默认k＝75)。测定所述多个基因组区域中每个基因组区域的可映射性分数可包括测定其第一碱基与所述多个基因组区域的所述基因组区域重叠的k聚体的反向参照基因组频率的平均值。
[0059]
一方面，方法400还可包括归一化所述读段覆盖度数据。归一化所述读段覆盖度数据可包括测定每条染色体相对于其他常染色体的全外显子组读段覆盖度比率。可以通过以下方式测定每条染色体(i)的全外显子组比率(γ)：
[0060][0061]
其中α为常染色体集合，r为读段覆盖度。
[0062]
在步骤404，可以测定所述多个样本中每条染色体的期望读段覆盖度数据。测定所述多个样本中每条染色体的期望读段覆盖度数据可包括应用线性回归模型测定每条染色体的全外显子组期望比率，其中将多个度量用作协变量。所述多个度量可包括测序质量控制度量(qc度量)。由于测序条件的可变性而产生的系统覆盖度偏倚通常称为“批次效应”。一方面，本方法和系统可以配置为校正批次效应。例如，不基于读段覆盖度曲线即高维空间来比较读段覆盖度数据，本方法和系统可以配置为基于测序质量控制(qc)度量来考虑低维度量空间。例如，测序qc度量可包括七个测序qc度量。测序qc度量可以包括来自测序工具(例如picard)的测序qc度量。在该低维空间中工作可以提高可扩展性。例如，可以提前为样本编索引(例如，使用任何适当的索引和/或搜索算法)。
[0063]
一方面，可以通过以下方式测定每条染色体(i)的全外显子组期望比率
[0064][0065]
其中染色体j、k定义为相对于染色体i的gc含量分布具有最小d统计值的两条常染色体并且εi是和γj，γk之间线性关系的随机分量。
[0066]
在步骤406，可以测定所述多个样本中至少一条染色体的读段覆盖度数据和期望读段覆盖度数据之间的偏差。测定所述多个样本中至少一条染色体的读段覆盖度数据和期望读段覆盖度数据之间的偏差可以包括为所述多个样本中的每条染色体，测定读段覆盖度数据和期望读段覆盖度数据之间的差异以生成多个残差并且对于所述多个样本中具有协变量x的单个样本，相对于平均估计值的标准误差对所述多个残差进行z分数归一化：
[0067][0068]
其中se为残差标准误，并且：
[0069][0070]
参见图6，其描绘了使用线性回归模型获得的结果，其中协变量包括qc度量和染色体，并且其中6显示了线性回归拟合后的观测值(γi)和期望值另一方面，可以使用不同的标准误差估计器，例如原始残差标准误(整个模型的一个值)或使用异方差一致标准误。
[0071]
方法400还可包括测定每条染色体基于z分数的p值以鉴定显著较大的残差，其表示染色体i的异常核型。显著较大的残差可包括p值小于0.05的残差。参见图6。
[0072]
在步骤408，可以将所述至少一条染色体鉴定为异常核型。可以输出鉴定的一个或多个异常核型。例如，可以将鉴定的一个或多个异常核型输出给用户(例如，经由用户界
面)。可以经由网络将鉴定的一个或多个异常核型传输到远程位置。可以提供鉴定的一个或多个异常核型作为另一可执行程序的输入。鉴定的一个或多个异常核型可以存储在存储位置，例如数据库或其他文件格式。图7-10中示出了示例输出。
[0073]
图7a-f是等位基因平衡图，其显示了9、13和20号染色体的部分染色体等位基因平衡事件。子图编号为染色体编号。阴影条701表示杂合snp等位基因平衡为0.5的期望正常变异范围。线702表示全染色体中值等位基因平衡。虚线703表示约20snp滚动窗口中的局部中值等位基因平衡。线704表示连续性纯合片段。图8是读段覆盖度图，其显示相同样本13和20号染色体读段的显著代表性不足。
[0074]
图9a-f是21三体样本(唐氏综合征)的等位基因平衡图。等位基因平衡图显示了样本中21号染色体上的异常和涵盖整个x染色体的连续性纯合片段，暗示核型正常的雄性样本仅有一条x染色体。阴影条(901)表示杂合snp等位基因平衡为0.5的期望正常变异范围。实线(902)表示全染色体中值等位基因平衡。虚线(903)表示约20snp滚动窗口中的局部中值等位基因平衡。线(904)表示连续性纯合片段。图10是相同样本的读段覆盖度图。
[0075]
使用本文公开的方法获得的信息可以由临床医生报告给患者，例如以便提供对现有诊断，例如自闭症或自闭症谱系病状的进一步临床见解。
[0076]
使用本文公开的方法获得的信息也可以由临床医生用于为患者明确已知或未知生育问题，例如在具有性染色体异常的患者中。
[0077]
本文公开的方法也可用于监测癌症检测和发展。
[0078]
本文公开的方法还可用于确定dna样本是否含有来自于两个个体的dna，例如，如果来自一个个体的dna样本被来自另一个体的dna污染，则可能发生所述情况。发生双胞胎死亡/人类嵌合体事件时，dna也可以来自两个个体，即多胎妊娠受孕，其中并非每个胎儿都存活并且已故双胞胎的dna并入到存活胎儿的dna中。在此类情况下，对于其中双胞胎dna不相同的基因组的所有区域而言，结果将是偏态、多峰等位基因平衡，所述基因组是异卵双胞胎基因组的约75％。将来自一个个体的血液或组织移植到另一个体中时，dna也可来自于两个个体。当获得非侵入性产前检测样本时，当母胎dna混合时，dna也可来自于两个个体。
[0079]
回到图4，方法400还可包括检测一个或多个异常值并将所述一个或多个异常值去除，不考虑鉴定为异常核型。检测一个或多个异常值可包括对于每条染色体而言在线性回归模型上标记所述多个样本中具有高于阈值的杠杆值(hi，其中1/n《hi《1)的一个或多个样本，其中杠杆值作为n和p的函数得以确定：
[0080][0081]
其中p为所述模型中协变量的数量，xi表示样本i的协变量向量，并且为样本群体中协变量平均值的向量。所述阈值可以为约3至约5。
[0082]
可以由在本文的karyoscan方法之前生成的基因组比对序列读段计算读段覆盖度数据，以便检测单个样本的单核苷酸多态性(snp)、插入和缺失(插入缺失)。仅有两个观测等位基因的snp(或一个观测纯合等位基因不同于参照基因组定义的在该特定样本的序列读段中未被观测到的第二等位基因)称为双等位基因snp。通过关注双等位基因snp，可以计算基因组中特定位点的等位基因平衡。
[0083]
另一方面，等位基因平衡分析可用于鉴定一种或多种核型。等位基因平衡是对支持每个等位基因的序列读段数量的量度。例如，如果杂合snp被100个序列读段覆盖，并且样本在该基因组区域中为二倍体，则可以预计50个读段为一个等位基因和50个读段为另一等位基因，从而产生0.5/0.5的等位基因平衡。因为两个等位基因的等位基因平衡总和为1并且关于约0.5对称，所以焦点在于次要等位基因平衡(例如，具有较少读段的等位基因，或者如果两个等位基因在覆盖度上完全相等则随机选择的等位基因)。实际上，观测到的等位基因平衡很少恰好为50％，但是会遵循概率分布，该概率分布反映了在给定真实比例p的情况下，在大小为n的样本(n＝比对序列读段的数量)中出现的每个等位基因的读段数量。理想地，二倍体样本中的杂合snp的p＝0.5，并且等位基因平衡可以用二项分布建模，期望值为0.5。
[0084]
在具有非二倍体区域(例如，21三体综)的样本中，非二倍体区域中的双等位基因杂合snp将不具有0.5的期望等位基因平衡。如果一条染色体是重复的，例如对于21三体综合征而言，则2/3的21号染色体拷贝将具有一个等位基因，并且1/3的21号染色体拷贝将具有另一个等位基因，从而产生～0.333的期望等位基因平衡。因此，通过对整条染色体上具有一定程度集中趋势的等位基因平衡分布建模，可以通过确保相应等位基因平衡收敛于大致0.333来从读段深度模型验证21三体调用。可以使用诸如染色体上的中值等位基因平衡估计值的度量。类似地，对于单体染色体，只能存在一个等位基因，并且不会鉴定出杂合snp。因此，等位基因平衡将为0或完全未观测到，并且仅可鉴定出纯合snp(半合子)。可以经由连续性纯合片段来鉴定这些区域。
[0085]
这两个实例都假定全染色体重复或缺失。然而，在等位基因平衡分布中也可以观测到部分染色体重复和缺失。为了区分部分染色体事件，可以使用集中趋势的局部估计并且鉴定该局部估计与染色体其余部分的偏差。实际上，等位基因平衡的方差与覆盖snp的读段数量成比例，并且必须使局部估计在足够数量的位点上平滑以减小由各个位点贡献的总方差。为了实现这种平滑，可以计算20个杂合双等位基因snp的窗口上的滚动中值。该窗口大小可以根据测序深度增大或减小，这是因为由于样本量增加，测序越深，每个特定位点的方差越低。类似地，可以鉴定出仅跨越部分染色体的连续性纯合片段。
[0086]
除了部分染色体事件外，在等位基因平衡分布偏差中还会反映镶嵌事件(全染色体或部分染色体)。镶嵌事件是在为测序样本提供dna的细胞群子集中发生的事件。镶嵌现象可能是体细胞突变(如癌症)或早期种系细胞分裂错误的结果。例如，如果仅在50％的测序细胞中发生全染色体缺失，则来自缺失染色体的杂合snp除了读段覆盖度的25％损耗外还将具有25％期望等位基因平衡。因此，等位基因平衡也可用于区分镶嵌事件。
[0087]
并非所有异常核型都产生不同数量的染色体。例如，当染色体具有两个来自同一亲本的拷贝而没有来自另一亲本的拷贝时，出现单亲二体(upd)。在读段覆盖度偏差中不会检测到这些事件，但可以从杂合等位基因平衡(如果事件为镶嵌型)或从连续性纯合片段(如果事件非镶嵌型)中鉴定出。
[0088]
也可以发生染色体覆盖度的异常，不会引起等位基因平衡的异常。例如，如果染色体复制成四个拷贝(四体)，则所得核型可以具有每个亲本来源的两条染色体，产生～50％的正常等位基因平衡。这在镶嵌和非镶嵌事件中会具有相同效应。
[0089]
图11是说明用于检测异常核型的示例性方法1100的流程图，其并入了读段覆盖度
和等位基因平衡分析。方法1100可以测定为了方便在此处进行描述并且在方法流程的描述中提及的一个或多个度量。方法1100可以测定变体等位基因平衡，其可以是通过计算最小值(备选等位基因读段数量，参照等位基因读段数量)/总读段数量而确定的变体特异性度量。一方面，方法1100可以利用来自一个或多个vcf文件的“ad”(等位基因深度)和“dp”(读段深度)标签来测定变体等位基因平衡。
[0090]
方法1100可以测定可调用的染色体长度，其可以是通过计算染色体上第一个和最后一个未过滤外显子之间的碱基对数量即重叠着丝粒碱基的数量而确定的染色体特异性度量。对着丝粒碱基的调整是看似很大的事件跨越着丝粒的原因，其中不存在读段覆盖。实际上，基因组着丝粒边界可以调整到最近的外显子边界。类似地，限制为第一个和最后一个未过滤的外显子是长端粒区没有外显子覆盖，以及染色体整条臂缺乏外显子覆盖(例如，许多近端着丝粒染色体)的原因。
[0091]
方法1100可以确定全染色体杂合等位基因平衡(称为chromhetab)，其是染色体特异性度量，能够过滤推定的杂合snp，由此变体等位基因平衡》0.02(可以根据测序深度将阈值调整为与0接近或相差更大)。chromhetab可以是表示染色体内所有未过滤变体间的全染色体杂合snp等位基因平衡的汇总统计值(例如中值)。例如，chromhetab可以通过计算染色体内所有未过滤变体的中值(变体等位基因平衡)来确定。提到关于特定snp、localhetab事件或roh事件的chromhetab可以是指发生snp或事件的染色体的chromhetab值。chromhetab可以是表示染色体内所有未过滤变体间的全染色体杂合snp等位基因平衡的汇总统计值(例如中值)。
[0092]
方法1100可以确定局部中值杂合等位基因平衡(称为localhetab)，其是变体特异性度量，能够过滤可能的杂合snp，由此变体等位基因平衡》0.02(可以根据测序深度将阈值调整为与0接近或相差更大)。可以通过使用20snp窗口和恒定末端计算整条染色体上变体等位基因平衡的动态中值来确定localhetab。一方面，确定localhetab可以包括确定染色体上所有未过滤变体的样本杂合snp等位基因平衡的平滑、亚染色体规模(例如局部)的汇总统计值(例如动态中值)。
[0093]
方法1100可以确定两个或更多个snp，全部具有localhetab《chromhetab(称为localhetab事件)的连续区域。方法1100可以通过localhetab事件内的第一个和最后一个snp来定义坐标(染色体起始和结束位置)。每条染色体可以有零到多个localhetab事件。方法1100可以通过计算localhetab事件的归一化“曲线下面积”来确定localhetab事件面积。例如，对于localhetab事件内的一对相邻snp，确定成对面积＝[chromhetab-平均(localhetab(snp1)，localhetab(snp2))]*(snp2位置-snp1位置-重叠着丝粒碱基对数量)。在最小形式中，localhetab事件可以恰好具有两个相邻snp。具有两个以上snp的localhetab事件可以视为n-1个相邻snp对的链，其中n＝事件中的snp数量。具有两个或更多个snp的localhetab事件可以通过计算(localhetab事件中所有n-1个相邻snp对的成对面积)总和/(可调用染色体长度*chromhetab)来确定。
[0094]
方法1100可以通过确定最小值(localhetab事件中所有snp的localhetab)来确定localhetab事件的等位基因平衡(ab)汇总统计值(称为localhetab事件ab)。因为localhetab是等位基因平衡的平滑(动态中值)估计值，所以最小值是对整个事件的良好估计。然而，替代度量(例如，平均值、中值、第1四分位数等)可能更适合于其他应用(例如更大
的snp窗口大小、更深的测序、全基因组测序等)。
[0095]
方法1100可以确定连续性纯合片段(称为roh)，其是观测到很少杂合性甚至无杂合性的染色体区域的变体特异性度量。roh是每个变体的二元(是/否)标记，但是可以具有支持度量(例如，置信度分数)。一方面，确定roh可以包括使用通过引用整体并入本文的narasimhan,v.等(2016)bioinformatics,32(11),1749
–
1751)描述的bcftools/roh方法。roh确定的示例选项包括但不限于：autozygous-to-hardy weinberg转移概率(-a选项)＝6.6e-09、hardy weinberg-to-autozygous转移概率(-h选项)＝5.0e-10、忽略插入缺失(-i选项)、限于外显子内的snp(即无侧翼区域snp)以及利用内部rgc(eve)变体频率。一方面，可以使用一种或多种替代方法。例如，如purcell s,neale b,todd-brown k等plink:a tool set for whole-genome association and population-based linkage analyses.american journal of human genetics.2007；81(3):559-575描述的plink，其通过引用整体并入本文。
[0096]
方法1100可以确定预测为roh(称为roh事件)的一个或多个snp的连续区域。事件坐标可以定义为roh事件内第一个和最后一个snp的染色体位置。
[0097]
回到图11，在方框1102处可以对所有样本的数据进行质量控制(qc)过滤。数据可以包括例如vcf文件(例如每个样本一个vcf文件)、覆盖深度文件和/或外部质量控制度量(例如由bam读段映射文件计算的picard度量)。vcf文件可包含基因序列变异的标志和基因型数据。覆盖深度文件可以包含许多读段的指示，所述读段包括给定核苷酸或核苷酸序列。qc过滤可以包括将一个或多个样本过滤标准应用于覆盖深度文件、vcf文件和/或外部质量控制度量。所述一个或多个样本过滤标准可包括例如：标准污染过滤(例如高杂合子与纯合子snp调用比)、基于低序列覆盖度(20x覆盖度或更高覆盖度下《75％的碱基)的过滤和/或基于低dna质量的过滤、其组合等。一方面，qc过滤可包括将一个或多个变体过滤标准应用于vcf文件。所述一个或多个变体过滤标准可以包括例如：仅分析双等位基因snp(去除多等位基因位点和插入缺失)、基于最低变体质量的过滤(qd》5，gt》30，通过vqsr过滤器[变体质量得分重新校准])、基于最小读段深度(dp》＝20)的过滤、和/或基于基因座质量的过滤(1.仅可映射性》90％的外显子，2.排除常见拷贝数》2的外显子(例如多拷贝cnv基因座)，3.排除具有可映射性问题的其他区域(例如hla基因))、其组合等。
[0098]
在方框1104处，可以对在方框1102处通过qc过滤的样本相关的数据进行性别分配。性别分配可以包括确定最小y染色体读段覆盖度比率(相对于x染色体读段覆盖度比率)以确定样本是雄性(高于阈值)还是雌性(低于阈值)。图13a是所有样本的x染色体相对于y染色体覆盖度比率的图和实线1306所示用于确定雄性(1302)和雌性(1304)样本的阈值。另外，可以使用y染色体覆盖度比率阈值(虚线1308)鉴定有y染色体重复的雄性样本。如果样本性别已知或已为样本报告性别，则可以使用现有分配来帮助确定适当的阈值。在方框1104处的性别分配之后，可以通过方法1100的剩余一个或多个方框处理来自一个或多个样本的每条染色体。
[0099]
如果将样本视为雄性，则方法1100将进行到方框1106。在方框1106处，方法1100可以确定y染色体的覆盖度是否大于阈值，例如0.0015。如果y染色体的覆盖度大于阈值，则方法1100可以在方框1108处确定存在y染色体重复并且进行到方框1138，这是因为y染色体可以独立于其他染色体进行处理。如果y染色体的覆盖度小于阈值，则方法1100可以在方框
1108处确定样本具有雄性样本的正常剂量的y染色体读段，因此在y染色体上未出现可检测的异常。
[0100]
回到方框1104，性别分配可以包括确定是期望样本具有一条还是两条x染色体(雄性或雌性)，在这种情况下，方法1100将进行到方框1110以处理样本的x染色体。在方框1110处，方法1100可以确定数据是否源自雄性。如果在方框1110处，确定数据是源自雄性，则方法1100将前进到方框1112和1114。如果在方框1110处，确定数据不是源自雄性，则方法1100将进行到方框1112、1114、1116和1118。回到方框1104，性别分配可包括确定数据包含常染色体，在这种情况下方法1100将进行到方框1112、1114、1116和1118。
[0101]
在方框1112处，方法1100可以检测读段覆盖度异常。方框1112可以如本文所述，参考图1和/或图4的一部分来进行。在方框1114处，方法1100可以检测chromhetab异常。在方框1116处，方法1100可以检测roh异常。在方框1118处，方法1100可以检测localhetab异常。
[0102]
方框1114、1116和1118涉及测定三个等位基因平衡度量(分别为chromhetab、roh和localhetab)。这三个等位基因平衡度量可用于检测不同类型的异常，但可能导致重叠的证据。例如，roh可用于鉴定组成型染色体缺失(全染色体或部分染色体)，因为在这些区域中不应观测到杂合性。类似地，roh可以鉴定大的单亲二体(upd)事件(复制中性，全染色体或部分染色体)，但对于鉴定重复无效。然而，localhetab和chromhetab度量也可以通过鉴定少量类似于推定杂合性的噪声(由于技术误差，例如测序错误而引起)而在roh事件内产生异常信号，其中变体等位基因平衡值接近于0；可以忽略这些信号来代替roh异常(参见图12，表示具有重叠localhetab事件的roh事件)。图12是样本4号染色体上的示例等位基因平衡图，其中检测到大的连续性纯合片段(1202)，其由于异常区域中纯合snp间有少量非零等位基因平衡而具有重叠localhetab事件(1204)。在全染色体重复或其他镶嵌型全染色体事件的情况下，chromhetab可以是最相关的度量；对于三体而言应大致等于1/3，或者是代表镶嵌事件的拷贝数和细胞分率的分率。对于部分染色体事件，localhetab可以是最相关的度量，因为它将检测事件起始和结束坐标。然而，大的部分染色体事件也会影响全染色体chromhetab度量，从而产生通过localhetab事件更好捕获的异常信号。
[0103]
因此，平衡由每个度量提供的证据(并相对于读段覆盖度异常信号对每一个进行解释)可以是使染色体异常的检测和表征自动化的组成部分。为了处理这些在灵敏度、特异性、规模和范围方面有差异的集成信号，可以为每个度量定义三个级别的推定异常信号，其中1级信号最显著，3级最不显著。层级评级用于标准化和集成这些异构度量，从而能够简单地决定哪些信号最相关。可以使用和定义其他数量的层级。
[0104]
回到方框1112，读段覆盖度异常的检测可以利用以下层级定义。1级可以包括读段覆盖度p值《阈值，例如0.05/(检验的染色体/样本对数量)。可以按族系误差率＝5％应用邦费罗尼校正(bonferroni-correction)。2级可以包括未通过1级且染色体特异性fdr校正p值(q值)《阈值，例如0.05。可以按每条染色体错误发现率＝5％应用benjamini-hochberg fdr校正。3级可以包括未通过1级或2级且读段覆盖度p值《阈值，例如0.05。一项或多项例外可适用于x染色体分析。例如，如果估计的染色体剂量分率的绝对值(大小)》5％，则可以将x染色体上的3级信号提升至2级。
[0105]
回到方框1114，因为变体等位基因平衡度量总会反映来自具有较少读段的双等位基因snp的等位基因的读段的分率，所以核型正常的二倍体样本的指定染色体上的期望
chromhetab值不可能恰好为50％，而是随着测序深度的增加接近50％。因此，相对于pcttargetbases50x质量控制度量(使用picard为每个样本计算的)，可以对所有样本指定染色体(雌性仅针对x染色体)的chromhetab值进行线性回归拟合(图13b，显示核型正常样本中chromhetab值相对于pcttargetbases50x值增加而增加，并且鉴定出不同显著性层级的异常信号[有色点])。线性回归模型一经拟合，就可以计算每个样本chromhetab值残差的z分数(如回归所定义的chromhetab观测值-chromhetab期望值)。z分数可以按照z＝(样本残差)/(回归模型的残差标准偏差)计算。z分数可以转换为p值。
[0106]
在方框1114处，chromhetab异常的检测可以利用以下层级定义。1级可以包括chromhetab残差p值《阈值，例如0.05/(检验的染色体/样本对数量)。可以按族系误差率＝5％应用邦费罗尼校正。2级可以包括未通过1级而通过染色体特异性fdr校正p值(q值)《阈值，例如0.05。可以按每条染色体fdr校正错误发现率＝5％应用benjamini-hochberg。3级可以包括未通过1级或2级且chromhetab残差p值《阈值，例如0.05。一项或多项例外可适用于x染色体分析。如果样本为雄性，则可以忽略chromhetab并且不进行检验，除非度量计算中包括》75个snp且chromhetab》0.15。这些过滤器允许包括雄性中的x染色体重复，只要它们具有比期望值(即，对于单条x染色体而言为零)大得多的chromhetab值并且具有足够数量的snp用于确信地调用chromhetab值，同时去除了来自核型正常的雄性样本的噪声(图14，显示了在雄性样本x染色体上检验到的chromhetab值和推定杂合snp数量，实线表示最小阈值)。在这种情况下，定义雄性样本可以是指(基于x和y读段覆盖度)分配为具有一个x染色体和一个y染色体的期望，假设是核型正常状态。通过这些过滤器的来自于x染色体的任何雄性chromhetab信号都可以分配为1级(与p值无关)。一方面，检测chromhetab异常可以包括鉴定chromhetab值显著小于核型正常样本的期望chromhetab值(或值的范围)的样本。
[0107]
回到方框1116，可以检测roh异常。小roh事件在核型正常样本中相对常见，并且例如在血亲样本中可能特别频繁。因此，可以定义roh事件的大小最小的阈值以便仅捕获大的染色体规模的事件。真纯合变体由于技术误差而具有非零变体等位基因平衡时，roh事件检测可能具有挑战性。因此，一些大的roh事件被分割成两个或更多个roh事件(图12)。因此，组合考虑染色体内的独立roh事件。roh异常的检测可以利用以下层级定义。可以过滤长度(不包括重叠着丝粒碱基)《5,000,000bp的roh事件。1级可包括来自未过滤的roh事件的非着丝粒roh碱基(全基因组)总数》＝20,000,000。2级可包括未通过1级的未过滤的roh事件。一项或多项例外可适用于x染色体分析。雄性样本x染色体上的所有roh信号都可以忽略。在这种情况下，定义雄性样本是指(基于x和y读段覆盖度)分配为具有一个x染色体和一个y染色体的期望，假设是核型正常状态。
[0108]
回到方框1118，可以检测localhetab异常。定性地，显著localhetab事件异常应具有大的localhetab事件面积度量并且应受大量所包括的snp支持。可以定义线性函数以根据经验拟合为关联localhetab事件面积和在localhetab事件中所包括的snp数量(“snp数量”)的外显子组数据集，其中层级定义在相同斜率系数上使用不同截距进行定义(例如，具有特定localhetab事件面积的事件所需的最小snp数量)。见图15，显示了面积》0.02(垂直线)的所有localhetab事件(点)以及它们基于对角线之间的区域分级(红点表示存在重叠roh事件，提供了localhetab事件检测到异常的支持证据)。localhetab异常的检测可以利
用以下层级定义。可以过滤localhetab事件面积《0.02的事件。1级可以包括snp数量+(localhetab事件面积*第一量，例如3000)》＝第二量，例如230。2级可以包括snp数量+(localhetab事件面积*第一量如3000)》＝第二量如170。3级可以包括snp数量+(localhetab事件面积*第一量如3000)》＝第二量如110。一项或多项例外可适用于x染色体分析。雄性样本x染色体上的所有localhetab信号均被忽略。在这种情况下，定义雄性样本是指(基于x和y读段覆盖度)分配为具有一个x染色体和一个y染色体的期望，假设是核型正常状态。一方面，检测localhetab异常可以包括其中localhetab值低于(例如，显著低于)染色体区域上的相应chromhetab值的情况，表明可能部分染色体异常。
[0109]
对于所公开的度量，在方框1112、1114、1116和1118处进行的分析有助于染色体异常预测。然而，度量可以用层级评级来注释，过滤去除非异常度量，并且在方框1136处进行聚合，之后在方框1138处进行核型预测。在方框1120处，方法1100可以报告通过来自方框1112、1114、1116和1118的每个度量鉴定的异常事件，并且将每个事件评分为标准化度量之间的比例并简化其聚合以进行异常核型(染色体异常)评估的层级(例如，1级、2级、3级等
…
)。一方面，对于所使用的所述一个或多个层(例如，1级、2级、3级等
…
)中的每一个而言，方法1100都可以对事件进行报告和/或评分。在方框1112处，读段覆盖度异常度量可用于评估染色体剂量，其余三个可用于评估等位基因平衡和接合性(chrommedab、roh和localhetab事件)。
[0110]
在方框1122处，方法1100可以确定事件反映了拷贝增益、拷贝损失还是拷贝中性。该评估可以主要基于读段覆盖度异常存在与否来进行，但是也可以考虑来自等位基因平衡相关度量的补充信息。例如，可以独立地将所有1级读段覆盖度异常预测为增益或损失，但是如果在同一染色体上也检测到支持等位基因平衡异常，则可以仅将2级和/或3级读段覆盖度异常视为增益或损失。如果未检测到读段覆盖度异常来调用增益或损失，则假定该事件为拷贝中性，并且如果检测到但是过滤了低质量读段覆盖度异常，则可以另外标记为不确定。
[0111]
如果在方框1122处，确定事件反映了拷贝增益，则方法1100以可进行到方框1124，以基于染色体上异常localhetab和chromhetab事件的比较来确定潜在染色体是全染色体还是部分染色体。例如，如果localhetab事件具有比重叠chromhetab事件更低的层级评级(例如，更显著)(或者如果未报告chromhetab事件)，则假定localhetab事件更显著，可以将该事件预测为部分染色体。相反，低级chromhetab事件将表明全染色体事件更有可能是这种情况。如果两个事件都发生在相同的层级评级，则该方法可以将异常报告为不确定和/或偏好一个事件为加权更重(例如，偏好localhetab并调用部分染色体事件)。此外，该方法还可以将由读段覆盖度11计算的染色体分率增益的估计值与来自每个等位基因平衡异常的类似估计值进行比较，并且按照其估计值与来自读段覆盖度的估计值的匹配紧密程度来对等位基因平衡事件进行加权。无论预测潜在染色体异常是全染色体(chromhetab)还是部分(localhetab)染色体，方法1100可以进行到方框1126并利用相应的等位基因平衡度量，通过确定杂合等位基因平衡估计值有多接近1/n来确定拷贝增益是否是镶嵌事件，其中n为预测的染色体拷贝数(例如对于单拷贝常染色体增益，为1/3)。可以应用在该期望比率(例如，1/3
±
0.02)左右的误差阈值，以对镶嵌现象进行二元(是或否)分类。如果未报告重叠localhetab和chromhetab事件，可以将染色体和镶嵌分率估计值分配为不确定和/或设置
默认值。
[0112]
如果在方框1122处，确定事件反映了拷贝损失，则方法1100以可进行到方框1128，以利用roh确定拷贝损失是否为镶嵌事件。如果拷贝损失不是镶嵌的(例如，检测到roh事件)，则方法1100可以利用roh通过评估roh事件所覆盖的可调用染色体的比例来确定潜在染色体异常是全部还是部分的。如果拷贝损失是镶嵌的(例如，未报告roh事件)，则方法1100可以利用并比较来自报告的chromhetab和localhetab事件的层级评级，以确定潜在染色体是全部还是部分。该评估类似于拷贝增益的评估(方框1124)，其中更显著的localhetab事件可指示部分染色体事件且更显著的chromhetab事件可指示全染色体事件，并且可以将基于等位基因平衡的染色体分率估计值与读段覆盖事件的染色体分率估计值进行比较。
[0113]
如果在方框1122处，确定事件为拷贝中性，则方法1100以可进行到方框1128，以利用roh确定拷贝中性事件是否是镶嵌的。如果拷贝中性事件不是镶嵌的(例如，报告了roh事件)，则方法1100可以利用roh来确定潜在染色体是全部还是部分。如果拷贝中性事件是镶嵌的，则方法1100可以利用chromhetab和localhetab来确定潜在染色体是全部还是部分。
[0114]
方框1126、1130和1134的输出流到方框1136，其中可以为每个异常报告以下一项或多项：1)拷贝中性、拷贝增益或拷贝损失的预测；2)全染色体或部分染色体事件的预测；3)镶嵌或非镶嵌的预测；4)最终层级评级，对于染色体上报告的所有事件而言可以等于最低(最显著的)层级评级，或者如果它们有多个中级事件，则可以另外修改为更高或更低权重的异常(例如，具有支持性2级localhetab事件的2级读段覆盖事件可视为1级)；5)对染色体报告的一些或所有事件的汇总、其层级评级以及它们是否被选为主要或支持事件(例如，对于具有1级读段覆盖度、roh和chromhetab事件的非镶嵌全染色体丢失，读段覆盖度和roh是主要事件，但chromhetab尽管为1级，考虑到它被重叠roh事件所胜过，却是支持事件)。方框1136接收零条或多条染色体的异常并针对一个样本将其进行聚集，然后进行到方框1138以进行最终核型预测。
[0115]
方框1138的输出表示核型预测，其中已经针对样本聚集了一些或所有染色体异常并且相对于期望核型进行了解释(鉴于来自方框1104的性别分配)。这可以表示为传统核型编码(例如，“47,xxy”)和/或异常列表及其支持信息。鉴于自动化核型预测的不确定性以及某些复杂核型(例如，等臂染色体)具有不易自动解释的独特模式的事实，可以通过方框1138为每个样本计算支持读段覆盖度和等位基因平衡诊断图，从而能够手动检查预测染色体异常及其支持证据。一方面，最终异常核型调用可以包括以下一项或多项：样本；染色体；起始/结束坐标；剂量变化与拷贝中性预测(增益、损失、中性、不确定)；全染色体与部分染色体事件预测(全部、部分、不确定)；预测镶嵌事件(是、否、不确定)；读段覆盖度的分率估计值(即染色体分率*镶嵌分率，其中单拷贝、非镶嵌染色体增益＝1，或损失＝-1)；等位基因平衡的分率估计值(基于异常等位基因平衡度量，如果存在一个以上则视为最相关)；该样本/染色体对的3级或更高层级的所有异常信号汇总；最终综合层级评级；允许手动检查核型和分类的支持读段覆盖度和等位基因平衡诊断图(例如图2、图3、图5、图6、图7、图8、图9和/或图10中所示的类型)；其组合等。
[0116]
在示例性方面，所述方法和系统可以在如图16所示和下面所述的计算机1601上实现。类似地，公开的方法和系统可以利用一个或多个计算机在一个或多个位置执行一个或
多个功能。图16是说明用于执行所公开的方法的示例性操作环境的方框图。该示例性操作环境仅为操作环境的实例，并非旨在对操作环境结构的使用范围或功能提出任何限制。也不应将操作环境解释为对示例性操作环境中示出的任一部件fig或其组合有任何依赖性或要求。
[0117]
本方法和系统可以用许多其他通用或专用计算系统环境或配置操作。可以适于和所述系统和方法一起使用的公知计算系统、环境和/或配置的实例包括但不限于个人计算机、服务器计算机、膝上型计算机设备和多处理器系统。另外的实例包括机顶盒、可编程消费性电子产品、网络pc、小型计算机、大型计算机、包括上述任何系统或设备的分布式计算环境等。
[0118]
公开的方法和系统的处理可以由软件部件执行。公开的系统和方法可以在由一个或多个计算机或其他设备执行的计算机可执行指令(例如程序模块)的一般背景下进行描述。通常，程序模块包括执行特定任务或实现特定抽象数据类型的计算机代码、例程、程序、对象、部件、数据结构等。还可以在基于网格的分布式计算环境中实践公开的方法，其中任务由通过通信网络链接的远程处理设备执行。在分布式计算环境中，程序模块可以位于包括存储器存储设备的本地和远程计算机存储介质中。
[0119]
此外，本领域技术人员将理解，本文公开的系统和方法可以通过呈计算机1601形式的通用计算设备来实现。计算机1601的部件可以包括但不限于一个或多个处理器1603、系统存储器1612以及将包括所述一个或多个处理器1603的各个系统部件耦合到系统存储器1612的系统总线1613。该系统可以利用并行计算。
[0120]
系统总线1613代表几种可能类型的总线结构中的一种或多种，包括使用各种总线结构中的任一种的存储器总线或存储器控制器、外围总线、加速图形端口或本地总线。总线1613和本说明书中指定的所有总线也可以通过有线或无线网络连接实现并且每个子系统，包括所述一个或多个处理器1603、大容量存储设备1604、操作系统1605、karyoscan软件1606、karyoscan数据1607、网络适配器1608、系统存储器1612、输入/输出接口1610、显示适配器1609、显示设备1611和人机界面1602，均可包含在一个或多个远程计算设备1614a、b、c中的物理分隔位置，通过这种形式的总线连接，实际上实现了全分布式系统。
[0121]
计算机1601通常包括各种计算机可读介质。示例性可读介质可以是计算机1601可存取的任何可用介质，并且包括，例如但不意味着限制，易失性和非易失性介质、可移动和不可移动介质。系统存储器1612包括易失性存储器形式的计算机可读介质，例如随机存取存储器(ram)，和/或非易失性存储器形式的计算机可读介质，例如只读存储器(rom)。系统存储器1612通常包含数据诸如karyoscan数据1607和/或程序模块诸如操作系统1605和karyoscan软件1606，其可由一个或多个处理器1603立即存取和/或当前由其操作。karyoscan数据1607可以包括读段覆盖度数据和/或期望读段覆盖度数据。
[0122]
另一方面，计算机1601还可以包括其他可移动/不可移动、易失性/非易失性计算机存储介质。举例而言，图16说明了大容量存储设备1604，其可以为计算机1601提供计算机代码、计算机可读指令、数据结构、程序模块和其他数据的非易失性存储。例如并且不意味着限制，大容量存储设备1604可以是硬盘、可移动磁盘、可移动光盘、磁带盒或其他磁存储设备、闪存卡、cd-rom、数字通用盘(dvd)或其他光存储器、随机存取存储器(ram)、只读存储器(rom)、电可擦可编程只读存储器(eeprom)等。
包括以用于存储诸如计算机可读指令、数据结构、程序模块或其他数据等信息的任何方法或技术实现的易失性和非易失性、可移动和不可移动介质。示例性计算机存储介质包括但不限于ram、rom、eeprom、闪存或其他存储器技术、cd-rom、数字通用盘(dvd)或其他光存储器、磁带盒、磁带、磁盘存储器或其他磁存储设备，或可用于存储所需信息并且可由计算机存取的任何其他介质。
[0129]
所述方法和系统可以采用人工智能技术，例如机器学习和迭代学习。此类技术的实例包括但不限于专家系统、案例推理、贝叶斯网络(bayesian networks)、基于行为的人工智能、神经网络、模糊系统、进化计算(例如遗传算法)、群体智能(例如蚂蚁算法)和混合智能系统(例如通过神经网络生成的专家推理规则或来自统计学习的产生式规则)。
[0130]
本文的karyoscan方法使用新型共归一化技术，该技术在其gc含量和测序性能的背景下对染色体进行评估，从而可以实现更准确的覆盖度归一化。这与针对检测较小基因组变化的方法不同，因为它们完全取决于局部gc含量偏倚。虽然针对较小变化的方法有时可以检测多件较大事件，但通常使用平滑函数(例如，隐马尔可夫模型(hidden markov model))来理解在染色体臂规模下分解的较大事件背景下的高分辨率拷贝数变化。此外，将等位基因频率数据集成到karyoscan调用中提供了独特的特征，包括检测到平衡化基因组变化，其在覆盖空间中不呈现任何信号，但是由于遗传变异的丧失可能表现出显著影响。
[0131]
与为分率cnv，例如体细胞癌突变或镶嵌事件(即，仅在体内细胞子集中)强制或提供整数值调用的方法相反，本文的karyoscan方法提供了分率估计值。
[0132]
提出以下实例以向本领域普通技术人员提供对于如何制备和评价本文要求保护的化合物、组合物、制品、设备和/或方法的完整公开和描述，并且旨在纯粹为示例性的而非旨在限制所述方法和系统的范围。已努力确保有关数字(例如量等)的准确性，但应对一些误差和偏差进行说明。
[0133]
公开的方法应用于来自regeneron genetics center的人外显子组变体数据库的～100,000个样本。总共有3,150个样本在至少一条试验染色体上在最高严格性水平上被标记为核型异常，其中472个为增益或损失(不是拷贝中性)。超过200个样本被标记为具有性染色体异常(x染色体或y染色体)，包括极其罕见的核型(48，xxxx)和(48，xxxy)。
[0134]
虽然已经结合优选实施方案和具体实例描述了所述方法和系统，但是并非旨在将范围限于所阐述的特定实施方案，因为本文的实施方案在所有方面都旨在为说明性而非限制性。
[0135]
除非另有明确说明，否则决非旨在将本文阐述的任何方法解释为要求其步骤以特定顺序执行。因此，在方法权利要求实际上没有列举其步骤所遵循的顺序，或者在权利要求或说明书中没有另外特别说明将步骤限于特定顺序的情况下，其决非旨在在任何方面都可以推断出顺序。这适用于任何可能的非明确的解释基础，包括：关于步骤安排或操作流程的逻辑问题；从语法组成或标点符号中得出的简单含义；说明书中描述的实施方案的数量或类型。
[0136]
对于本领域技术人员显而易见的是，在不脱离范围或精神的情况下，可以进行各种修改和变化。考虑到本文公开的说明书和实践，其他实施方案对于本领域技术人员而言将显而易见。其意图是仅将说明书和实例视为示例性，真正的范围和精神由以下权利要求指出。
技术特征：1.一种方法，其包括：测定多个样本中每条染色体的读段覆盖度数据，其中每条染色体包含多个基因组区域；测定所述多个样本中每条染色体的期望读段覆盖度数据；测定所述多个样本中至少一条染色体的所述读段覆盖度数据和所述期望读段覆盖度数据之间的偏差；并且将所述至少一条染色体鉴定为异常核型。2.根据权利要求1所述的方法，其中测定多个样本中每条染色体的读段覆盖度数据包括测定外显子组区域上gc含量在一定范围内并且可映射性分数高于阈值的读段深度的总和，其中每条染色体包含多个基因组区域。3.根据权利要求1所述的方法，其还包括过滤所述读段覆盖度数据。4.根据权利要求3所述的方法，其中过滤所述读段覆盖度数据包括基于所述多个基因组区域中一个或多个基因组区域内的鸟嘌呤-胞嘧啶(gc)含量水平过滤所述读段覆盖度数据。5.根据权利要求4所述的方法，其中基于所述多个基因组区域中一个或多个基因组区域内的鸟嘌呤-胞嘧啶(gc)含量水平过滤所述读段覆盖度数据包括：测定所述多个基因组区域中每一个的gc含量水平；并且排除所述多个基因组区域中gc含量水平在一定范围之外的一个或多个基因组区域。6.根据权利要求3所述的方法，其中过滤所述读段覆盖度数据包括基于所述多个基因组区域中一个或多个基因组区域的可映射性分数来过滤所述多个基因组区域中的所述一个或多个基因组区域。7.根据权利要求6所述的方法，其中基于所述多个基因组区域中一个或多个基因组区域的可映射性分数来过滤所述多个基因组区域中的所述一个或多个基因组区域包括：测定所述多个基因组区域中每个基因组区域的可映射性分数；并且如果所述多个基因组区域中一个或多个基因组区域的可映射性分数低于预定阈值，则排除所述多个基因组中的所述一个或多个基因组区域。8.根据权利要求1所述的方法，其还包括归一化所述读段覆盖度数据。9.根据权利要求8所述的方法，其中归一化所述读段覆盖度数据包括测定每条染色体相对于其他常染色体的全外显子组读段覆盖度比率。10.根据权利要求9所述的方法，其中通过以下方式测定每条染色体(i)的全外显子组比率(γ)：其中α为常染色体集合，r为读段覆盖度。11.根据权利要求1所述的方法，其中测定所述多个样本中每条染色体的期望读段覆盖度数据包括应用线性回归模型测定每条染色体的全外显子组期望比率，其中将多个度量用作协变量。
12.根据权利要求11所述的方法，其中所述多个度量包括测序质量控制度量(qc度量)，并且通过以下方式测定每条染色体(i)的全外显子组期望比率方式测定每条染色体(i)的全外显子组期望比率其中染色体j、k定义为相对于染色体i的gc含量分布具有最小d统计值的两条常染色体。13.根据权利要求11所述的方法，其中测定所述多个样本中至少一条染色体的所述读段覆盖度数据和所述期望读段覆盖度数据之间的偏差包括：对于所述多个样本中的每条染色体，测定所述读段覆盖度数据和所述期望读段覆盖度数据之间的差异以生成多个残差；并且对于所述多个样本中具有协变量x的单个样本而言，相对于平均估计值的标准误差对所述多个残差进行z分数归一化：其中s
e
为残差标准误差，并且：14.根据权利要求13所述的方法，其还包括测定每条染色体基于z分数的p值以鉴定显著较大的残差，其表示染色体i的异常核型。15.根据权利要求14所述的方法，其中显著较大的残差包括p值小于0.05的残差。16.根据权利要求14所述的方法，其还包括：检测一个或多个异常值；并且将所述一个或多个异常值去除，不考虑鉴定为异常核型。17.根据权利要求16所述的方法，其中检测一个或多个异常值包括对于每条染色体而言在线性回归模型上标记所述多个样本中具有高于阈值的杠杆值(h
i
，其中1/n<h
i
<1)的一个或多个样本，其中杠杆值作为n和p的函数得以确定：其中p为所述模型中协变量的数量，n为建模样本的数量，x
i
表示样本i的协变量向量，并且x为样本群体中协变量平均值的向量。18.根据权利要求17所述的方法，其中所述阈值为约3至约5。
技术总结公开了用于检测异常核型的方法和系统。示例方法可包括测定读段覆盖度数据，杂合SNP的等位基因平衡分布，以及未观测到杂合性的染色体区段。然后所述方法和系统可以测定可以指示一个或多个异常核型的一个或多个度量。一个或多个异常核型的一个或多个度量。一个或多个异常核型的一个或多个度量。
技术研发人员：E
受保护的技术使用者：瑞泽恩制药公司
技术研发日：2017.02.13
技术公布日：2022/11/1
