CN112562787B - 一种基于ngs平台的基因大片段重排检测方法 - Google Patents

Abstract

本发明涉及一种基于NGS平台的基因大片段重排检测方法，具体通过比较样本与阴性样本在目标区域的测序深度的差异来确定基因大片段重排的区域。提高了检测分辨率和灵敏度，有效的解决了NGS中常规的CNV检测方法分辨率低无法检测LGR的问题。

Description

一种基于NGS平台的基因大片段重排检测方法

技术领域

本申请涉及生信分析领域，具体涉及一种基于NGS平台的基因大片段重排检测方法。

技术背景

基因的大片段重排(Large Genomic Rearrangements，LGR)指一个或多个外显子的重复或缺失。重排类型大部分为缺失，也存在二倍、三倍重复等，这些变异往往引起读码框偏移，导致蛋白结构与功能的异常。

常见的致病性变异多为单碱基变异(Single Nucleotide Variation，SNV)或短片段插入缺失(Small Insertion and Deletion，INDEL)引起的移码突变和无义突变，这些变异类型可以在下一代测序(Next Generation Sequencing，NGS)中稳定检测。此外，LGR也是一种重要的致病性变异类型，很多研究表明其在肿瘤以及遗传疾病领域具有重要的临床价值。例如，BRCA基因致病性变异使女性发生乳腺癌的风险提高5倍，发生卵巢癌的风险提高10～30倍。且BRCA基因的LGR致病性变异的发生频率较高，在乳腺癌高风险人群中的检出频率为2.4％，占所有BRCA致病性变异的9.9％，在普通人群中的检出频率为0.5％，占所有BRCA致病性变异的5.9％。在临床上，携带BRCA基因致病性变异的乳腺癌和卵巢癌患者在接受铂类治疗中获益，具有更好的总生存率(Overall Survival，OS)、无进展生存率(Progressive-Free Survival，PFS)和总缓解率(Overall Response Rate，ORR)。以上研究显示，当常规的NGS检测结果(SNV，INDEL)为阴性时，还应进行LGR的检测，防止出现致病性的LGR漏检问题，以便制定相应的遗传管理措施和治疗方案。然而，目前常规的NGS检测方法由于分辨率较低尚不足以检测LGR，因此开发一种基于NGS技术的LGR检测方法具有重要意义。

多重链接依赖探针扩增技术(Multiplex Ligation-dependent ProbeAmplification，MLPA)是最常用的检测LGR的方法，也是目前LGR检测的金标准。其基本原理是DNA片段PCR扩增，并通过毛细管电泳进行分离和定量，通过将样品的峰图与正常(阴性)样本的峰图进行比较，确定样本中是否存在LGR。MLPA是一种高通量、廉价的检测LGR的方法。但其存在一些缺陷：(1)实验样本存在污染会降低检测的灵敏度；(2)无法检测单细胞水平的异常；(3)无法检测未知的点突变；(4)探针杂交区域附近存在单核苷酸多态性位点(Single Nucleotide Polymorphism，SNP)时，会影响探针的连接效率，进一步影响检测结果。

拷贝数变异(Copy Number Variation，CNV)一般指长度1Kbp以上的基因组大片段的重复或缺失。LGR与CNV在定义上相似，二者的区别在于检测分辨率不同。CNV检测主要关注基因组的大片段的重复或缺失，其最小单位为基因；而LGR检测关注单个基因在外显子水平的重复或缺失，其最小单位为基因的外显子。与CNV相比，LGR的分辨率更高。在NGS中，常规的CNV检测方法由于分辨率较低，无法适用于LGR的检测。因此，亟需开发了一种基于NGS技术的LGR检测方法，用于检测肿瘤样本中发生的LGR。

有鉴于此，提出本发明。

发明内容

本申请的第一目的在于提供一种基因大片段重排检测的数据处理方法，用于解决目前NGS检测方法分辨率低无法检测LGR的问题。

为实现上述目的，本申请采用的技术方案如下：

本申请提供一种基因大片段重排检测的数据处理方法，其特征在于，所述方法包括如下步骤：

步骤1)将目标区域(Bed)划分成区间(Bin)；

步骤2)所有区间测序深度数据预处理，包括标准化处理、误差校正；

步骤3)计算测序深度数据中位值(log2)和中位偏差(spread)；

步骤4)计算待测样本与阴性样本的测序深度比值(logratio_bin)及权重(weight)；

步骤5)确定基因重复/缺失区域：基于区间logratio，采用循环二元分割方法确定断点；并计算区域拷贝数(cn)：

cn＝2^{logratio_region}×2

若断点区域cn≥2.5，则确定为重复区域；若断点区域cn≤1.5，则确定为缺失区域。

在另一优选例中，所述步骤1)中，区间Bin的取值范围为100～200bp。

在另一优选例中，所述步骤3)中，中位值(log2)、中位偏差(spread)可采用双权中位值(biweight_location)、双权中位偏差(biweight_midvariance)：

x_i为某个区间的测序深度值。

在另一优选例中，所述步骤4)中：

logratio_bin＝log2_待测样本-log2_阴性样本

log2_阴性样本为若干阴性样本数据按照步骤1)、2)、3)计算所得。

在另一优选例中，所述步骤5)中，循环二元分割方法为对目标区域所有区间的测序深度比值以位置从小到大排序，计算最大T统计量，若T统计量对应的P值在统计学上显著，则为断点；优选的，所述最大T统计量的方法为：

在另一优选例中，所述步骤2)中，标准化处理步骤为，所有区间测序深度数据进行以2为底的对数转换，统计中位值，并将每个转换后数据减去中位值；优选的，所述步骤2)中，误差校正包括消除GC占比、序列重复、区域密度其中任一或组合。

在另一优选例中，所述方法还包括：

步骤6)注释基因外显子，确定重复/缺失区域所在基因或者外显子。

本申请还提供一种计算机装置/设备/系统，包括存储器、处理器及存储在存储器上的计算机程序，所述处理器执行所述计算机程序以实现上述任一所述方法的步骤。

本申请还提供一种计算机可读存储介质，其上存储有计算机程序/指令，该计算机程序/指令被处理器执行时实现上述任一所述方法的步骤。

本申请还提供一种计算机程序产品，包括计算机程序/指令，该计算机程序/指令被处理器执行时实现上述任一所述方法的步骤。

相比于现有技术，本申请的有益效果在于：

1)本发明方法有效的解决了MLPA方法中无法检测未知突变的缺陷，可对NGS目标捕获区域的任意基因进行LGR的检测；

2)本发明方法通过优化检测算法，提高了检测分辨率和灵敏度，有效的解决了NGS中常规的CNV检测方法分辨率低无法检测LGR的问题；

3)本发明方法同时考虑了目标区域的测序深度和SNP位点突变频率的变化，二者相互验证，进一步保证LGR检测结果的准确性。

附图说明

图1基于测序深度比值检测LGR；

图2将目标区域(Bed)划分为区间(Bin)；

图3Bin拷贝数+SNP突变频率分布散点图；

图4Sample 5NF1基因(阴性，NGS)；

图5Sample_5NF1基因(阴性，MLPA)；

图6Sample_5NF2基因Exon1-4杂合缺失(阳性，NGS)；

图7Sample_5NF2基因Exon1-4杂合缺失(阳性，MLPA)。

具体实施方式

下面将结合实施例对本申请的实施方案进行详细描述，但是本领域技术人员将会理解，下列实施例仅用于说明本申请，而不应视为限制本申请的范围。实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市场购买获得的常规产品。

以下术语或定义仅仅是为了帮助理解本发明而提供。这些定义不应被理解为具有小于本领域技术人员所理解的范围。

除非在下文中另有定义，本发明具体实施方式中所用的所有技术术语和科学术语的含义意图与本领域技术人员通常所理解的相同。虽然相信以下术语对于本领域技术人员很好理解，但仍然阐述以下定义以更好地解释本发明。

如本发明中所使用，术语“包括”、“包含”、“具有”、“含有”或“涉及”为包含性的(inclusive)或开放式的，且不排除其它未列举的元素或方法步骤。术语“由...组成”被认为是术语“包含”的优选实施方案。如果在下文中某一组被定义为包含至少一定数目的实施方案，这也应被理解为揭示了一个优选地仅由这些实施方案组成的组。

在提及单数形式名词时使用的不定冠词或定冠词例如“一个”或“一种”，“所述”，包括该名词的复数形式。

本发明中的术语“大约”、“大体”表示本领域技术人员能够理解的仍可保证论及特征的技术效果的准确度区间。该术语通常表示偏离指示数值的±10％，优选±5％。

此外，说明书和权利要求书中的术语第一、第二、第三、(a)、(b)、(c)以及诸如此类，是用于区分相似的元素，不是描述顺序或时间次序必须的。应理解，如此应用的术语在适当的环境下可互换，并且本发明描述的实施方案能以不同于本发明描述或举例说明的其它顺序实施。

本发明中的术语“核酸”或“核酸序列”指包含核糖核酸、脱氧核糖核酸或其类似物单元的任何分子、优选聚合分子。所述核酸可为单链的或双链的。单链核酸可为变性双链DNA的一条链的核酸。或者，单链核酸可为不来源于任何双链DNA的单链核酸。

实施例1发明构建

本发明核心思路是通过比较两个样本在基因组上目标区域的测序深度的差异，计算两个样本在每个目标区域的测序深度比值，确定存在重复或缺失的区域，以此获得彼此相对的拷贝数，如图1所示，具体试验设计及流程如下。

(1)划分区间

如图2所示，将目标区域(Bed)划分成长度相近的区间(Bin)。划分的Bin越长，检测的敏感性越低；划分的Bin越短，检测的特异性越低。优选的，所述Bin取值范围为100～200bp。

(2)构建基线

如表1所示，选择若干阴性(正常)样本使用NGS方法获得每个样本的每个Bin的测序深度。但每个Bin的测序深度会受到样本测序量、GC占比、序列重复和目标区域密度的影响。

表1阴性样本的每个Bin的测序深度

chromosome	start	end	depth
				chr1	11169315	11169445	4015.88
chr1	11169675	11169815	3469.62
				…	…	…	…
chr1	11174350	11174444	5589.18
				chr1	11174444	11174539	4837.35

为消除每个样本的测序量对Bin的测序深度的影响，针对每个样本，对所有Bin的测序深度进行标准化处理，具体步骤如下：

a.将每个Bin的测序深度进行以2为底的对数转换；

b.使用上一步的结果计算所有Bin的测序深度中位值；

c.将每个Bin经以2为底对数转换后的测序深度减去上一步得到的中位值。

为消除GC占比、序列重复和目标区域密度对测序深度的影响，针对每个样本，对标准化之后的测序深度进行误差校正平滑处理，以GC占比为例，具体步骤如下：

a.计算每个Bin的GC占比，并将Bin按照GC占比从小到大的顺序进行排序；

b.以滑动窗口取中位值为GC占比引起的误差(Bias)；

c.将每个Bin的测序深度减去每个Bin的Bias。

完成每个样本的深度标准化和误差校正后，计算每个Bin在所有正常样本中的测序深度的中位值和中位偏差，具体步骤如下：

a.计算测序深度的中位值；

M为中位值，MAD为绝对中位偏差。

b.计算测序深度的中位偏差；

M为中位值，MAD为绝对中位偏差。

基线如表2所示。Log2_阴性样本表示上述经过标准化、误差校正的每个Bin在若干阴性样本中的测序深度的中位值，spread表示上述经过标准化、误差校正的每个Bin在若干阴性样本中的测序深度的中位偏差。

表2若干阴性样本构建的基线

chromosome	start	end	Log2<sub>阴性样本</sub>	depth	gc	spread
							chr1	11169315	11169445	0.185291	4468.56	0.407692	0.0454296
chr1	11169675	11169815	0.089919	4250.19	0.442857	0.0461256
							…	…	…	…	…	…	…
chr1	11174350	11174444	0.216256	5375.45	0.531915	0.0319252
							chr1	11174444	11174539	-0.038949	4496.49	0.536842	0.0389313

(3)计算深度比值

基线构建完成后，计算待测样本与阴性样本的测序深度比值，具体步骤如下：

a.使用NGS方法获得每个待测样本Bin的测序深度，并使用与构建基线中相同的标准化方法对测序深度进行标准化；

b.使用与构建基线中相同的误差校正方法，校正GC占比、序列重复和目标区域密度对测序深度的影响；

c.使用Bin的长度、基线中Bin的测序深度中位值log2和中位偏差spread计算每个Bin的权重weight；

d.计算样本与阴性样本在每个Bin的测序深度比值；

logratio_bin＝log2_待测样本-log2_阴性样本

测序深度比值的结果如表3所示，logratio是样本与阴性样本在每个Bin上的测序深度比值，weight表示每个Bin的权重。

表3与阴性样本比较后的测序深度比值结果

chromosome	start	end	logratio_bin	depth	weight
						chr1	11169315	11169445	0.189527	4015.88	0.792957
chr1	11169675	11169815	-0.235901	3469.62	0.840212
						…	…	…	…	…	…
chr1	11174350	11174444	0.0116619	5589.18	0.712979
						chr1	11174444	11174539	0.0515672	4837.35	0.738904

(4)确定重复/缺失区域

获得样本与阴性样本在目标区域的测序深度比值后，使用循环二元分割(Circular Binary Segmentation，CBS)方法确定重复/缺失的区域。如图1所示，关键是找到重复/缺失区域与正常区域的断点(Change Points)，断点是为了在目标区域的中，找出Bin拷贝数差异较小的区域。但是区域内部Bin拷贝数差异较小并不代表该区域发生重复/缺失，只有该区域的拷贝数达到一定的阈值才能判断该区域是否发生重复/缺失。。

基本原理为：假设X₁，X₂，…X_n表示N个连续并且按照染色体、位置从小到大排序的Bin的测序深度比值，在找到下一个断点之前，如果存在X_v(1＜v＜n)，使得X₁，X₂，…X_v的分布F₀不同于X_v+1，X_v+2，…X_n的分布F₁，称X_v(1＜v＜n)为X₁，X₂，…X_n内的一个断点。

具体步骤如下：

a.将所有Bin按照参考基因组的染色体、位置从小到大进行排序；

b.对于目标区域的所有Bin，假设X₁，X₂，…X_n表示N个连续并且按照染色体、位置从小到大排序的Bin的测序深度比值，计算最大T统计量，若T统计量对应的P值在统计学上显著(例如0.01)，则认为X_i，X_j为两个断点。

c.计算每个区域所有Bin的加权平均测序深度比值和拷贝数，若该区域的拷贝数≥2.5，则表示该区域发生了重复，若该区域的拷贝数≤1.5，则表示该区域发生了缺失。

cn＝2^{logratio_region}×2

重复/缺失区域的结果如表4所示，cn表示每个区域的拷贝数，logratio_region表示每个区域的测序深度比值，probes表示每个区域包含的所有Bin的数量，weight表示每个区域的所有Bin的权重之和。

表4确定重复/缺失区域后的结果

chromosome	start	end	Cn	logratio_region	depth	probes	weight
								chr17	41196291	41229870	1.94	-0.043884	433.664	238	163.806
chr17	41230350	41250138	0.93	-1.11064	218.474	156	110.881
								chr17	41250977	41277601	2.02	0.011853	371.578	185	121.845

(5)注释基因外显子

将每个区域与参考基因组的注释文件比较，取交集，获得对应的基因和外显子。

表5重复/缺失区域的注释结果

(6)结果展示

如图4所示，将所有Bin和SNP按照染色体、位置从小到大排序，依次画散点图。横坐标表示染色体的位置，纵坐标分别表示拷贝数和突变频率，圆点分别表示Bin和SNP。

实施例2灵敏度与特异性检验

根据已发表综述中的LGR[11]，使用阴性样本模拟BRCA1和BRCA2基因不同外显子的缺失作为阳性样本，然后通过本分析方法对阴性样本和阳性样本进行数据分析，比较每个样本的数据检出结果是否与已知模拟结果一致，验证分析方法的灵敏度与特异性。统计结果如表6所示，灵敏度与特异性均为100％。

表6 LGR模拟数据分析结果

阳性样本	70	阴性样本	70
				检出阳性	70	检出阴性	70
灵敏度	70/70(100％)	特异性	70/70(100％)

实施例3临床样本的准确性检验

为验证NGS检测LGR的准确性，对临床样本分别使用NGS和MLPA方法检测LGR，比较检测结果的一致性。如表7所示，临床样本NGS与MLPA检测结果的一致性为100％。部分示例结果图如图4-7所示：Sample_5的NF1未发生大片段重排，NF2基因发生Exon1-4杂合缺失。图4和图6表示NGS的检测结果，图5和图7表示MLPA的检测结果，其横坐标表示基因的外显子，纵坐标表示拷贝数比值(待检样本拷贝数/阴性样本拷贝数)。

表7临床样本NGS与MLPA检测结果

Sample	Gene	NGS	MLPA
				Sample_1	NF1	阴性	阴性
Sample_1	NF2	阴性	阴性
				Sample_2	NF1	阴性	阴性
Sample_2	NF2	阴性	阴性
				Sample_3	NF1	阴性	阴性
Sample_3	NF2	阴性	阴性
				Sample_4	NF1	阴性	阴性
Sample_4	NF2	阴性	阴性
				Sample_5	NF1	阴性	阴性
Sample_5	NF2	Exon 1-4杂合缺失	Exon 1-4杂合缺失
				Sample_6	NF1	阴性	阴性
Sample_6	NF2	Exon 2-9杂合缺失	Exon 2-9杂合缺失
				Sample_7	NF1	阴性	阴性
Sample_7	NF2	阴性	阴性
				Sample_8	NF1	阴性	阴性
Sample_8	NF2	阴性	阴性
				Sample_9	BRCA1	Exon 8-12杂合缺失	Exon 8-12杂合缺失
Sample_10	BRCA2	Exon 2-3杂合缺失	Exon 2-3杂合缺失
				Sample_11	BRCA1	Exon 3杂合缺失	Exon 3杂合缺失
Sample_12	BRCA1	Exon 2-18杂合缺失	Exon 2-18杂合缺失

综上所述，本发明的技术要领及特点，其目的在于让熟知此技术的业内人士能够了解本发明的内容并能以此实施。本发明的内容并不局限在上述的实施例中，凡根据本发明技术思想实质所作的等效变化或修饰，都在本发明的保护范围之内。

Claims (7)

1.一种基因大片段重排检测的数据处理方法，其特征在于，所述方法包括如下步骤：

步骤1)将目标区域Bed划分成区间Bin；

步骤2)所有区间测序深度数据预处理，包括标准化处理、误差校正；

步骤3)计算测序深度数据中位值log2和中位偏差spread；

步骤4)计算待测样本与阴性样本的测序深度比值logratio_bin及权重weight；

步骤5)确定基因重复或缺失区域：基于区间logratio，采用循环二元分割方法确定断点，并计算区域拷贝数cn：

cn＝2^{IogLatio_region}×2

若断点区域cn≥2.5，则确定为重复区域；若断点区域cn≤1.5，则确定为缺失区域；

所述步骤1)中，区间Bin的取值范围为100～200bp；

所述步骤4)中：

logratio_bin＝log2_待测样本-log2_阴性样本

log2_阴性样本为若干阴性样本数据按照步骤1)、2)、3)计算所得；

所述步骤5)中，循环二元分割方法为对目标区域所有区间的测序深度比值以位置从小到大排序，计算最大T统计量，若T统计量对应的P值在统计学上显著，则为断点；所述最大T统计量的方法为：

2.如权利要求1所述方法，其特征在于，所述步骤3)中，中位值log2、中位偏差spread采用双权中位值biweight_location、双权中位偏差biweight_midvariance：

x_i为某个区间的测序深度值。

3.如权利要求1-2任一所述方法，其特征在于，所述步骤2)中，标准化处理步骤为，所有区间测序深度数据进行以2为底的对数转换，统计中位值，并将每个转换后数据减去中位值；所述步骤2)中，误差校正包括消除GC占比、序列重复、区域密度其中任一或组合。

4.如权利要求3所述方法，其特征在于，所述方法还包括：

步骤6)注释基因外显子，确定重复或缺失区域所在基因或者外显子。

5.一种计算机设备，其特征在于，包括存储器、处理器及存储在存储器上的计算机程序，所述处理器执行所述计算机程序以实现权利要求1-4任一所述方法的步骤。

6.一种计算机可读存储介质，其上存储有计算机程序或指令，其特征在于，该计算机程序或指令被处理器执行时实现权利要求1-4任一所述方法的步骤。

7.一种计算机程序产品，包括计算机程序或指令，其特征在于，该计算机程序或指令被处理器执行时实现权利要求1-4任一所述方法的步骤。

Images