CN112687336B - 确定upd类型的方法、计算设备和存储介质 - Google Patents

Abstract

本发明提供了一种确定UPD类型的方法、计算设备和计算机可读存储介质。该方法包括：基于关于先证者、先证者的父样本和先证者的母样本的测序序列数据，分别确定先证者、先证者的父样本和先证者的母样本的同一基因序列中的各个SNP位点信息；基于该SNP位点信息，确定每个SNP位点的第一UPD类型；基于每个SNP位点的第一UPD类型确定UPD区域；对该先证者的SNP位点进行AF过滤以确定该先证者的突变频率大于突变阈值的高频突变SNP位点；基于该先证者的高频突变SNP位点确定该先证者的ROH片段；以及基于该先证者的第一UPD类型和该先证者的ROH片段确定该先证者的第二UPD类型。

Description

确定UPD类型的方法、计算设备和存储介质

技术领域

本发明概括而言涉及生物信息检测领域，并且具体地，涉及一种用于确定UPD类型的方法、计算设备和计算机存储介质。

背景技术

单亲二倍体（Uniparental Disomy，UPD）是指一个个体的两条同源染色体都来自同一亲本。当UPD来自同一亲本的同一条染色体时，称为单亲同二体（iUPD或isoUPD）。当UPD来自同一亲本的一对同源染色体时，称为单亲异二体（hUPD或hetUPD）。此外，根据UPD来源亲本的不同，又分为母源UPD（matUPD）和父源UPD（patUPD）。由于减数分裂过程中交叉互换（crossover）造成染色体的某个区域产生的UPD又称为片段型UPD（segUPD）。

不是所有的UPD都有致病性，当UPD所覆盖的染色体区域存在隐性遗传的致病位点或印迹效应（imprinting effects）时，则可能导致各种疾病。例如，印迹区域15q11.2-q13的父源性缺失可导致Prader-Willi综合征（PWS），而当该区域发生母源性缺失时则可导致Angelman综合征（AS），因此确定片段型UPD对于确定是否存在致病风险非常重要。

随着技术的发展，UPD的检测手段包括微卫星标记（又称短串联重复序列STRs）分析、特异性甲基化检测、单核苷酸多态性阵列分析（SNP array analysis）、全外显子测序以及全基因组测序等。通过核型鉴定，能够发现因着丝粒融合导致的UPD，如罗氏易位。根据串联重复序列，也可以检测染色体来源，进而确定UPD，但受限于SSR（Simple SequenceRepeat，简单重复序列）的特殊位置，能够检测的UPD类型非常有限，不能检测片段型UPD。目前比较流行的UPD检测方法是芯片检测，基于SNP-microarray的数据开发了SNPtrio和UPDtools，能够很方便的看到大片段的纯合区域（runs of homozygosity，ROH），辅助确定iUPD的类型，当有家系样本时则可以鉴定出UPD。但是随着下一代测序（Next GenerationSequencing，NGS）技术的发展，测序成本下降。在遗传检测领域，全外显子组测序（WholeExome Sequencing，WES）已经发展为一种常规的检测技术。在该领域，先后出现了H3M2、UPDio等检测ROH和UPD的软件。然而，在检测UPD时，UPDio只能检测大于10M碱基长度的UPD，并且无法给出片段型UPD的具体位置。

在NGS测序中，针对WES测序数据，单个样本检测UPD只能借助大片段ROH来定位iUPD，而对于hUPD则无法检测。当存在trio家系时，利用孟德尔遗传错误（MIE）的方法可以找到UPD发生的染色体。例如，文献“Yauy, K., de Leeuw, N., Yntema, H.G. et al.Accurate detection of clinically relevant uniparental disomy from exomesequencing data. Genet Med (2019) doi:10.1038/s41436-019-0704-x”中介绍了利用H3M2识别ROH，使用中位数绝对偏差（MAD）对每条染色体的ROH大小进行归一化，通过MAD值的分布，将MAD值转换为P值，随后使用Bonferroni方法对P值进行校正，采用-log10 p >100作为阈值识别出潜在的iUPD。使用UPDio检测UPD，采用-log10 p > 48作为阈值。这个思路分析的iUPD依赖于ROH区域的大小，而H3M2检测到的ROH数量较多，导致iUPD的假阳性偏高。另外UPDio只能通过MIE来推断整条染色体发生UPD的P值，当发生了segUPD时，没有明确给出UPD的位置，需要用户通过原始数据判断。

发明内容

针对上述问题中的至少一个，本发明提供了一种基于ROH确定UPD类型的方法，通过将ROH检测与传统的UPD检测相结合能够更加准确地检测UPD类型，尤其是片段型UPD。并且，在一些实施例中，通过滑动窗的方式来确定UPD区域和ROH片段，提高了UPD类型检测的精度，从而有助于进一步准确确定UPD是否引起印记效应。

根据本发明的一个方面，提供了一种确定UPD类型的方法。该方法包括：基于关于先证者、先证者的父样本和先证者的母样本的测序序列数据，分别确定先证者、先证者的父样本和先证者的母样本的同一基因序列中的各个SNP位点信息；基于该SNP位点信息，确定每个SNP位点的第一UPD类型；基于每个SNP位点的第一UPD类型确定UPD区域；对该先证者的SNP位点进行AF过滤以确定该先证者的突变频率大于突变阈值的高频突变SNP位点；基于该先证者的高频突变SNP位点确定该先证者的ROH片段；以及基于该先证者的第一UPD类型和该先证者的ROH片段确定该先证者的第二UPD类型。

根据本发明的另一个方面，提供了一种计算设备。该计算设备包括：至少一个处理器；以及至少一个存储器，该至少一个存储器被耦合到该至少一个处理器并且存储用于由该至少一个处理器执行的指令，该指令当由该至少一个处理器执行时，使得该计算设备执行根据上述方法的步骤。

根据本发明的再一个方面，提供了一种计算机可读存储介质，其上存储有计算机程序代码，该计算机程序代码在被运行时执行如上所述的方法。

在一种实施例中，确定每个SNP位点的第一UPD类型包括：确定该先证者、该先证者的父样本和该先证者的母样本的每个SNP位点的BAF；基于该先证者、该先证者的父样本和该先证者的母样本的每个SNP位点的BAF分别确定该SNP位点的基因型；以及基于该先证者、该先证者的父样本和该先证者的母样本的每个SNP位点的基因型之间的关系确定该SNP位点的第一UPD类型。

在一种实施例中，确定该先证者、该先证者的父样本和该先证者的母样本的每个SNP位点的BAF包括：基于每个SNP位点的变异等位基因深度和总深度确定该SNP位点的BAF。

在一种实施例中，基于该先证者、该先证者的父样本和该先证者的母样本的每个SNP位点的BAF分别确定该SNP位点的基因型包括：确定该SNP位点的BAF是否小于第一阈值；如果确定该SNP位点的BAF小于该第一阈值，确定该SNP位点的基因型为野生型；如果确定该SNP位点的BAF大于或者等于该第一阈值，确定该SNP位点的BAF是否小于第二阈值，该第二阈值大于该第一阈值；如果确定该SNP位点的BAF小于该第二阈值，确定该SNP位点的基因型为杂合突变型；以及如果确定该SNP位点的BAF大于或者等于该第二阈值，确定该SNP位点的基因型为纯合突变型。

在一种实施例中，基于每个SNP位点的第一UPD类型确定UPD区域包括：确定在一个UPD滑动窗中，支持UPD的SNP位点与有效SNP位点的数量之比是否大于第一比率并且属于该第一UPD类型的连续有效SNP位点的数量是否大于预定值；如果确定在该UPD滑动窗中，支持UPD的SNP位点与有效SNP位点的数量之比大于该第一比率并且属于该第一UPD类型的连续有效SNP位点的数量大于该预定值，确定该UPD滑动窗的下一UPD滑动窗中，支持UPD的SNP位点与有效SNP位点的数量之比是否大于该第一比率并且属于该第一UPD类型的连续有效SNP位点的数量是否大于该预定值；以及如果确定该UPD滑动窗的下一UPD滑动窗中，支持UPD的SNP位点与有效SNP位点的数量之比大于该第一比率并且属于该第一UPD类型的连续有效SNP位点的数量大于该预定值，将该UPD滑动窗和该下一UPD滑动窗所覆盖的SNP位点设为一个UPD区域。

在一种实施例中，基于该先证者的高频突变SNP位点确定该先证者的ROH片段包括：获取基于相同探针捕获的多个样本所构建的群体标记库；确定该先证者的SNP位点；在该群体标记库中轮询该先证者的SNP位点，以确定该先证者的SNP位点的突变类型；以及基于该先证者的每个SNP位点的突变类型，以ROH滑动窗方式检测该先证者的所有SNP位点中的多个ROH片段。

在一种实施例中，获取基于相同探针捕获的多个样本所构建的群体标记库包括：获取相同探针捕获的多个样本；检测该多个样本中每个样本的每个SNP位点的突变频率；选择突变频率大于该突变阈值的SNP位点作为该高频突变SNP位点；检测该高频突变SNP位点中的杂合SNP位点的杂合突变频率和纯合SNP位点的纯合突变频率；以及将杂合突变频率大于杂合突变阈值的杂合SNP位点以及纯合突变频率大于纯合突变阈值的纯合SNP位点放入该群体标记库。

在一种实施例中，基于该先证者的每个SNP位点的突变类型，以ROH滑动窗方式检测该先证者的所有SNP位点中的多个ROH片段包括：确定一个ROH滑动窗中，突变类型为纯合突变型的SNP位点与该ROH滑动窗中的所有SNP位点的比值是否大于预定比值；如果确定该比值大于该预定比值，确定该ROH滑动窗的下一ROH滑动窗中基因型为纯合突变型的SNP位点与该下一ROH滑动窗中的所有SNP位点的比值是否大于该预定比值；以及如果该下一ROH滑动窗中基因型为纯合突变型的SNP位点与该下一ROH滑动窗中的所有SNP位点的比值大于该预定比值，将该ROH滑动窗和该下一ROH滑动窗所覆盖的SNP位点放入一个ROH片段。

在一种实施例中，基于该先证者的第一UPD类型和该先证者的ROH片段确定该先证者的第二UPD类型包括：确定该先证者的第一UPD类型中的母源UPD与该先证者的ROH片段的SNP位点之间的第一重合度；如果该第一重合度大于预定阈值，则确定该先证者的第二UPD类型为母源iUPD；确定该先证者的第一UPD类型中的父源UPD与该先证者的ROH片段之间的第二重合度；如果该第二重合度大于该预定阈值，则确定该先证者的第二UPD类型为父源iUPD；以及如果该第一重合度和该第二重合度均小于或等于该预定阈值，则确定该先证者的第二UPD类型为hUPD。

在一种实施例中，该方法还包括：确定该先证者的第二UPD类型的SNP位点占染色体的比例；如果该比例大于预定比例，则确定该先证者第三UPD类型为染色体级别UPD；如果该比例小于或等于该预定比例，确定该先证者的第三UPD类型为片段型UPD。

附图说明

通过参考下列附图所给出的本发明的具体实施方式的描述，将更好地理解本发明，并且本发明的其他目的、细节、特点和优点将变得更加显而易见。

图1示出了根据本发明的实施例的用于确定UPD类型的系统的示意图。

图2示出了根据本发明的实施例的用于确定UPD类型的方法的流程图。

图3A示出了根据本发明一种实施例的用于确定SNP位点的第一UPD类型的步骤的流程图。

图3B示出了根据本发明实施例的一个实例的用于确定第一阈值和第二阈值的BAF分布图。

图4示出了根据本发明一种实施例的用于确定UPD区域的步骤的流程图。

图5示出了根据本发明一种实施例的确定先证者的ROH片段的步骤的流程图。

图6示出了根据本发明一种实施例的用于确定先证者的第二UPD类型的步骤的流程图。

图7示出了适合实现本发明的实施例的计算设备的结构方框图。

具体实施方式

下面将参照附图更详细地描述本发明的优选实施例。虽然附图中显示了本发明的优选实施例，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。

在本文中使用的术语“包括”及其变形表示开放性包括，即“包括但不限于”。除非特别申明，术语“或”表示“和/或”。术语“基于”表示“至少部分地基于”。术语“一个实施例”和“一些实施例”表示“至少一个示例实施例”。术语“另一实施例”表示“至少一个另外的实施例”。术语“第一”、“第二”等等可以指代不同的或相同的对象。

如前所述，在当前的UPD检测中，能够检测的UPD类型非常有限，不能准确检测片段型UPD，并且检测精度仅限于大于10M碱基长度的UPD。

针对上述问题中的至少一个，本发明提供了一种确定UPD类型的方案，其中基于关于先证者、先证者的父样本和先证者的母样本的测序序列数据，分别确定先证者、先证者的父样本和先证者的母样本的同一基因序列中的各个SNP位点信息；基于该SNP位点信息，确定每个SNP位点的第一UPD类型；基于每个SNP位点的第一UPD类型确定UPD区域；对该先证者的SNP位点进行AF过滤以确定该先证者的突变频率大于突变阈值的高频突变SNP位点；基于该先证者的高频突变SNP位点确定该先证者的ROH片段；以及基于该先证者的第一UPD类型和该先证者的ROH片段确定该先证者的第二UPD类型。

在上述方案中，通过将ROH检测与传统的UPD检测相结合能够更加准确地检测UPD类型，尤其是片段型UPD。并且，在一些实施例中，通过滑动窗的方式来确定UPD区域和ROH区域，提高了UPD类型检测的精度，从而有助于进一步准确确定UPD是否引起印记效应。

图1示出了根据本发明的实施例的用于确定UPD类型的系统100的示意图。如图1中所示，系统100包括SNP位点确定单元110、计算设备120和数据库130。其中，计算设备120可以进一步包括第一UPD类型确定单元121、UPD区域确定单元122、高频突变确定单元123、ROH片段确定单元124和第二UPD类型确定单元125。此外，在一些实施例中，计算设备120还可以包括第三UPD类型确定单元126。

计算设备120可以包括至少一个处理器和与该至少一个处理器耦合的至少一个存储器，该存储器中存储有可由该至少一个处理器执行的指令，该指令在被该至少一个处理器执行时执行如下所述的方法200的至少一部分。第一UPD类型确定单元121、UPD区域确定单元122、高频突变确定单元123、ROH片段确定单元124、第二UPD类型确定单元125和第三UPD类型确定单元126可以分别实现为单独的硬件（如芯片），或者可以分别由上述指令的一部分来实现为软件形式。计算设备120的具体结构例如可以如下结合图7所述。

图2示出了根据本发明的实施例的用于确定UPD类型的方法200的流程图。方法200例如可以由图1中所示的系统100执行，更具体地，可以由图1中所示的系统100中的计算设备120执行。

如图2中所示，在步骤210，计算设备120基于关于先证者、先证者的父样本和先证者的母样本的测序序列数据，分别确定先证者、先证者的父样本和先证者的母样本的同一基因序列的各个SNP位点信息。

“先证者”是指在对某个遗传性状进行家系调查时，其家系中第一个被发现具有这种遗传性状（例如某种疾病）的那个人，根据先证者可以逐步追溯家系中其他成员的这种遗传性状出现的情况。在本文中，以先证者及其父亲（父样本）和母亲（母样本）作为研究样本。

SNP是single nucleotide polymorphism（单核苷酸多态性）的缩写，其是指在基因组水平上由单个核苷酸的变异所引起的DNA序列多态性。它是人类可遗传的变异中最常见的一种，占所有已知多态性的90%以上。SNP在人类基因组中广泛存在，平均每300个碱基对中就有1个，估计其总数可达300万个甚至更多。SNP是一种二态的标记，由单个碱基的转换或颠换所引起，也可由碱基的插入或缺失所致。SNP既可能在基因序列内，也可能在基因以外的非编码序列上。在本文中，以基因序列内的SNP为例进行描述。因此，确定先证者、先证者的父样本和先证者的母样本的同一基因序列的SNP位点即确定先证者的一个基因序列相对于其父亲（父样本）和母亲（母样本）的该同一基因序列发生SNP变异的碱基位置。在本文中，以一条染色体作为基因序列的一个实例进行描述，然而本领域技术人员可以理解，基因序列可以是一条染色体也可以是一条染色体的一个片段（即染色体片段）。

在步骤210中，确定各个SNP位点信息可以是由SNP位点确定单元110确定的，并且计算设备120从该SNP位点确定单元110接收所确定的各个SNP位点。该SNP位点确定单元110可以包含第三方的SNP检测软件，如VeritaTrekker^®或诸如GATK等其他开源软件，其可以根据接收的样本（如先证者、父样本或母样本）的测序序列数据确定该样本的SNP位点信息。该测序序列数据可以是通过本领域当前已知的或者未来开发的测序方法得到的，本发明对此不作限定。此外，虽然图1中将SNP位点确定单元110显示为位于计算设备120之外，但是SNP位点确定单元110也可以集成在计算设备120内。

在步骤220，计算设备120（如第一UPD类型确定单元121）基于步骤210确定的SNP位点信息确定每个SNP位点的第一UPD类型。

图3A示出了根据本发明一种实施例的用于确定SNP位点的第一UPD类型的步骤220的流程图。图3B示出了根据本发明实施例的一个实例的用于确定第一阈值和第二阈值的BAF分布图。

如图3A中所示，步骤220可以包括子步骤222，其中计算设备120可以确定先证者、先证者的父样本和先证者的母样本的每个SNP位点的变异等位基因频率（B-allelefrequency，BAF），BAF用来表示一个位点测序到变异碱基的读长比例，可以用于进一步评估位点的纯合状态。

在一种实施例中，一个SNP位点的BAF可以基于该SNP位点的变异等位基因深度和总深度来确定。例如，可以通过如下公式（1）确定一个SNP位点的BAF：

BAF_i=Dep_a/Dep_r+a （1）

其中，BAF_i表示SNP位点i的变异等位基因频率，Dep_a表示SNP位点i的变异等位碱基深度，Dep_r+a为表示SNP位点i的总深度，即覆盖该位点的所有读长数，包含覆盖参考碱基和变异等位碱基。可以理解，一个SNP位点的BAF是一个处于[0, 1]之间的值。

接下来，在子步骤224，计算设备120可以基于先证者、先证者的父样本和先证者的母样本的每个SNP位点的BAF分别确定该SNP位点的基因型。这里，SNP位点的基因型可以包括野生型（wild type，WT）、杂合突变型（heterozygous，HET）和纯合突变型（homozygosity，HOM）中的任一种。

具体地，在一种实施例中，子步骤224可以进一步包括：确定每个SNP位点i的BAF_i是否小于第一阈值，并且在确定SNP位点i的BAF_i小于该第一阈值时，确定该SNP位点i的基因型为野生型WT。这里，对于依赖于种系（gemline）变异的ROH/UPD分型来说，根据如图3B所示的实验结果，可以将该第一阈值设置为0.2左右，BAF_i低于该第一阈值则认为该位点i纯合未突变。其中，图3B中所示的BAF分布图是在500个样品的17267249个HET位点下检测到的。更一般地，在基因型为WT时，BAF接近或者等于0，在基因型为HOM时，BAF接近1，在基因型为HET时，BAF在0.5左右。

另一方面，如果确定SNP位点i的BAF_i大于或等于该第一阈值，则可以进一步确定该SNP位点i的BAF_i是否小于第二阈值，其中第二阈值大于该第一阈值。如果确定该SNP位点i的BAF_i小于该第二阈值，则可以确定该SNP位点i的基因型为杂合突变型HET。这里，根据如图3B所示的实验结果，可以将该第二阈值设置为0.8。

另一方面，如果确定SNP位点i的BAF_i大于或等于该第二阈值，则可以确定该SNP位点i的基因型为纯合突变型HOM。

也就是说，可以通过每个SNP位点的变异等位基因频率BAF处于哪个区间（小于第一阈值、第一阈值和第二阈值之间还是大于第二阈值）确定该SNP位点的基因型。

接下来，在子步骤226，计算设备120可以基于先证者、先证者的父样本和先证者的母样本的每个SNP位点的基因型之间的关系确定该SNP位点的第一UPD类型。

先证者、先证者的父样本和先证者的母样本的对应SNP位点的基因型分别有三种可能性（WT、HET和HOM），这三者一共有27种遗传组合。表1示出了先证者、父样本和母样本的对应SNP位点的基因型与UPD类型之间的对应关系。

从表1可以看出，基于先证者、父样本和母样本的SNP位点的基因型之间的关系，首先可以确定对应的SNP位点为无效SNP位点（no.eff）还是有效SNP位点（除了no.eff之外的UPD类型），对于有效SNP位点，所确定的UPD类型（这里为了与下面进行区分，称为第一UPD类型）包括五种，其中upd[mat]表示母源UPD，upd[pat]表示父源UPD，notupd表示不支持UPD，isoupd[mat]表示母源单亲同二体（母源iUPD），isoupd[pat]表示父源单亲同二体（父源iUPD）。

至此，对于有效SNP位点的UPD类型，能够准确区分出了iUPD的两种来源，即母源iUPD和父源iUPD，以及不支持UPD，但是对于母源UPD和父源UPD，仍无法准确确认其属于iUPD还是hUPD，也不能确认基因序列中的片段型UPD，因此这里的第一UPD类型可以认为是对SNP位点的初步分类结果。

继续图2，在步骤220之后，在步骤230，计算设备120（例如UPD区域确定单元122）可以基于每个SNP位点的第一UPD类型确定UPD区域。

图4示出了根据本发明一种实施例的用于确定UPD区域的步骤230的流程图。

如前面子步骤226中结合表1所述，基于先证者、父样本和母样本的SNP位点的基因型之间的关系可以确定对应的SNP位点是无效SNP位点（no.eff）还是有效SNP位点（除了no.eff之外的UPD类型），并且进一步地，第一UPD类型包括母源UPD、父源UPD、母源iUPD、父源iUPD和不支持UPD。这里，相对于不支持UPD的SNP位点，将除了不支持UPD的有效SNP位点之外的有效SNP位点称为支持UPD的SNP位点。

如图4中所示，在步骤230的子步骤232，计算设备120确定在一个UPD滑动窗中，支持UPD的SNP位点与有效SNP位点的数量之比是否大于第一比率并且属于第一UPD类型的连续有效SNP位点的数量是否大于预定值。这里，UPD滑动窗是一个具有特定长度的窗口，也可以称为UPD检测步长。在本发明的一些实施例中，UPD滑动窗的大小被设置为2M碱基长度。

该第一比率用于筛选掉UPD滑动窗中不支持UPD的SNP位点过多的那些UPD滑动窗，这种滑动窗中由于不支持UPD的SNP位点的占比过大而对于UPD区域划分没有意义。在一种实施例中，该第一比率可以设置为[0.8, 1)之间的一个值。优选地，该第一比率可以设置为0.9。

该预定值用于作为UPD滑动窗中对于某一UPD类型（母源UPD、父源UPD、母源iUPD或者父源iUPD）的连续有效SNP位点的判断阈值。当某一UPD类型的连续有效SNP位点的数量过少时，认为该UPD滑动窗的这种UPD类型的SNP位点不能构成一个UPD区域。在一种实施例中，可以根据滑动窗口的大小设置该预定值。例如，考虑到检测灵敏度，可以将该预定值设置为5。

如果在该UPD滑动窗中，支持UPD的SNP位点与有效SNP位点的数量之比大于该第一比率并且属于该第一UPD类型的连续有效SNP位点的数量大于该预定值，在子步骤234中，继续确定该UPD滑动窗的下一UPD滑动窗是否满足上述条件，即，该下一UPD滑动窗中支持UPD的SNP位点与有效SNP位点的数量之比是否大于该第一比率并且属于该第一UPD类型的连续有效SNP位点的数量是否大于该预定值。

如果确定该下一UPD滑动窗也满足上述条件，在子步骤236中，则将该UPD滑动窗和该下一UPD滑动窗所覆盖的SNP位点设为一个UPD区域。

另一方面，如果确定该下一UPD滑动窗不满足上述条件，在子步骤237中，可以将该UPD滑动窗所覆盖的SNP位点单独设为一个UPD区域。

子步骤234和236可以重复执行，直至遇到一个不满足上述条件的UPD滑动窗为止，这样，属于同一UPD类型的所有相邻UPD滑动窗所覆盖的SNP位点都被放入一个UPD区域。尤其是，在对一个基因序列重复执行上述操作之后，可能产生多个UPD区域，其中一些UPD区域中可能具有确定的细分UPD类型，如父源iUPD或母源iUPD，其一方面指明了来源亲本（父源或母源），另一方面指明了来自同一亲本的同一条染色体（iUPD）。而另一些UPD区域的细分UPD类型并不确定，如父源UPD或母源UPD，其仅仅指明了来源亲本（父源或母源），而不能指明是来自同一亲本的同一条染色体（iUPD）还是来自同一亲本的一对同源染色体（hUPD）。因此，可以根据该UPD的位点纯合度来进一步分型，一般的纯合度较高，即存在ROH时，分型为iUPD，其余为hUPD，如下所详述。

继续图2，在方法200的另一条支路上，在步骤240，计算设备120（例如高频突变确定单元123）对先证者的SNP位点进行AF过滤以确定先证者的、突变频率大于突变阈值的高频突变SNP位点。

具体地，在一种实施例中，在步骤240中，首先计算先证者的SNP位点的突变频率。例如，可以通过测量SNP位点的等位基因频率（allele frequency，AF）来确定SNP位点的突变频率。测量SNP位点的AF的方法可以使用现有技术中已知的或者未来开发的各种方法，而不会限制本发明的范围。

在一种实施例中，可以将该给定突变阈值设置为10%。进一步地，可以分别检测杂合SNP位点的杂合突变频率（如AF_het）和纯合SNP位点的纯合突变频率（如AF_hom），并且确定杂合SNP位点的杂合突变频率（如AF_het）和纯合SNP位点的纯合突变频率（如AF_hom）是否分别大于相应的杂合突变阈值和纯合突变阈值，从而确定高频突变杂合SNP位点和高频突变纯合SNP位点。在一种实施例中，可以将杂合突变频率和纯合突变频率分别设置为5%。

在步骤250，计算设备120（例如ROH片段确定单元124）可以基于先证者的高频突变SNP位点确定先证者的ROH片段。

图5示出了根据本发明一种实施例的确定先证者的ROH片段的步骤250的流程图。

如图5中所示，步骤250可以包括子步骤252，其中计算设备120可以获取基于相同探针捕获的多个样本所构建的群体标记（Marker）库。

更具体地，在子步骤252中，计算设备120可以获取相同探针捕获的多个样本。这里，样本的数量可以是以“千”为量级，例如一千个。使用相同探针捕获保证了构建的群体标记库中的样本的检测标准一致。

然后，计算设备120可以检测多个样本中每个样本的每个SNP位点的突变频率，并且选择突变频率大于给定突变阈值的SNP位点作为高频突变SNP位点。这里，确定高频突变SNP位点的方式与上述步骤240中类似，并且这里的给定突变阈值也与上述步骤240中相同。进一步地，可以分别检测高频突变SNP位点中的杂合SNP位点的杂合突变频率（如AF_het）和纯合SNP位点的纯合突变频率（如AF_hom），并且将杂合突变频率大于杂合突变阈值的杂合SNP位点（也称为高频突变杂合SNP位点）以及纯合突变频率大于纯合突变阈值的纯合SNP位点（也称为高频突变纯合SNP位点）放入群体标记库中。也就是说，群体标记库中的SNP位点需要满足突变频率大于上述突变阈值，杂合突变频率大于上述杂合突变阈值并且纯合突变频率大于上述纯合突变阈值。

群体标记库可以是在方法200执行过程中产生的，也可以是预先产生并且存储在数据库130中。在后者的情况下，在子步骤252中，计算设备120可以从数据库130获取该群体标记库。注意，虽然图1中将数据库130显示为独立于计算设备120，但是在其他实现中，数据库130也可以与计算设备120集成在一起。

继续图5，在步骤250的另一条支路上，在子步骤254确定先证者的SNP位点。这里，确定先证者的SNP位点的方法与上述步骤210中确定SNP位点的方法相同，因此不再赘述。

然后，在子步骤256，在子步骤252获取的群体标记库中轮询子步骤254确定的先证者的SNP位点，以确定该先证者的SNP位点的突变类型，即，确定该先证者的SNP位点是否是高频突变SNP位点，以及是高频突变杂合SNP位点（即杂合突变型HET）还是高频突变纯合SNP位点（即纯合突变型HET）。

接下来，在子步骤258中，计算设备120基于先证者的每个SNP位点的突变类型，以ROH滑动窗方式检测该先证者的所有SNP位点中的多个ROH片段。

更具体地，在一些实施例中，子步骤258可以进一步包括：计算设备120确定一个ROH滑动窗中，突变类型为纯合突变型的SNP位点与该ROH滑动窗中的所有SNP位点的比值是否大于预定比值。这里，ROH滑动窗是一个具有特定长度的窗口，也可以称为ROH检测步长。在本发明的一些实施例中，ROH滑动窗的大小被设置为与UPD滑动窗的大小相同，例如被设置为2M碱基长度，每个滑动窗内约包含50个标记点。此外，该预定比值用于筛选出大部分SNP位点是纯合突变型的ROH滑动窗，因此该预定比值是一个接近1的值，例如可以是[0.8,1]之间的一个值。优选地，该预定比值例如是0.90。

如果确定该比值大于该预定比值，计算设备120继续确定该ROH滑动窗的下一ROH滑动窗中突变类型为纯合突变型的SNP位点与所述下一ROH滑动窗中的所有SNP位点的比值是否大于该预定比值。

如果该下一ROH滑动窗中突变类型为纯合突变型的SNP位点与该下一ROH滑动窗中的所有SNP位点的比值也大于该预定比值，将该ROH滑动窗和该下一ROH滑动窗所覆盖的SNP位点设为一个ROH片段。

另一方面，如果确定该下一ROH滑动窗中突变类型为纯合突变型的SNP位点与该下一ROH滑动窗中的所有SNP位点的比值不大于该预定比值，可以将该ROH滑动窗所覆盖的SNP位点单独设为一个ROH片段。

通过这种方式，可以基于先证者的高频突变SNP位点，并且借助于该群体标记库共同确定该先证者的ROH片段。

上述过程可以重复执行，直至遇到一个不满足上述预定比值的ROH滑动窗为止，这样，突变类型为纯合突变型的所有相邻ROH滑动窗所覆盖的SNP位点都被设置在一个ROH片段中。

在对先证者的基因序列重复执行上述操作之后，可能产生多个ROH片段，这些ROH片段分别包含突变类型为纯合突变型的SNP位点。

继续图2，在步骤260中，计算设备120（例如第二UPD类型确定单元125）可以基于步骤220确定的先证者的第一UPD类型和步骤250确定的先证者的ROH片段确定先证者的第二UPD类型。

如前面步骤220中所述，有效SNP位点的第一UPD类型可以是母源UPD、父源UPD、母源iUPD、父源iUPD和不支持UPD。其中，对于母源UPD和父源UPD的SNP位点，仅通过步骤220无法确认其属于iUPD还是hUPD。另一方面，在步骤250中，ROH片段中包含的SNP位点基本上都是纯合突变型的SNP位点。因此，在步骤260中，可以通过对二者进行整合来进一步将步骤220中确定的母源UPD和父源UPD的SNP位点确认为iUPD还是hUPD。

图6示出了根据本发明一种实施例的用于确定先证者的第二UPD类型的步骤260的流程图。

如图6中所示，步骤260包括子步骤262，其中计算设备120确定先证者的第一UPD类型中的母源UPD与先证者的ROH片段的SNP位点之间的第一重合度。

在子步骤264，如果计算设备120确定该第一重合度大于预定阈值，则确定该先证者的第二UPD类型为母源iUPD。

在子步骤266，计算设备120可以确定先证者的第一UPD类型中的父源UPD与该先证者的ROH片段之间的第二重合度。

在子步骤268，如果计算设备120确定该第二重合度大于该预定阈值，则确定该先证者的第二UPD类型为父源iUPD。

这里，理论上iUPD和ROH是重叠的，但限于捕获的样本中，边界检测存在难度，为了提升UPD分型率，该预定阈值可以设置为0.8，用于判断是父源UPD和母源UPD与纯合区域ROH的重合程度。

另一方面，如果计算设备120确定第一重合度和第二重合度都小于或等于该预定阈值，则可以确定该先证者的第二UPD类型为hUPD。

通过这种方式，可以将步骤220确定的父源UPD和母源UPD进一步确定为是iUPD还是hUPD，从而进一步准确地确定了UPD类型。

然而，这种方式仍然无法识别出片段型UPD（segUPD）。

为此，在一些实施例中，方法200还可以包括步骤270，其中计算设备120（例如通过第三UPD类型确定单元126）可以根据步骤220确定的第一UPD类型和步骤260确定的第二UPD类型，进一步确定先证者的第三UPD类型，即确定是否是片段型UPD。

具体地，在一些实施例中，步骤270可以包括：计算设备120确定先证者的第二UPD类型的SNP位点占染色体的比例。如果该比例大于预定比例，则可以确定先证者的第三UPD类型为染色体级别UPD。这里，该预定比例是一个接近1的值，例如可以在[0.8, 1)之间，优选地为0.9。也就是说，如果第二UPD类型的SNP位点（即某种UPD区域）占据染色体（这里假设该基因序列是一条染色体）的绝大部分，则可以将该先证者的UPD类型确定为染色体级别UPD。

否则，将该先证者的第三UPD类型确定为片段型UPD（segUPD）。

通过上述步骤210至270，可以精确地确定先证者的UPD类型是染色体级别的iUPD、hUPD还是segUPD。

图7示出了适合实现本发明的实施例的计算设备700的结构方框图。计算设备700例如可以是如上所述的计算设备120。

如图7中所示，计算设备700可以包括一个或多个中央处理单元（CPU）710（图中仅示意性地示出了一个），其可以根据存储在只读存储器（ROM）720中的计算机程序指令或者从存储单元780加载到随机访问存储器（RAM）730中的计算机程序指令，来执行各种适当的动作和处理。在RAM 730中，还可存储计算设备700操作所需的各种程序和数据。CPU 710、ROM 720以及RAM 730通过总线740彼此相连。输入/输出（I/O）接口750也连接至总线740。

计算设备700中的多个部件连接至I/O接口750，包括：输入单元760，例如键盘、鼠标等；输出单元770，例如各种类型的显示器、扬声器等；存储单元780，例如磁盘、光盘等；以及通信单元790，例如网卡、调制解调器、无线通信收发机等。通信单元790允许计算设备700通过诸如因特网的计算机网络和/或各种电信网络与其他设备交换信息/数据。

上文所描述的方法200例如可由计算设备700的CPU 710执行。例如，在一些实施例中，方法200可被实现为计算机软件程序，其被有形地包括于机器可读介质，例如存储单元780。在一些实施例中，计算机程序的部分或者全部可以经由ROM 720和/或通信单元790而被载入和/或安装到计算设备700上。当计算机程序被加载到RAM 730并由CPU 710执行时，可以执行上文描述的方法200的一个或多个操作。此外，通信单元790可以支持有线或无线通信功能。

本领域技术人员可以理解，图7所示的计算设备700仅是示意性的。在一些实施例中，计算设备700可以包含比图7所示更多或更少的部件。

芯片验证的结果显示，本发明的方案比相同功能的H3M2等软件假阳性率更低，与芯片结果的一致性更高，可以达到90%。并且，利用家系位点分析MIE来鉴定segUPD的方法使得UPD的检测精度可以达到2M碱基长度，这对于是否引起印记效应的解读将更加精准。

此外，通过捕获样本来构建群体标记库的方法来鉴定ROH，将其与UPD区域相结合可以更加精确的区分SNP位点的纯合状态的原因，是捕获问题还是本身纯合。此外，利用滑动窗口的方法检测UPD能够提升UPD的检测精度。

以上结合附图对根据本发明的基于ROH确定UPD类型的方法200以及可用于执行方法200的计算设备700进行了描述。然而本领域技术人员可以理解，方法200的步骤或子步骤的执行并不局限于图中所示和以上所述的顺序，而是可以以任何其他合理的顺序来执行。例如，图2所示的方法200的步骤220和230可以与步骤240和250并行执行或以其他顺序执行。此外，计算设备700也不必须包括图7中所示的所有组件，其可以仅仅包括执行本发明中所述的功能所必须的其中一些组件，并且这些组件的连接方式也不局限于图中所示的形式。

本发明可以是方法、装置、系统和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于执行本发明的各个方面的计算机可读程序指令。

在一个或多个示例性设计中，可以用硬件、软件、固件或它们的任意组合来实现本发明所述的功能。例如，如果用软件来实现，则可以将所述功能作为一个或多个指令或代码存储在计算机可读介质上，或者作为计算机可读介质上的一个或多个指令或代码来传输。

本文公开的装置的各个单元可以使用分立硬件组件来实现，也可以集成地实现在一个硬件组件，如处理器上。例如，可以用通用处理器、数字信号处理器（DSP）、专用集成电路（ASIC）、现场可编程门阵列（FPGA）或其它可编程逻辑器件、分立门或者晶体管逻辑、分立硬件组件或用于执行本文所述的功能的任意组合来实现或执行结合本发明所描述的各种示例性的逻辑块、模块和电路。

本领域普通技术人员还应当理解，结合本发明的实施例描述的各种示例性的逻辑块、模块、电路和算法步骤可以实现成电子硬件、计算机软件或二者的组合。

本发明的以上描述用于使本领域的任何普通技术人员能够实现或使用本发明。对于本领域普通技术人员来说，本发明的各种修改都是显而易见的，并且本文定义的一般性原理也可以在不脱离本发明的精神和保护范围的情况下应用于其它变形。因此，本发明并不限于本文所述的实例和设计，而是与本文公开的原理和新颖性特性的最广范围相一致。

Claims (12)

1.一种用于确定单亲二倍体类型的方法，包括：

基于关于先证者、先证者的父样本和先证者的母样本的测序序列数据，分别确定先证者、先证者的父样本和先证者的母样本的同一基因序列中的各个单核苷酸多态性位点信息；

基于所述单核苷酸多态性位点信息，确定每个单核苷酸多态性位点的第一单亲二倍体类型，其中所述第一单亲二倍体类型包括父源单亲二倍体和母源单亲二倍体；

基于每个单核苷酸多态性位点的第一单亲二倍体类型确定单亲二倍体区域；

对所述先证者的单核苷酸多态性位点进行AF过滤以确定所述先证者的、突变频率大于突变阈值的高频突变单核苷酸多态性位点；

基于所述先证者的高频突变单核苷酸多态性位点确定所述先证者的纯合区域片段；以及

基于所述先证者的第一单亲二倍体类型和所述先证者的纯合区域片段确定所述先证者的第二单亲二倍体类型，其中所述第二单亲二倍体类型包括母源单亲同二体、父源单亲同二体和单亲异二体。

2.如权利要求1所述的方法，其中确定每个单核苷酸多态性位点的第一单亲二倍体类型包括：

确定所述先证者、所述先证者的父样本和所述先证者的母样本的每个单核苷酸多态性位点的变异等位基因频率；

基于所述先证者、所述先证者的父样本和所述先证者的母样本的每个单核苷酸多态性位点的变异等位基因频率分别确定所述单核苷酸多态性位点的基因型；以及

基于所述先证者、所述先证者的父样本和所述先证者的母样本的每个单核苷酸多态性位点的基因型之间的关系确定所述单核苷酸多态性位点的第一单亲二倍体类型。

3.如权利要求2所述的方法，其中确定所述先证者、所述先证者的父样本和所述先证者的母样本的每个单核苷酸多态性位点的变异等位基因频率包括：

基于每个单核苷酸多态性位点的变异等位基因深度和总深度确定所述单核苷酸多态性位点的变异等位基因频率。

4.如权利要求2所述的方法，其中基于所述先证者、所述先证者的父样本和所述先证者的母样本的每个单核苷酸多态性位点的变异等位基因频率分别确定所述单核苷酸多态性位点的基因型包括：

确定所述单核苷酸多态性位点的变异等位基因频率是否小于第一阈值；

如果确定所述单核苷酸多态性位点的变异等位基因频率小于所述第一阈值，确定所述单核苷酸多态性位点的基因型为野生型；

如果确定所述单核苷酸多态性位点的变异等位基因频率大于或者等于所述第一阈值，确定所述单核苷酸多态性位点的变异等位基因频率是否小于第二阈值，所述第二阈值大于所述第一阈值；

如果确定所述单核苷酸多态性位点的变异等位基因频率小于所述第二阈值，确定所述单核苷酸多态性位点的基因型为杂合突变型；以及

如果确定所述单核苷酸多态性位点的变异等位基因频率大于或者等于所述第二阈值，确定所述单核苷酸多态性位点的基因型为纯合突变型。

5.如权利要求1所述的方法，其中基于每个单核苷酸多态性位点的第一单亲二倍体类型确定单亲二倍体区域包括：

确定在一个单亲二倍体滑动窗中，支持单亲二倍体的单核苷酸多态性位点与有效单核苷酸多态性位点的数量之比是否大于第一比率并且属于所述第一单亲二倍体类型的连续有效单核苷酸多态性位点的数量是否大于预定值；

如果确定在所述单亲二倍体滑动窗中，支持单亲二倍体的单核苷酸多态性位点与有效单核苷酸多态性位点的数量之比大于所述第一比率并且属于所述第一单亲二倍体类型的连续有效单核苷酸多态性位点的数量大于所述预定值，确定所述单亲二倍体滑动窗的下一单亲二倍体滑动窗中，支持单亲二倍体的单核苷酸多态性位点与有效单核苷酸多态性位点的数量之比是否大于所述第一比率并且属于所述第一单亲二倍体类型的连续有效单核苷酸多态性位点的数量是否大于所述预定值；以及

如果确定所述单亲二倍体滑动窗的下一单亲二倍体滑动窗中，支持单亲二倍体的单核苷酸多态性位点与有效单核苷酸多态性位点的数量之比大于所述第一比率并且属于所述第一单亲二倍体类型的连续有效单核苷酸多态性位点的数量大于所述预定值，将所述单亲二倍体滑动窗和所述下一单亲二倍体滑动窗所覆盖的单核苷酸多态性位点设为一个单亲二倍体区域。

6.如权利要求1所述的方法，其中基于所述先证者的高频突变单核苷酸多态性位点确定所述先证者的纯合区域片段包括：

获取基于相同探针捕获的多个样本所构建的群体标记库；

确定所述先证者的单核苷酸多态性位点；

在所述群体标记库中轮询所述先证者的单核苷酸多态性位点，以确定所述先证者的单核苷酸多态性位点的突变类型；以及

基于所述先证者的每个单核苷酸多态性位点的突变类型，以纯合区域滑动窗方式检测所述先证者的所有单核苷酸多态性位点中的多个纯合区域片段。

7.如权利要求6所述的方法，其中获取基于相同探针捕获的多个样本所构建的群体标记库包括：

获取相同探针捕获的多个样本；

检测所述多个样本中每个样本的每个单核苷酸多态性位点的突变频率；

选择突变频率大于所述突变阈值的单核苷酸多态性位点作为所述高频突变单核苷酸多态性位点；

检测所述高频突变单核苷酸多态性位点中的杂合单核苷酸多态性位点的杂合突变频率和纯合单核苷酸多态性位点的纯合突变频率；以及

将杂合突变频率大于杂合突变阈值的杂合单核苷酸多态性位点以及纯合突变频率大于纯合突变阈值的纯合单核苷酸多态性位点放入所述群体标记库。

8.如权利要求6所述的方法，其中基于所述先证者的每个单核苷酸多态性位点的突变类型，以纯合区域滑动窗方式检测所述先证者的所有单核苷酸多态性位点中的多个纯合区域片段包括：

确定一个纯合区域滑动窗中，突变类型为纯合突变型的单核苷酸多态性位点与所述纯合区域滑动窗中的所有单核苷酸多态性位点的比值是否大于预定比值；

如果确定所述比值大于所述预定比值，确定所述纯合区域滑动窗的下一纯合区域滑动窗中基因型为纯合突变型的单核苷酸多态性位点与所述下一纯合区域滑动窗中的所有单核苷酸多态性位点的比值是否大于所述预定比值；以及

如果所述下一纯合区域滑动窗中基因型为纯合突变型的单核苷酸多态性位点与所述下一纯合区域滑动窗中的所有单核苷酸多态性位点的比值大于所述预定比值，将所述纯合区域滑动窗和所述下一纯合区域滑动窗所覆盖的单核苷酸多态性位点放入一个纯合区域片段。

9.如权利要求1所述的方法，其中基于所述先证者的第一单亲二倍体类型和所述先证者的纯合区域片段确定所述先证者的第二单亲二倍体类型包括：

确定所述先证者的第一单亲二倍体类型中的母源单亲二倍体与所述先证者的纯合区域片段的单核苷酸多态性位点之间的第一重合度；

如果所述第一重合度大于预定阈值，则确定所述先证者的第二单亲二倍体类型为母源单亲同二体；

确定所述先证者的第一单亲二倍体类型中的父源单亲二倍体与所述先证者的纯合区域片段之间的第二重合度；

如果所述第二重合度大于所述预定阈值，则确定所述先证者的第二单亲二倍体类型为父源单亲同二体；以及

如果所述第一重合度和所述第二重合度均小于或等于所述预定阈值，则确定所述先证者的第二单亲二倍体类型为单亲异二体。

10.如权利要求1所述的方法，还包括：

确定所述先证者的第二单亲二倍体类型的单核苷酸多态性位点占染色体的比例；

如果所述比例大于预定比例，则确定所述先证者第三单亲二倍体类型为染色体级别单亲二倍体；以及

如果所述比例小于或等于所述预定比例，确定所述先证者的第三单亲二倍体类型为片段型单亲二倍体。

11.一种计算设备，包括：

至少一个处理器；以及

至少一个存储器，所述至少一个存储器被耦合到所述至少一个处理器并且存储用于由所述至少一个处理器执行的指令，所述指令当由所述至少一个处理器执行时，使得所述计算设备执行根据权利要求1至10中任一项所述的方法的步骤。

12.一种计算机可读存储介质，其上存储有计算机程序代码，所述计算机程序代码在被运行时执行如权利要求1至10中任一项所述的方法。

Images