CN113113081B - 基于CNV-seq测序数据检测多倍体和基因组纯合区域ROH的*** - Google Patents

Abstract

本发明公开了基于CNV‑seq测序数据检测多倍体和基因组纯合区域ROH的***。发明人研究发现，不同类型的样本或基因组区域对应的基因型的种类和占比是不同的，在在极低平均测序深度下，虽然每个SNP位点的基因型不可知，但是同一测序深度下不同类型的样本/基因组区域观察到的杂合信号P(1,1)的概率不同。通过合并计算N个SNP杂合信息、计算SNP混合杂合度指标以判断样本的倍型信息和是否存在基因组纯合区域ROH。打破了现有低平均测序深度数据不能用于判断ROH和多倍体的局限，大幅减少了ROH和多倍体检测的成本。

Description

基于CNV-seq测序数据检测多倍体和基因组纯合区域ROH的 ***

技术领域

本发明涉及产前诊断细胞遗传学检测技术，特别是基于低测序深度全基因组测序数据的胎儿多倍体、基因组纯合区域(ROH，region of homozygosity)的检测***和建模方法。

背景技术

染色体异常是导致自然流产、胎儿多发畸形和出生缺陷的重要临床因素。染色体异常测常见类型包括：染色体平衡易位(染色体拷贝数未改变但发生倒位、易位、成环等)、染色体拷贝数变异(copy number variation；简称CNV)、异倍体(三倍体、四倍体等)、基因组纯合区域(absence of heterozygosity；简称ROH)等。约20％的可检测妊娠以流产终结，其中染色体异常是导致孕早期流产的最常见因素。超过50％的孕早期流产的样本存在不同类型的染色体异常[1-3]，其中染色体非整倍体、多倍体、染色体微缺失/微重复、基因组纯合区域 ROH分别约占70％～80％、12％～15％、5％和1～2％。在超声提示结构异常的胎儿中，有5％～ 10％的样本存在明确致病或者可能致病的染色体异常[4]。因此，染色体异常检测是产前诊断的重要手段，有助于查明流产原因和降低出生缺陷率。

正常的人类体细胞为二倍体细胞，拥有两份基因组(两套染色体)，一份来自父亲，另一份来自母亲。胚胎发生三倍体异常，通常是由于亲代的同源染色体在减数***过程中不分离或双精授精，导致胚胎比正常的二倍体胚胎多出一套染色体。多倍体胚胎通常会在孕早期流产，占孕早期流产中染色体异常样本的10％左右。

基因组纯合区域ROH指的是基因组区域中呈现的等位基因杂合性缺失现象。对于大部分的二倍体细胞如人类体细胞，拥有两份基因组，一份来自父亲，另一份来自母亲。在基因组某个SNP位点，如果来自父本和母本的碱基不同时，则称为该SNP位点杂合(heterozygous)。如果因为某种机制(如缺失、减数***错误、近亲等)导致在该区域中只有一类基因组(来自父本或母本)，则该区域为基因组纯合区域ROH。ROH中最有代表性的一类为单亲二体 (uniparental disomy；UPD)，指来自父母一方的染色体片段被另一方的同源部分取代，或一个个体的两条同源染色体都来自于同一个亲本。如果整条染色体发生UPD，则称为单亲二倍体。基因组上存在的一系列带有亲源烙印而出现不同生物学功能的基因，则被称为印记基因。在基因组的大部分区域，发生ROH不会由于印记基因缺失导致临床表型，但是存在印记基因的区域例如染色体6、7、11、14、15、20号染色体发生ROH则可能会导致严重的表型，其中最具代表性的区域为15q11-13，该区域upd(15)mat和upd(15)pat分别对应父源印记基因缺失Prader—Willi综合征(PWS综合征)和母源印记基因缺失Angelman综合征(AS综合征)。此外，发生ROH会增加该区域发生隐性遗传基因纯合变异的可能性，例如女性的两条X染色体如发现了ROH，则可能会导致该女性X连锁的隐性遗传病发生。

目前临床上应用最为广泛的全基因组范围内检测染色体异常检测的方法为核型分析、染色体微阵列分析(CMA，chromosomal microarray analysis)以及基于低通量全基因组测序技术的染色体拷贝数异常检测(copy number variation sequencing；简称CNV-seq)方法。这三种方法分别各有特点和优劣[3-5]，如表1所示。

表1:全基因组水平染色体异常检测方法对比

(1)核型分析目前仍然是产前一线诊断方法，其独特的优势在于可以直接观察到染色体的形态，可识别染色体平衡易位，如倒位、成环等，但是核型分析无法检测<5Mb的拷贝数异常，同时实验操作繁琐容易失败。

(2)染色体微阵列分析(CMA，chromosomal microarray analysis)。根据探针设计的不同，CMA技术可分为两大类：基于微阵列的比较基因组杂交(array-based comparativegenomic hybridization，aCGH)技术和单核苷酸多态性微阵列(single nucleotidepolymorphism array，SNP array)技术。aCGH技术只包括非多态性探针因此只能检测染色体拷贝数异常； SNP array在非多态性探针的基础上增加了SNP探针，因此除了染色体拷贝数异常外，还能够检测出大多数的ROH和三倍体。由于目前市场上主流的CMA方法都是基于SNP array技术，因此本文所述CMA方法如无特殊说明均基于SNP array。染色体微阵列CMA针对SNP 位点进行了扫描，可识别出每个SNP位点的基因型，如AA、AB和BB。CMA方法利用SNP提供的基因型信息可以发现多倍体异常以及基因组中纯合子区域ROH[2,6]。但是，CMA方法的主要缺陷在于拷贝数异常检测的准确性比CNV-seq低，容易出现拷贝数异常的假阳性结果，无法检测低于30％嵌合比例的拷贝数异常，同时CMA的检测成本也被CNV-seq高1～3倍。

(3)CNV-seq，基于高通量测序技术的全基因组测序方法来检测基因组拷贝数变异的染色体分析技术。CNV-seq在极低测序深度(0.1～1X)的基础下即可准确进行全基因组水平的染色体拷贝数异常检测，可检测长度低至100kb、嵌合比例低至10％的染色体拷贝数异常。 CNV-seq对于拷贝数异常的检测的准确性比CMA更高[4,7]，并且由于CNV-seq的检测成本显著低于CMA方法，因此CNV-seq在产前筛查领域也具有重要的地位。但是，现有的CNV-seq 方法只能检测拷贝数异常，无法检测多倍体异常和基因组纯合区域ROH。这一缺陷在针对流产物检测的时候尤为突出，因为孕早期流产的染色体异常样本中，12～15％存在多倍体、2％～ 4％存在ROH异常，这些样本使用CNV-seq的检测方法会直接导致漏检。

人类基因组是二倍体，一套DNA来自父亲，另外一套DNA来自母亲，所有双等位基因SNP位点(biallele SNP)位点只有三种基因型AA、AB、BB三种可能。但是胚胎，尤其是早期的胚胎中，有一部分是存在多倍体异常(常见三倍体和四倍体)和ROH纯合状态异常， ROH区域只有AA、BB两种基因型；三倍体样本则有AAA、AAB、ABB、BBB四种基因型；四倍体样本则有AAAA、AAAB、AABB、ABBB、BBBB四种基因型。

染色体微阵列分析CMA或针对SNP的高深度测序法例如捕获测序可以准确地判断每个 SNP位点的基因型，进而简单地统计各种基因型存在的占比即可判断样本是否存在多倍体和 ROH异常。这种操作虽然结果可靠，但是成本相对过高。

全基因组测序的成本随着测序深度的增加而急剧增加，如何使用低测序深度的数据满足检测需要，是一项极具挑战的工作。因为在极低测序深度下的全基因组高通量测序数据中(如 CNV-seq数据)，绝大部分SNP位点的深度都是1X，少数SNP位点的测序深度为2X或3X (如图1所示)，在这种情况下，无法对每个SNP位点的真实基因型做出判断。因此，CNV-seq 方法一直局限于拷贝数异常检测，行业认为CNV-seq方法在检测拷贝数异常方面的准确性和分辨率优于染色体微阵列CMA技术，但是无法判断多倍体和ROH。2019年，中国医师协会医学遗传分会发表《低深度全基因组测序技术在产前诊断中的应用专家共识》(参见中华医学会医学遗传学分会临床遗传学组,中国医师协会医学遗传医师分会遗传病产前诊断专业委员会,中华预防医学会出生缺陷预防与控制专业委员会遗传病防控学组.低深度全基因组测序技术在产前诊断中的应用专家共识[J].中华医学遗传学杂志,2019,36(4):293-296.)，认为基于低深度全基因组测序技术的CNV-seq无法检测多倍体和ROH异常，建议临床考虑结合 STR方法来解决多倍体检测的问题。

参考文献：

1.Smits,M.A.J.,et al.,Cytogenetic testing of pregnancy loss tissue:ameta-analysis. Reprod Biomed Online,2020.

2.Sahoo,T.,et al.,Comprehensive genetic analysis of pregnancy loss bychromosomal microarrays:outcomes,benefits,and challenges.Genet Med,2017.19(1):p.83-89.

3.Shearer,B.M.,et al.,Reflex fluorescent in situ hybridizationtesting for unsuccessful product of conception cultures:a retrospectiveanalysis of 5555samples attempted by conventional cytogenetics andfluorescent in situ hybridization.Genet Med,2011.13(6):p.545-52.

4.Zhou,B.,et al.,Whole-genome sequencing analysis of CNV using low-coverage and paired-end strategies is efficient and outperforms array-basedCNV analysis.J Med Genet,2018. 55(11):p.735-743.

5.Levy,B.,et al.,Genomic imbalance in products of conception:single-nucleotide polymorphism chromosomal microarray analysis.Obstet Gynecol,2014.124(2Pt 1):p.202-9.

6.Wang,J.C.,et al.,Regions of homozygosity identified byoligonucleotide SNP arrays: evaluating the incidence and clinical utility.EurJ Hum Genet,2015.23(5):p.663-71.

7.Wang,J.,et al.,Prospective chromosome analysis of 3429amniocentesissamples in China using copy number variation sequencing.Am J Obstet Gynecol,2018.219(3):p.287e1-287 e18.

8.Olshen,A.B.,et al.,Circular binary segmentation for the analysis ofarray-based DNA copy number data.Biostatistics,2004.5(4):p.557-72.。

发明内容

本发明的目的在于克服现有技术的至少一种不足，提供一种基因组纯合区域ROH的检测***。

本发明开创性地提出SNP混合杂合度指标SMHS(snp_mixed_heterozygosity_score)，表明在极低测序深度下(0.12X覆盖深度)，虽然每个SNP位点的基因型不可知，但是通过合并计算N个SNP杂合信息，计算snp_mixed_heterozygosity_score可以判断样本的倍型信息(二倍体、多倍体)和是否存在基因组纯合区域ROH。

本发明所采取的技术方案是：

本发明的第一个方面，提供：

一种基于低深度测序数据的基因组纯合区域ROH的检测***，包括：

测序数据存储装置：用于存储低深度全基因组测序数据；

数据分析装置：用于对低深度全基因组测序数据进行ROH分析，所述ROH分析的流程包括：

基于现有SNP信息，将低深度全基因组测序数据进行比对，获取目标SNP位点的等位基因allele，所述allele信息包括：每个SNP位点总的测序深度、每个SNP位点A allele和B allele的测序深度；所述A allele为目标SNP一种基因型，其在人群中的参考频率为p，所述 B allele为该目标SNP的其他基因型，其在人群中的参考频率为1-p；

将多个在人类基因组中位置相邻的SNP序列合并成一个窗口bin，选取测序深度相同且不低于2的SNP序列结果，统计每个bin在该测序深度下的SNP数n以及其中包含杂合信息 SNP数k，求得样本i在第j个bin的SNP混合杂合度指标SMHS_ij；

计算不同测序深度下的SNP混合杂合度指标SMHS_ij，将连续t个SMHS_ij低于该测序深度下期望值的bins作为候选ROH区段；

基于候选ROH区段的基因拷贝数信息，排除拷贝数＝1的候选ROH区段，得到最终的ROH结果；

结果输出装置：用于输出数据分析结果。

在一些实例中，所述SNP为只存在2种等位基因的双等位基因SNP位点。

在一些实例中，所述SNP的最小等位基因频率MAF值为0.01～0.5。

在一些实例中，选择测序深度为2的SNP结果并基于测序深度分别计算样本的SMHS_ij值

式(1)中，n为当前测试样本i在第j个bin中SNP被测序深度为2的个数，k为同时检测到A allele和B allele的个数；

筛选存在连续t个SMHS_ij<0.01的bins作为候选ROH区段，t个SMHS_ij<0.01的bins的长度之和要大于样本测序深度λ对应的最小ROH检测长度L。

在一些实例中，结合最低SNP位点数和SNP位点在基因组中分布的密度，推算出当前测序深度λ下可检出的ROHs的最小长度L。

在一些实例中，最小能检出的ROHs长度L(Mb)为：

式中，λ为样本的平均测序深度。

在一些实例中，根据泊松分布可以计算不同测序深的SNP的概率进而计算测序深度为定值的SNP个数n。

在一些实例中，测序深度为2时，则窗口bin对应的SNP个数N需满足：

式中，λ为样本的平均测序深度。

一个bin中的SNP数，可以根据ROH检出分辨率的要求和测序深度调整。在一些实例中，合并的一个bin中，含有200～2000个SNP位点，或长度为0.1～1Mb。

出于成本控制考虑，在一些实例中，所述低深度全基因组测序数据的平均测序深度低于 30X。

出于检测的准确性考虑，在一些实例中，所述低深度全基因组测序数据的平均测序深度优选不低于0.12X。

本发明的第二个方面，提供：

一种基于低深度测序数据的多倍体异常的检测***，包括：

测序数据存储装置：用于存储低深度全基因组测序数据；

数据分析装置：用于对低深度全基因组测序数据进行多倍体分析，所述多倍体分析的流

程包括：

ROH分析：按本发明第一个方面所述的ROH分析流程对样本i进行ROH分析；

基因拷贝数异常分析：识别样本i中的拷贝数异常区域；

计算样本i的SNP混合杂合度指标SMHS_i：

式(2)中，n为样本i中去除了ROH区域和拷贝数异常区域后剩余的被测序深度为2的SNP的个数，k为其中同时检测到A allele和B allele的个数；

利用验证的多例已知二倍体样本作为参考数据集推算总体样本的SMHS_i的均值μ和标准差S，求得样本i的杂合率Z值：

SMHS_zscore_i大于3判断为多倍体异常；

结果输出装置：用于输出数据分析结果。

在一些实例中，使用人类短串联重复序列STR方法验证已知二倍体样本。

在一些实例中，所述已知二倍体样本的例数不低于30例，优选为60～200例、100例。

出于成本控制考虑，在一些实例中，所述低深度全基因组测序数据的平均测序深度低于 30X。

出于检测的准确性考虑，在一些实例中，所述低深度全基因组测序数据的平均测序深度优选不低于0.12X。

本发明的有益效果是：

本发明的一些实例，打破了现有理论的局限，在不改变现有的CNV-seq实验方案、测序量的前提下，对其功能进行了拓展，利用极低测序深度的SNP信息进行多倍体和基因组纯合区域ROH检测，开发出一整套基于低测序深度全基因组测序CNV-seq数据的胎儿多倍体、基因组纯合区域ROH的检测***。

本发明的一些实例，大幅减少了ROH和多倍体检测的成本。

附图说明

图1是低测序深度下的SNP位点的测序深度分布情况；

图2是不同类型样本SNP位点在测序深度等于2时检测到杂合信号的期望概率；

图3是参考数量样本量对标准误和样本方差的比值的影响；

图4是样本A18T1402 CMA方法对应的基因组拷贝数和SNP信息；

图5是样本A18T1402本发明一些实例的方法对应的基因组拷贝数和SNP信息；

图6是样本EM000457T CMA方法对应的基因组拷贝数和SNP信息；

图7是样本EM000457T本发明一些实例的方法对应的基因组拷贝数和SNP信息；

图8是样本EM000378B CMA方法对应的基因组拷贝数和SNP信息；

图9是样本EM000378B本发明一些实例的方法对应的基因组拷贝数和SNP信息；

图10是样本A17T4760 CMA方法对应的基因组拷贝数和SNP信息；

图11是样本A17T4760本发明一些实例的方法对应的基因组拷贝数和SNP信息。

具体实施方式

方便理解起见，统一将MAF值对应的allele称为“A”allele，“A”allele的人群频率为p；另外一种allele称为“B”allele，“B”allele的人群频率为q，q＝1-p。如果比对结果显示某个SNP位点同时存在“A”allele和“B”allele，则判为杂合位点，反之成为纯合位点。需要注意的是，本文在极低测序深度下判断的杂合位点指的是在当前测序深度下对该位点是否为杂合的一种判断，而非真实情况下的基因型，因为极低测序深度下真实状态的AB 杂合型有很大概率被检测成纯合状态，本文统一用(1,1)表示两个等位基因均被检测到的情况。

本发明的理论基础：

在极低测序深度下的全基因组高通量测序数据中(如CNV-seq数据)，绝大部分测序的 SNP位点的深度都是1X，少数SNP位点的测序深度为2X或3X，如图1所示。

根据泊松分布可知在平均测序深度为X的时候，对应测序深度为某个定值的SNP占比。随着测序深度的增加，测序深度≥2的有效SNP位点显著增加，检测准确性逐步提升。由于不同测序深度的SNP位点，检测到杂合信号的概率不同，因此建议在不同测序深度下的SNP 的位点分别统计检验来判断样本是否存在ROH或多倍体异常。不同平均测序深度下SNP测序深度的预期占比如下表所示：

本技术使用的测序深度平均约0.2X左右，测序深度≥3的SNP占比仅约0.11％，因此优选depth＝2的测序深度的SNP进行分析。

不同类型的样本或基因组区域(如二倍体、三倍体、ROH纯合区域)对应的基因型的种类和占比是不同的，SNP-array在已知每个SNP位点的基因型的前提下，可以简单的计算当前样本/基因组区域出现的基因型的种类即可判断样本的类型。但是在在极低测序深度下(平均测序深度≈0.2X)，虽然每个SNP位点的基因型不可知，但是通过理论推导可知，同样在测序深度等于2的时候，不同类型的样本/基因组区域观察到的杂合信号P(1，1)的概率不同。利用二项分布计算不同样本类型(ROH、二倍体、三倍体、四倍体)的所有基因型以及对应的概率，并进一步计算不同基因型在2X覆盖深度下能同时检测到A和B位点的期望概率(即低测序深度下检测到的杂合率)，从原理上充分证实了不同类型的样本在低覆盖深度下的杂合率期望值是存在显著差异的。

不同类型样本SNP位点在测序深度等于2时检测到杂合信号的概率(即SMHS)如下表或图2所示：

三倍体和四倍体胎儿在发育过程中，由于减数***一期的时候有少量同源染色体会发生交叉互换(crossing-over)现象，导致三倍体样本有少量基因组区域并非来源于不同的祖源，所以实际上三倍体的杂合率会比期望值稍低，但显著高于二倍体样本(图2)。

在理论基础上，本发明开创性地提出SNP混合杂合度指标(snp_mixed_heterozygosit_score，简称SMHS)，表明在极低测序深度下(平均测序深度≈0.2X)，虽然每个SNP位点的基因型不可知，但是通过合并计算N个SNP杂合信息，计算SNP混合杂合度指标SMHS_ij(SMHS为snp mixed heterozygosity ratio缩写)，可以用于确定ROH。

选取测序深度相同且不低于2的SNP序列结果，将N个在基因组位置上相邻的SNP序列合并成一个窗口bin，统计每个bin在该测序深度下的SNP数n以及其中包含杂合信息SNP数k，求得样本i在第j个bin的SNP混合杂合度指标SMHS_ij；

当前测试样本i在bins j，有n个SNP被测序深度为2，其中k个同时检测到A和B等位基因，则样本i在bin j的SNP混合杂合度指标SMHS_ij：

计算不同测序深度下的SNP混合杂合度指标SMHS_ij，将连续t个SMHS_ij低于该测序深度下期望值的bins作为候选ROH区段。连续t个SMHS_ij<0.01的bins的长度之和要大于样本测序深度λ对应的最小ROH检测长度L；优选的，结合最低SNP位点数和SNP位点在基因组中分布的密度，推算出当前测序深度λ下可检出的ROHs的最小长度L。

需要注意的是，测序深度不同的时候，SMHS_ij的期望值是不同的，因此需选定测序深度，然后相同测序深度下的SMHS_ij值进行统计检验来判断样本是否存在ROH或多倍体异常。本发明使用的数据，测序深度平均0.2X左右，测序深度≥3的SNP占比仅约0.11％，因此优选 depth＝2的测序深度的SNP进行分析。

多倍体检测计算

进行ROH分析和基因拷贝数异常分析之后，计算样本i的SNP混合杂合度指标SMHS_i：

式(2)中，n为样本i中去除了ROH区域和拷贝数异常区域后剩余的被测序深度为2的SNP的个数，k为其中同时检测到A和B allele的个数；

根据参考数据集样本杂合率SMHS_i的均值μ和标准差S，求得样本i的杂合率Z值：

SMHS_zscore_i大于3判断为多倍体异常。

利用人类短串联重复序列(包含特异性的遗传STR标记)简称STR方法验证的100例已知二倍体样本作为参考数据集推算总体样本的SMHS_i的均值和方差。利用shapiro.test检验样本的杂合率值SMHS_i是否符合正态分布，检验结果为“W＝0.98325,p-value＝0.1285”，表明SMHS_i符合正态分布，适用Z检验方法。因此可以根据参考数据集样本杂合率SMHS_i的均值μ和标准差S，求得样本i的杂合率Z值。

已知二倍体样本的例数n对据参考数据集样本杂合率SMHS_i的均值μ和标准差S推算的影响：

由于抽样误差，参考数据集计算出来的样本均数和样本均数的标准差只能是逼近但是无法等同于总体均数和总体标准差。样本均数的标准差也称为样本均数的标准误

反映了样本均数和总体均数的差异。样本均数的标准误为：

即样本均数的标准误和样本量n的平方呈反比，如图3所示，样本量在1～30的时候，随着样本量的增加，标准误和样本方差的比值急速下降，当样本量>30时，比值降低的速度明显趋于平缓，即增加样本量对于降低标准误的效能在逐步降低，因此通常情况下构建参考数据集要求样本量>＝30例即可。当n足够大，如n>60时，可以用样本统计量直接作为总体参数的估计值。此时样本方差S逼近总体方差σ；而样本均数

逼近总体均数μ。在本发明的一些实施例中，采用样本量n＝100构建参考数据集，样本均数标准误为样本方差的10％，远小于总体方差。

最终，本发明一些实例在不改变现有的CNV-seq实验方案、测序量的前提下，增加了多倍体和基因组纯合区域ROH检测功能，并进一步开发出基于低测序深度全基因组测序(CNV-seq)数据的胎儿多倍体、基因组纯合区域ROH的一整套解决方案、检测***和建模方法。

详细的理论依据

1)推算染色体ROH区域检测到杂合信号的概率：

ROH区域AB位点的概率：P(genotype＝AB)＝0，因此ROH当测序深度＝2的时候，二倍体样本同时检测到A和B allele的概率P(AB)的期望值为0。

2)推算二倍体样本检测到杂合信号的概率即SMHS

二倍体样本中，biallele SNP位点有“AA”、“AB”和“BB”三种组合，根据Hardy-weinberge equilibrium三种基因型的概率分别为P(AA)＝p²,P(AB)＝2pq,and P(BB)＝q²。

当测序深度为2的时候，三种基因型检测到杂合信号的概率为：

P((1,1)|AA)＝0

P((1,1)|BB)＝0

根据全概率公式，二倍体样本检测到杂合信号的概率即P(1,1|diploid)为：

3)推算三倍体样本检测到杂合信号的概率：

三倍体样本中，biallele SNP位点有“AAA”、“AAB”、“ABB”和“BBB”四种组合，概率分别为P(AAA)＝p³,P(AAB)＝3p²q,P(ABB)＝3pq²,andP(BBB)＝q³.

当测序深度为2的时候，三种基因型检测到杂合信号的概率为：

P((1,1)|AAA)＝0

P((1,1)|BBB)＝0

根据全概率公式，三倍体样本(三条染色体来源于不同祖源染色体的区间)检测到杂合信号的概率P(1,1|triploid)为：

4)推算四倍体样本检测到杂合信号的概率：

四倍体样本中，biallele SNP位点有“AAAA”、“AAAB”、“AABB”、“ABBB”和“BBBB”五种组合

四倍体样本中只有基因型AAAB、AABB和ABBB的点有可能检测到杂合信号

P((1，1)|AAAA)＝0

P((1，1)|BBB)＝0

根据全概率公式，四倍体样本(四条染色体来源于不同祖源染色体的区间)检测到杂合信号的概率P(1，1|tetraploid)为：

下表汇总不同类型样本SNP位点在测序深度等于2时检测到杂合信号的概率：

下面结合实例，进一步说明本发明的技术方案。

样本收集

收集150例产前诊断标本，样本类型包括羊水、绒毛、脐血、外周血、流产组织等；对照方法的结果包括：

STR方法验证的二倍体样本100例作为参考数据集(reference)，用于计算正常样本SNP 杂合率的均值(mean)和标准差(SD)。STR方法分析人类短串联重复序列(ShortTandem Repeat, STR)，这些STR序列的核心一般由2～6个碱基组成，该核心序列在不同个体间的串联重复次数不同而呈现出长度多态性，是一种高度多态性的遗传标记，适用于个体识别、亲权鉴定和遗传分析等。本实施例所用的STR方法为阅微基因的“Microreader^TM21Direct ID System”***，共检测21个STR位点，采用5色荧光标记、多重扩增检测STR位点，最终使用配套的遗传分析仪器检测STR信号并进行分型。STR方法可鉴别样本是否属于多倍体异常，但是由于所用STR标记过少无法提供ROH的结果，因此STR验证的样本非多倍体样本仅用于构建参考数据集。

染色体微阵列分析CMA方法验证样本40例。由于CMA方法不仅可以判断样本是否多倍体，同时明确提供样本拷贝数变异或ROH区域的详细位置，因此该方法验证的样本可用于准确评估方法的准确性。CMA方法验证的样本包括：三倍体样本14例，二倍体含ROH的样本24例，二倍体不含ROH的样本2例。本实施例中所用的CMA方法为affymetrix公司的“Cytoscan^R750K”CMA，包括550,000个非多态性探针以及200,000个可用于基因分型的 SNP探针，共计750,000个探针。

超低深度(≈0.2X)全基因组测序

提取样本的基因组DNA并进行DNA片段化，片段化的方法主要有酶切法和超声打断法，在这里优选酶切方法进行DNA片段化。然后对片段化的DNA进行文库构建，主要包括末端修复、加接头、PCR富集、文库定量、文库质检等过程。文库质检通过后进行全基因组测序。

比对和过滤

测序结束后，每个样本会产生大量读长(reads)。

比对和过滤也可以采用已有的方法或软件进行，去除未比对上参考基因组的读长(unmapped reads)；去除mapping质量较差(比对质量得分MAQ≤10)的比对结果；删除同时比对到基因组中多个位置的reads。

利用BWA比对软件将读长比对到人类参考基因组序列(参考基因组序列版本hg19，比对参数bwa mem-t 20-B 1-O 1-L 1–T-M reference)。为了提高比对结果的准确性，降低噪音，需要对比对结果进行过滤，包括以下步骤：去除未比对上参考基因组的读长(unmapped reads)；去除mapping质量较差(比对质量得分MAQ≤10)的比对结果；删除同时比对到基因组中多个位置的reads；比对到基因组同一位置的多条reads(duplicatereads) 仅保留比对质量得分最高的一条记录，从而避免来自同一个游离DNA片段的reads被重复计算多次。

实施例中的140例样本中，平均每个样本测序获得的reads数(raw reads number)为5.90M，过滤低质量的比对结果和重复序列后唯一比对读长数(unique reads number)为4.45M，基因组覆盖度为19.24％即0.19X。

识别染色体拷贝数异常区域

识别染色体拷贝数异常区域可以使用已有的方法进行，或按如下操作进行：

根据过滤后的unique reads的比对结果，统计基因组每20kb滑动窗口的readsnumber数并进行GC校正。染色体以20KB为步长，计算GC％和reads number。计算每0.1％GCbin 的reads number中位数M_i以及所有GC bin的均值M_global，校正系数W＝M_i-M_global，之后将原来的reads数减去校正系数得到校正后的reads number。

进一步将20KB窗口合并成合适大小的窗口bins(例如100kb、1Mb bins、每条染色体)，并统计落在bins的reads number(RN_i)，然后用常染色体所有窗口的reads number对RN_i进行归一化，求得样本i在第j个bin的reads ratio值(RR_ij)。

计算每条染色体的RR_ij值，并进而计算Zscore，将Zscore>3的染色体判断为三体；Zscore<-3的染色体判断为单体。

使用环状二元分割CBS(Circular Binary Segmentation)算法(R包DNAcopy)识别除整倍体以外的更小的拷贝数变异[8]。

SNP位点信息的获取

dbSNP数据库汇总所有SNP位点以及SNP位点对应的ref allele和alt allele在人群中的参考频率(或基于其他已有数据确定参考频率)。方便比较起见，在本实例中，选择只存在两种等位基因的SNP位点(biallele SNP)以及MAF(minor allele frequency)在某个区间的 SNP位点。在本实例中，优选MAF值在0.01～0.5区间的SNP位点，共计约1200万个SNP。

利用samtools mpileup软件分析比对结果，获取目标SNP位点的allele信息，具体可以得知每个SNP位点总的测序深度以及“A”allele和“B”allele的测序深度。当然，也可以使用其他方法确定目标SNP位点的allele信息。

实施例中的140例样本中，平均每个样本基因组覆盖度(coverage)为19.24％、唯一比对读长均数(unique reads number)为4.45M、其中比对到SNP的reads均数为1.15M、有读长覆盖的SNP均数102.84万、测序深度等于2的SNP均数10.52万。如下表所示：

计算窗口内的SNP混合杂合度指标SMHS_ij

因为杂合位点在depth＝1的时候无法观察到杂合现象，所以过滤掉depth＝1的结果，选取 depth>＝2的SNP结果。本文中实施例的平均测序深度为0.2X，depth>＝3的SNP占比只有约 1‰，且depth＝3的时候各类型样本杂合率的期望值和depth＝2不同，因此优选depth＝2的SNP 进行分析。但是随着测序深度的增加，建议增加depth>＝3深度的位点，但是需分别计算每个 depth深度下的杂合率snp_mixed_heterozygosity_ratio，之后进行归纳统计。

进一步将N个SNP合并成合适大小的bin(例如2000个、5000个、整条染色体)，然后统计每个窗口测序深度在某个测序深度下，优选depth＝2的SNP数snpnum_i以及其中包含杂合信息(单个SNP位点同时发现A、B两种allele)的SNP数heter_snpnum_i，求得样本i 在第j个bin的SNP混合杂合度指标SMHS_ij。

假设当前测试样本i在bins j，有n个SNP被测序深度为2，其中k个同时检测到A和B等位基因，则样本i在bin j的SNP混合杂合度指标SMHS_ij的公式为：

实施例中的140例样本中，平均每个样本基因组覆盖度(coverage)为19.24％，唯一比对读长均数(unique reads number)为4.45M，其中比对到SNP的reads均数为1.15M，有读长覆盖的SNP均数102.84万，测序深度等于2的SNP均数10.52万，测序深度等于2的杂合SNP 均数1.41万。如下表所示：

识别基因组纯合区域ROH

基因组纯合区域(absence of heterozygosity；ROH)描述的是基因组一段区域内SNP只有纯合子而没有杂合子的状态，也通常称为杂合性缺失(loss of heterozygosity；LOH)。理论上，基因组纯合区域ROH区域所有的SNP都是纯合子，观察不到杂合的SNP位点，即杂合率heter_snpratio_i的期望值为零。而二倍体AB基因型概率P_AB＝2*p*q，二倍体2XSNP位点观察到杂合率的期望值为p*q。因此ROH和二倍体样本的杂合率的期望值时显著差异的。

根据全概率公式，单亲二倍体、二倍体、三倍体样本的SNP杂合率的期望值和SNP在人群中的杂合频率有关。

染色体ROH区域检测到杂合信号的概率：

ROH区域AB位点的概率：P(genotype＝AB)＝0，因此ROH当测序深度＝2的时候，二倍体样本同时检测到A和B allele的概率P(AB)的期望值为0。

P(1，1|ROH)＝0

二倍体样本检测到杂合信号的概率：

由于低测序深度下单个SNP位点的信息不足以推断杂合度，那么需要合并N个SNP的测序结果来计算SNP位点的混合杂合度SMHS。考虑到基因组中SNP位点的分布不是绝对均匀的，为了保证每个分析窗口(bins)都得到足够数量的SNP位点，建议根据SNP位点数代替基因组长度来划分窗口。

在实际应用中，需要结合测序深度选择窗口大小以确保窗口内的有效SNP位点数(深度大于等于2)大于最小值，因为SMHS_ij值的标准差和窗口内目标测序深度的有效SNP位点数的开方呈反比。正常二倍体区域的SMHS_ij值的期望值约0.13，那么在3倍SD的情况下，单个窗口的标准差需要控制在SD<0.13/3即0.04333以内。根据Linderberg中心极限定理，窗口内的SMHS_ij值服从N(μ，μ(1-μ)/n)正态分布，因此可推出每个窗口内落在目标测序深度的SNP数应大于60个

本技术共筛选到587万个MAF值在0.01～0.5之间的biallele SNP位点，其中位于常染色体的SNP位点572万个。常染色体总长2684.57Mb，因此优选每2000个SNP位点为一个分析窗口bin，相当于1Mb的基因组长度。

下面探索测序深度对应的适合的窗口bin包含的SNP个数的关系。

假设，窗口bin中包含N个SNP位点，测序深度为λ，根据泊松分布可以计算不同测序深的SNP的概率进而计算测序深度为定值的SNP个数n。

下面的公式计算测序深度等于2的时候，对应的SNP个数n和概率P0:

由于已知要求n>＝60，则窗口bin对应的SNP个数N需满足：

本技术共筛选到587万个MAF值在0.01～0.5之间的biallele SNP位点，其中位于常染色体的SNP位点572万个，相当于平均每Mb基因组中包含2130个SNP位点。已知样本的测序深度为λ，则最小能检出的ROHs长度L(Mb)为：

下表计算除了在已知测序深度λ对应的深度等于2的SNP概率，以及满足n>＝60时所需的最低SNP位点数，最终结合最低SNP位点数和SNP位点在基因组中分布的密度可推算出当前测序深度λ下可检出的ROHs的最小长度L。

根据《产前遗传学诊断拷贝数变异(CNV)和纯合状态(ROH)数据解读及报告规范的专家共识》指南：涉及chr6、chr7、chr11、chr14、chr15、chr20其中某一条染色体存在>＝5Mb(位于染色体末端)或>＝10Mb(非染色体末端)而其他染色体上未见ROH时，建议报告并提示进行UPD检测，原则上对临床意义不明确的ROH不建议报告。因此，如果需要检测>＝5Mb 的ROHs，则最小的测序深度应该为0.12X，ROH区间至少需包含>＝9396个SNP位点；当测序深度增加到0.28X时，可检测>＝1Mb的ROHs，ROH区间包含>＝2028个SNP位点。

测序深度对可检测ROH长度的影响

假设当前测试样本i在bins j，有n个SNP被测序深度为2，其中k个同时检测到A和B等位基因，则样本i在bin j的SNP混合杂合度指标SMHS_ij：

因此，识别基因组纯合区域ROH任务可以简化为识别N个连续的SMHS_ij≈0的bins。由于真实的测序数据会存在一定的测序错误，在这里设定候选ROH位点的阈值为SMHS_ij<＝0.05。首先筛选SMHS_ij<＝0.05区域，然后向前向后进行检索，如果下一个符合条件的位点和当前位点的距离<＝2则将相应的区域合并，之后再进行向前向后检索，在所获得的位置信息稳定后停止检索。最后筛选存在连续t个SMHS_ij<0.01的染色体区段(bins)作为ROH区段。t的最小值根据样本的测序深度来确定，根据上表所示t乘窗口bin的长度的总长需大于样本测序深度λ对应的最小ROH检测长度L。二倍体样本在拷贝数＝1即杂合缺失的位置SNP都显示为 ROH，因此核对该位置的拷贝数信息，如果是杂合缺失(即微缺失)位置的ROH直接出具缺失del的结果；只对拷贝数正常即拷贝数等于2的区域出具ROH的报告。

根据《产前遗传学诊断拷贝数变异(CNV)和纯合状态(ROH)数据解读及报告规范的专家共识》指南：涉及chr6、chr7、chr11、chr14、chr15、chr20其中某一条染色体存在>＝5Mb(位于染色体末端)或>＝10Mb(非染色体末端)而其他染色体上未见ROH时，建议报告并提示进行UPD检测，原则上对临床意义不明确的ROH不建议报告。

由于多倍体样本中理论上无ROH的区域，临床检测多倍体的ROH区域无意义，而STR验证的二倍体并不确定ROH所在的位置信息。因此使用CMA方法验证的26例二倍体样本评估两种方法对ROH检测的一致性。

表1展示了比较了基于SNP-混合杂合率和CMA方法的ROH检测结果。表中，

整体来看，两种方法对ROH区域检测的一致性高达100％，但是由于SNP-array和本实例中SNP覆盖的具***置有所差异，分析出来的具***置会有些差别。例如样本EM000457T，使用CMA方法的结果为7p12.2p11.1(49087710-58019983)hmz, 7q11.21q11.22(62569501-70689665)hmz,(9)x3，而使用本实例的SNP-混合杂合度的方法的结果为47,XY,+9,7p12.2-q11.23(49751619-72884378)hmz。两种方法都判断 7p12.2-q11.23(49751619-72884378)这个区段存在ROH，但是CMA方法在该区间有大片段区域无SNP探针覆盖，所以将ROH区域分成了两段。两种方法检测的重叠率在70％～100％之间，不一致的原因主要是部分ROH区域位于CMA方法SNP探针的检测盲区(未设计相应的SNP探针)而被分成两段。重叠性的指标overlap计算公式如下

样本A18T1402，CMA方法提供的结果为：arr(1-22,X)×2hmz，对应的基因组拷贝数和 SNP信息如图4所示；

样本A18T1402，基于本实例SNP杂合率分析的CNV-Seq检测结果为单亲二倍体，检测结果和CMA方法完全一致。对应的基因组拷贝数和SNP信息如图5所示；

样本EM000457T，CMA方法提供的结果为：arr(1-8,10-22)×2,(X,N)×1,(9)×3，7p12.2p11.1(49087710_58019983)hmz,7q11.21q11.22(62569501_70689665)hmz，对应的基因组拷贝数和SNP信息如图6所示；

样本EM000457T，基于本专利SNP杂合率分析的CNV-Seq检测结果为 47,XY,+9,7p12.2-q11.23(49751619-72884378)hmz，检测结果和CMA方法完全一致。对应的基因组拷贝数和SNP信息如图7所示；

样本EM000378B CMA方法提供的结果为：

arr(1-22)×2,(X,N)×1,1p32.2p31.3(56716570_66496108)hmz

1p35.3p34.3(28400877_38829065)hmz

2q12.2q21.1(106045084_132094747)hmz

4q24q27(106352104_122499793)hm,4q31.3q32.3(153113725_166664618)hmz

5p14.2p13.1(23709942_40437845)hmz

5q14.1q21.3(77259148_106713065)hmz

6q21q22.31(108975645_123967785)hmz

9p24.3p21.3(1095921_23587811)hmz

12p12.2p11.1(20607265_34761150)hmz

12q11q13.11(38190102_47460200)hmz

16p12.2p11.2(22133216_31935367)hmz

16q11.2q22.3(46504466_73641621)hmz

17p12p11.2(13768916_22170994)hmz

17q11.1q12(25309336_33005958)hmz。对应的基因组拷贝数和SNP信息如图8所示；

样本EM000378B，基于本实例SNP杂合率分析的CNV-Seq检测结果为 46,XY,1p32.2-p31.3(57065523-66716272)hmz

1p35.3-p34.2(28834322-40407973)hmz

2q12.1-q21.2(104469417-133378380)hmz

4q31.3-q32.3(153569494-165727163)hmz

4q24-q27(104653213-121400901)hmz

5q14.1-q21.3(78680234-108109454)hmz

5p14.1-p12(25146518-42747892)hmz

6q21-q22.31(110591731-123094973)hmz

7p22.2-p21.3(3705246-9804912)hmz

9pter-p21.3(1657050-25130611)hmz

12p12.3-q12(19917606-46029214)hmz

16p12.2-q22.3(23596334-73025281)hmz

17p12-q12(13974486-33838288)hmz。

基于SNP混合杂合度方法的检测结果和CMA方法完全一致。此样本检出多条染色体均存在大片段ROH怀疑父母存在亲缘关系所致。对应的基因组拷贝数和SNP信息如图9所示。

识别多倍体样本

根据全概率公式，单亲二倍体、二倍体、三倍体样本的SNP杂合率的期望值和SNP在人群中的杂合频率有关。筛选dbSNP数据库中最小等位基因频率(MAF,minor allelefrequency)>0.01的1200万个SNP位点。本文中统一将最小等位基因频率MAF值对应的allele 称为“A”allele，“A”allele的人群频率为p；另外一种allele称为“B”allele，“B”allele 的人群频率为q，q＝1-p。

染色体ROH区域检测到杂合信号的概率：

ROH区域AB位点的概率：P(genotype＝AB)＝0，因此ROH当测序深度＝2的时候，二倍体样本同时检测到A和B allele的概率P(AB)的期望值为0。

P(1，1|ROH)＝0

二倍体样本检测到杂合信号的概率：

三倍体样本检测到杂合信号的概率：

四倍体样本检测到杂合信号的概率：

假设当前测试样本i，去除了通过步骤5和7识别的拷贝数异常区域和ROH区域，剩余 n个SNP被测序深度为2，其中k个同时检测到A和B等位基因，则样本i的SNP混合杂合度指标SMHS_i：

进一步利用100个已知的二倍体样本计算样本杂合率SMHS_i的均值和标准差则可求得每个样本的杂合率Z值。

样本的SNP混合杂合度指标SMHS_iZ值：

最终将杂合率Z值>3的样本判断为多倍体异常。

100例正常样本的杂合率均值和标准差分别为0.1312和0.003049。

14例多倍体样本，基于低测序深度的SNP杂合率方法全部检出，对于多倍体的判断结果和CMA方法达到100％一致。由于CNV-seq方法对于拷贝数检测的灵敏性和准确性更高，因此样本EO900670T和EO901093T都比CMA方法额外检出了1个<1Mb的拷贝数异常。具体如表2所示。

表2展示了SNP-混合杂合模型和CMA方法两种方法检测多倍体的一致性。提供了每个样本的样本ID、样本类型、两种方法的多倍体检测结果。表2详细地提供了每个样本基于SMHM混合杂合度模型检测的杂合率值SMHS、SMHS的Z值、多倍体判断结果以及最终检测结果。

样本A17T4760CMA方法提供的结果为：arr(1-22)×3,(X)×2,(Y)×1，对应的基因组拷贝数和SNP信息如图10所示；

样本A17T4760，基于本实例SNP杂合率分析的CNV-Seq检测结果为69,XXY，检测结果和CMA方法完全一致。对应的基因组拷贝数和SNP信息如图11所示。

表3展示了实施例中150例样本倍型分析结果。其中包括STR方法验证的二倍体样本100 例、CMA方法验证的单亲二倍体样本2例、CMA方法验证的二倍体样本24例和CMA方法验证的三倍体样本16例。150例实施例样本，基于混合杂合度模型SMHM检测的倍型分析结果均和已知倍型一致。

表4展示了已知ROH结果的26例二倍体样本的详细分析结果。包括每个ROH区域所在的染色体、起始位置、终止位置、SMHS杂合率、对应的拷贝数信息以及模型判断结果。

表3样本倍型分析结果

注：CMA指基于SNP-array的染色体微阵列分析方法CMA

表4样本ROH分析结果

Claims (14)

1.一种基于低深度测序数据的基因组纯合区域ROH的检测***，包括：

测序数据存储装置：用于存储低深度全基因组测序数据；

数据分析装置：用于对低深度全基因组测序数据进行ROH分析，所述ROH分析的流程包括：

基于现有SNP信息，将低深度全基因组测序数据进行比对，获取目标SNP位点的等位基因allele信息，所述allele信息包括：每个SNP位点总的测序深度、每个SNP位点A allele和B allele的测序深度；所述A allele为目标SNP一种基因型，其在人群中的参考频率为p，所述B allele为该目标SNP的其他基因型，其在人群中的参考频率为1-p；

将多个在人类基因组中位置相邻的SNP序列合并成一个窗口bin，选取测序深度相同且不低于2的SNP序列结果，统计每个bin在该测序深度下的SNP数n以及其中包含杂合信息SNP数k，求得样本i在第j个bin的SNP混合杂合度指标SMHS_ij；

选择测序深度为2的SNP结果并基于测序深度分别计算样本的SMHS_ij值

式(1)中，n为当前测试样本i在第j个bin中SNP被测序深度为2的个数，k为同时检测到Aallele和B allele的个数；

筛选存在连续t个SMHS_ij<0.01的bins作为候选ROH区段，连续t个SMHS_ij<0.01的bins的长度之和要大于样本测序深度λ对应的最小ROH检测长度L；

基于候选ROH区段的基因拷贝数信息，排除拷贝数＝1的候选ROH区段，得到最终的ROH结果；

结果输出装置：用于输出数据分析结果。

2.根据权利要求1所述的检测***，其特征在于：所述SNP为只存在2种等位基因的双等位基因SNP位点。

3.根据权利要求1所述的检测***，其特征在于：所述SNP的最小等位基因频率MAF值为0.01～0.5。

4.根据权利要求1所述的检测***，其特征在于：结合最低SNP位点数和SNP位点在基因组中分布的密度，推算出当前测序深度λ下可检出的ROHs的最小长度L。

5.根据权利要求1所述的检测***，其特征在于：合并的一个bin中，含有200～2000个SNP位点，或长度为0.1～1Mb。

6.根据权利要求1所述的检测***，其特征在于：所述低深度全基因组测序数据的平均测序深度低于30X。

7.根据权利要求6所述的检测***，其特征在于，所述低深度全基因组测序数据的平均测序深度不低于0.12X。

8.一种基于低深度测序数据的多倍体异常的检测***，包括：

测序数据存储装置：用于存储低深度全基因组测序数据；

数据分析装置：用于对低深度全基因组测序数据进行多倍体分析，所述多倍体分析的流程包括：

ROH分析：按权利要求1～7任一项所述的ROH分析流程对样本i进行ROH分析；

基因拷贝数异常分析：识别样本i中的拷贝数异常区域；

计算样本i的SNP混合杂合度指标SMHS_i：

式(2)中，n为样本i中去除了ROH区域和拷贝数异常区域后剩余的被测序深度为2的SNP的个数，k为其中同时检测到A allele和B allele的个数；

利用验证的多例已知二倍体样本作为参考数据集推算总体样本的SMHS_i的均值μ和标准差S，求得样本i的SMHS_zscpre_i：

SMHS_zscore_i大于3判断为多倍体异常；

结果输出装置：用于输出数据分析结果。

9.根据权利要求8所述的检测***，其特征在于：使用人类短串联重复序列STR方法验证已知二倍体样本。

10.根据权利要求8或9所述的检测***，其特征在于：所述已知二倍体样本的例数不低于30例。

11.根据权利要求10所述的检测***，其特征在于，所述已知二倍体样本的例数为60～200例。

12.根据权利要求11所述的检测***，其特征在于，所述已知二倍体样本的例数为100例。

13.根据权利要求8所述的检测***，其特征在于：所述低深度全基因组测序数据的平均测序深度低于30X。

14.根据权利要求13所述的检测***，其特征在于，所述低深度全基因组测序数据的平均测序深度不低于0.12X。

Images