CN108920899A - 一种基于目标区域测序的单个外显子拷贝数变异预测方法 - Google Patents

Abstract

本发明涉及一种基于目标区域测序的单个外显子拷贝数变异预测方法，包括和测序数据处理和拷贝数变异预测2个步骤，其中拷贝数变异预测步骤包括统计覆盖到目标区域的测序序列总数目和碱基总数目，确定对照外显子区域，标准化每个对照样本和实验样本的待分析外显子区域覆盖度，计算对照样本中待分析外显子区域标准化后的覆盖度的平均值、标准差和变异系数，预测待分析外显子区域拷贝数变化步骤。本发明不用基于全基因组测序，直接利用外显子水平的覆盖度信息进行分析外显子水平的拷贝数变异，分析方法简单，不经过复杂的GC校正和建模。

Description

一种基于目标区域测序的单个外显子拷贝数变异预测方法

技术领域

本发明涉及生物医学领域，具体涉及一种基于目标区域测序的单个外显子拷贝数变异预测方法。

背景技术

从2003年人类基因组计划结束起，基因组测序技术突飞猛进，基于高通量测序技术的单核苷酸多态性(SNP)检测技术，已经成熟和普及。高通量测序技术可以实现数以万计的DNA分子同时进行边合成边测序反应，极大的提高的测序通量。基因检测的测序成本下降速度甚至快于计算机领域的摩尔定律。基于测序技术的应用，生物学的研究从传统的单基因、单位点的研究进入了组学研究时代，由此产生了一系列的研究成果和和具有社会价值的临床应用。

拷贝数变异是结构变异的一种，研究已经证实拷贝数变异与人类疾病的发生相关，例如智力缺陷、自闭症、精神分裂、癌症发生等。不同于基因的单碱基变异，外显子水平的拷贝数变异(Exon Copy Number Variant)是一类不常见但非常重要的突变类型，约10％的BRCA1癌症是由外显子拷贝数变异引起的。典型的外显子拷贝数变异可能导致蛋白紊乱甚至丧失功能。过往检测拷贝数变异的方法是利用多重连接探针扩增技术(multiplexligation-dependent probe amplification，MLPA)，染色体芯片(CMA)或荧光PCR。随着测序技术的发展，利用NGS的数据分析基因拷贝数变化是越来越受到关注和有效的确定基因拷贝数的方法。

拷贝数变异分析一般基于全基因组测序，但成本较目标区域测序高，并且分辨率较低，只能得到大片段的拷贝数变异(通常大于1Mbp)，不能用于检测外显子水平的拷贝数变异(外显子拷贝数变异的序列长度在100bp左右)。

利用目标区域测序数据进行分析是一种目的性强(有明确的分析目标基因)、成本节约的方法。但目标区域测序存在各区域捕获效率不一致等问题，覆盖度均一性较全基因组差。随着技术的发展，针对外显子水平的分析工具开始出现，多数软件都利用了外显子区域的覆盖度信息，用参考基因组GC含量校正，然后根据不同算法来识别拷贝数变异。这些工具存在这分析流程复杂，需要的对照样本总量在30个以上或者需要配对的对照样本。

发明内容

本发明所要解决的技术问题在于提供一种基于目标区域测序的单个外显子拷贝数变异预测方法，不用基于全基因组测序，直接利用外显子水平的覆盖度信息进行分析外显子水平的拷贝数变异，分析方法简单，不经过复杂的GC校正和建模。

为解决上述技术问题，本发明提供的一种基于目标区域测序的单个外显子拷贝数变异预测方法，包括以下步骤：

(1)测序数据处理：

a.对样本的目标区域进行捕获建库，然后在二代测序仪中进行高通量测序，得到测序原始数据；

b.通过开源比对工具，对测序得到的序列与参考基因组进行比对，得到测序序列在参考基因组上的位置和序列比对质量；

c.根据序列在染色体上的位置，对位于相同起始位置和相同终止位置的序列只保留一条，并对序列按照染色体顺序和起始位置顺序进行排序；

(2)拷贝数变异预测：

a.统计覆盖到目标区域的测序序列总数目和碱基总数目：

以目标区域为基础，以每个区域为单位，统计每个覆盖到该区域的测序序列总数目和碱基总数目；

b.确定对照外显子区域：

每个样本选取至少30个以上的包含在目标区域中且为不容易发生拷贝数变化的区域作为对照外显子区域，每个样本的对照外显子区域与待分析外显子区域一致；

c.标准化每个对照样本和实验样本的待分析外显子区域覆盖度；

d.计算对照样本中待分析外显子区域标准化后的覆盖度的平均值、标准差和变异系数；

e.根据变异系数、实验样本待分析外显子区域标准化后的覆盖度偏离对照样本标准化后的覆盖度平均值与标准差的关系，预测待分析外显子区域拷贝数变化。

优选地，所述步骤(2)c步骤中标准化单个样本的待分析外显子区域的覆盖度的具体方法如下：

对于特定样本，其对照外显子区域为m个，待分析外显子区域为n个，第j个对照外显子区域定义为ControlRegion_j,j∈{1,2,3…m}，其覆盖的reads数目定义为region_reads_j；第i个待分析外显子区域定义为TestRegion_i,i∈{1,2,3…n}，其覆盖的reads数目定义为region_reads_i；

该样本覆盖到对照外显子区域中的总reads数为：Total_reads_control＝∑region_reads_j,j∈{1,2,3...m}，该样本待分析外显子区域的覆盖度标准化基数数值等于对照外显子区域的总reads数Total_reads_control除以对照区域总个数m：该样本的待分析外显子区域TestRegion_i的标准化之后的覆盖度为：

对照样本数量为p个，第c个实验样本定义为smpc,t∈{1,2,3,…p}，根据上述方法标准化其每个待分析外显子区域的覆盖度，记为smp_{c_i}；

实验样本数量为q个，第t个实验样本定义为smpt,t∈{1,2,3,…q}，根据上述方法标准化其每个待分析外显子区域的覆盖度，记为smp_{t_i}；

所述步骤(2)d步骤中计算对照样本中待分析外显子区域标准化后的覆盖度的平均值、标准差和变异系数的具体方法如下：

对于p个对照样本，待分析外显子区域TestRegion_i,i∈{1,2,3...n}的标准化之后的覆盖度的平均值为：

标准差为：

变异系数为：

所述步骤(2)e步骤中根据变异系数、实验样本待分析外显子区域标准化后的覆盖度偏离对照样本标准化后的覆盖度平均值与标准差的关系，预测待分析外显子区域拷贝数变化的具体方法如下：

如果CV_i＞0.2，这提示该区域变异较大，预测不可靠；如果CV_i≤0.2，假设外显子区域的覆盖度在样本之间的符合正态分布，则正态分布的均值为标准差为σ_i；则对于特定实验样本smpt,t∈{1,2,3,...q}，其待分析外显子区域数值偏离对照样本覆盖度平均值为：

如果|diff _{smpt_i}|≤2σ_i，则认为没有发生拷贝数变化；如果|diff _{smpt_i}|＞2σ_i，diff smp_{t_i}分大于0和小于0两种情况：

第一种情况diff smp_{t_i}＞0，可能为拷贝数增加情况，则预测为重复变异，并根据正态分布的密度函数计算样本落入smp_{t_i}至+∞区间内的概率；第二种情况diff smp_{t_i}＜0,可能为拷贝数减少情况，则预测为缺失变异，并根据正态分布的密度函数计算样本落入-∞至区smpt_i间内的概率。

优选地，所述测序原始数据默认为通用的FASTQ格式，FASTQ文件中的原始序列经过去接头、去除低质量碱基的质控方法，得到用于比对的高质量序列数据。

优选地，所述目标区域为BED文件格式，至少记录目标区域所在的染色体、所在染色体的目标区域起始位置和目标区域染色体终止位置、单个区域的长度及区域编号或者所在基因等信息。

优选地，所述开源比对工具为BWA，bowtie2。

本发明的有益效果为：

1、利用目标区域测序数据分析外显子水平的拷贝数变异，而非基于全基因组测序，成本低。

2、直接利用外显子水平的覆盖度信息进行分析，而非大片段拷贝数变异，不经过复杂的GC校正和建模，减少外显子水平的拷贝数变异分析复杂性。

3、不需要配对样本，仅需要较少的对照样本，分析实验样本的外显子拷贝数变异。

4、利用分析脚本而非软件进行分析，方式简单，没有很复杂的分析流程。

5、利用已有数据，不需要经过专门的拷贝数变异检测实验，成本低。

附图说明

下面结合附图和实施方式对本发明作进一步详细的说明。

图1为本发明的拷贝数变异预测流程图。

图2为本发明测序数据处理的流程图。

具体实施方式

基于目标区域测序的单个外显子拷贝数变异预测方法，包括以下步骤：

(1)测序数据处理(其流程如图2所示)：

a.对样本的目标区域进行捕获建库，然后在二代测序仪中进行高通量测序，得到测序原始数据。测序原始数据默认为通用的FASTQ格式，FASTQ文件中的原始序列经过去接头、去除低质量碱基的质控方法，得到用于比对的高质量序列数据。捕获目标区域的探针可以自行设计或者用现成的试剂盒。质控软件可以使用FASTQC或trimmomatic。以trimmomatic为例：

java-jar$trimmonmatic PE-threads 8\

$fq1$fq2\

$out/$smp.clean.R1.fastq.gz$out/$smp.unpaired_R1.fastq.gz\

$out/$smp.clean.R2.fastq.gz$out/$smp.unpaired_R2.fastq.gz\

ILLUMINACLIP:$adapter:2:30:10LEADING:3TRAILING:3SLIDINGWINDOW:4:15MINLEN:36

b.通过开源比对工具，对测序得到的序列与参考基因组进行比对，得到测序序列在参考基因组上的位置和序列比对质量。比对工具可以使用BWA，bowtie2。以BWA为例：

$bwa mem-M-t 8-R"@RG\tID:$smp\tSM:$smp\tLB:$smp\tPL:ILLUMINA"$reffa$fq1$fq2>${dir}/${smp}.sam

c.为了保证结果的准确性，根据序列在染色体上的位置，对位于相同起始位置和相同终止位置的序列只保留一条，为了便于文件的快速检索，对序列按照染色体顺序和起始位置顺序进行排序。可以使用Picard软件进行序列去重和排序：

java-jar-Djava.io.tmpdir＝$HOME/tmp$picard SortSam\

INPUT＝${dir}/${smp}.sam\

OUTPUT＝${dir}/${smp}.sorted_reads.bam\

SORT_ORDER＝coordinate

java-jar-Djava.io.tmpdir＝$HOME/tmp$picard MarkDuplicates\

INPUT＝${dir}/${smp}.sorted_reads.bam\

OUTPUT＝${dir}/${smp}.dedup_reads.bam\

METRICS_FILE＝${dir}/${smp}.dedup_reads.bam.metrics.txt

java-jar$picard BuildBamIndex INPUT＝${dir}/${smp}.dedup_reads.bam

(2)拷贝数变异预测(其流程如图1所示)：

a.统计覆盖到目标区域的测序序列总数目和碱基总数目：

以目标区域为基础，以每个区域为单位，统计每个覆盖到该区域的测序序列总数目和碱基总数目。目标区域，为BED文件格式，至少记录目标区域所在的染色体，所在染色体的目标区域起始位置和目标区域染色体终止位置，单个区域的长度、区域编号或者所在基因等信息。该步骤利用软件为Bedtools。

bedtools coverage-b sample.dedup_reads.bam-a target.bed

其输出的格式为八列，如：

chr 31190464 31190530 NAME 52 66 66 1.000

第一列表示染色体编号

第二列表示区域的起始位置

第三列表示区域的终止位置

第四列为该区域的名称

第五列表示覆盖到该区域的测序序列数目

第六列表示该区域被测序序列覆盖的碱基数目，

第七列表示该区域的长度，

第八列表示该区域被测序序列覆盖的碱基占区域总碱基的比例。

b.确定对照外显子区域：

每个样本选取至少30个以上的包含在目标区域中且为不容易发生拷贝数变化的区域作为对照外显子区域，不论对照样本还是待分析样本，每个样本的对照外显子区域与待分析外显子区域一致；

c.标准化每个对照样本和实验样本的待分析外显子区域覆盖度：

对于特定样本，其对照外显子区域为m个，待分析外显子区域为n个，第j个对照外显子区域定义为ControlRegion_j,j∈{1,2,3…m}，其覆盖的reads数目定义为region_reads_j；第i个待分析外显子区域定义为TestRegioni,i∈{1,2,3…n}，其覆盖的reads数目定义为region_reads_i；

该样本覆盖到对照外显子区域中的总reads数为：Total_readscontrol＝∑region_reads_j,j∈{1,2,3...m}，该样本待分析外显子区域的覆盖度标准化基数数值等于对照外显子区域的总reads数Total_reads_control除以对照区域总个数m：该样本的待分析外显子区域TestRegion_i的标准化之后的覆盖度为：

对照样本数量为p个，第c个实验样本定义为smp_c,t∈{1,2,3,…p}，根据上述方法标准化其每个待分析外显子区域的覆盖度，记为smp_{c_i}；

实验样本数量为q个，第t个实验样本定义为smp_t,t∈{1,2,3,…q}，根据上述方法标准化其每个待分析外显子区域的覆盖度，记为smp_{t_i}；

d.计算对照样本中待分析外显子区域标准化后的覆盖度的平均值、标准差和变异系数；

对于p个对照样本，待分析外显子区域region_i的标准化之后的覆盖度的平均值为：

标准差为：

变异系数为：

e.根据变异系数、实验样本待分析外显子区域标准化后的覆盖度偏离对照样本标准化后的覆盖度平均值与标准差的关系，预测待分析外显子区域拷贝数变化：

如果CV_i＞0.2，这提示该区域变异较大，预测不可靠；

如果CV_i≤0.2，假设外显子区域的覆盖度在样本之间的符合正态分布，则正态分布的均值为标准差为σ_i；则对于特定实验样本smp_t,t∈{1,2,3,...q}，其待分析外显子区域数值偏离对照样本覆盖度平均值为：

如果|diff _{smpt_i}|≤2σ_i，则认为没有发生拷贝数变化；

如果|diff _{smpt_i}|＞2σ_i，diff _{smpt_i}分大于0和小于0两种情况：

第一种情况diff _{smpt_i}＞0，可能为拷贝数增加情况，则预测为重复变异，标为DUP，并根据正态分布的密度函数计算样本落入smp_{t_i}至+∞区间内的概率；第二种情况diffsmp_{t_i}＜0，可能为拷贝数减少情况，则预测为缺失变异，标为DEL，并根据正态分布的密度函数计算样本落入-∞至区smp_{t_i}间内的概率。

具体实施例1

对已知有外显子水平拷贝数变异的三个阳性样本进行分析，三个阳性样本的外显子拷贝数变异信息如下。

将三个阳性样本和五个阴性对照样本进行外显子组测序，得到测序数据。对测序数据进行质控，比对到hg19参考基因组上，用picard对比对后的reads进行去重和排序。用到的软件为trommomatic、bwa、picard。各样本的统计信息如下：

样本	Raw Bases	Duplication	Q20	Q30
					NA05123	18389070300	27.87％	97.58％	93.50％
NA09981	17933438100	27.00％	97.52％	93.42％
					NA23159	18144067200	26.40％	97.47％	93.30％
NA05169	23433262500	33.61％	98.22％	95.06％
					Control1	18983356500	19.88％	98.7％	97.8％
Control2	10140261600	19.15％	95.87％	93.52％
					Control3	23808953400	24.85％	98.53％	96.17％
Control4	53167072800	34.17％	98.39％	95.80％
					Control5	24256914900	24.89％	98.54％	96.18％

使用软件bedtools统计每个外显子的覆盖度，然后标准化每个样本的待分析外显子的覆盖度。将五个对照样本的覆盖度信息组成对照样本组，逐个分析阳性样本。检测结果如下：

样本	阳性外显子拷贝数变异数目	检测结果
			NA05123	18个外显子拷贝数变异	16
NA09981	3个外显子拷贝数变异	1
			NA23159	1个外显子拷贝数变异	1
NA05169	40个外显子拷贝数变异	40

检测的总灵敏度为(16+1+1)/(18+3+1)＝81.82％。对于没有检测出的外显子拷贝数变异，通过查看了该外显子的测序覆盖度发现该区域在目标区域捕获测序技术下，覆盖度较差，所以没有检测出来。

以样本NA23159为例，该样本的第7号外显子发生了重复变异。该区域在NA23159和质控样本的覆盖度为：

基于正态分布，预测拷贝数变异的统计数值：

项目	数值	备注
			NA23159在7号外显子上的覆盖度	0.958
质控样本在7号外显子上的覆盖度平均值	0.724
			质控样本在7号外显子上的覆盖度方差	0.077
质控样本在7号外显子上的覆盖度变异系数	0.106	小于20％
			NA23159在7号外显子上的覆盖度偏离值	0.234	大于三倍方差
NA23159在7号外显子上为该覆盖度的概率	0.001
			该外显子拷贝数变异预测结果	DUP	发生重复(Duplication)突变

结果表明，NA23159在该外显子区域的覆盖度值偏离对照样本均值较大，算法预测为拷贝数增加即重复变异。

具体实施例2

对于缺失变异的预测，我们同样取得了很好的结果，因为发生缺失的样本，在其外显子上的测序序列覆盖度近乎为0。以NA05169样本为例，该样本有40个外显子发生了拷贝数缺失变异，预测结果如下：

本方法将NA05169样本的40个缺失变异全部检测出。

综上所述，本方法对拷贝数重复和缺失变异的预测，没有用GC含量进行校正，也没有像其它软件一样进行复杂的建模进行预测，仅利用了5个对照样本，即将该重复和缺失变异预测出，显示了算法良好的应用性能。数据来源为已有的目标区域测序得到的数据，不需要增加额外的实验成本。

上述说明是示例性的而非限制性的。通过上述说明本领域技术人员可以意识到本发明的许多种改变和变形，其也将落在本发明的实质和范围之内。

Claims (5)

1.一种基于目标区域测序的单个外显子拷贝数变异预测方法，其特征在于：包括以下步骤：

(1)测序数据处理：

a.对样本的目标区域进行捕获建库，然后在二代测序仪中进行高通量测序，得到测序原始数据；

b.通过开源比对工具，对测序得到的序列与参考基因组进行比对，得到测序序列在参考基因组上的位置和序列比对质量；

c.根据序列在染色体上的位置，对位于相同起始位置和相同终止位置的序列只保留一条，并对序列按照染色体顺序和起始位置顺序进行排序；

(2)拷贝数变异预测：

a.统计覆盖到目标区域的测序序列总数目和碱基总数目：

以目标区域为基础，以每个区域为单位，统计每个覆盖到该区域的测序序列总数目和碱基总数目；

b.确定对照外显子区域：

每个样本选取至少30个以上的包含在目标区域中且为不容易发生拷贝数变化的区域作为对照外显子区域，每个样本的对照外显子区域与待分析外显子区域一致；

c.标准化每个对照样本和实验样本的待分析外显子区域覆盖度；

d.计算对照样本中待分析外显子区域标准化后的覆盖度的平均值、标准差和变异系数；

e.根据变异系数、实验样本待分析外显子区域标准化后的覆盖度偏离对照样本标准化后的覆盖度平均值与标准差的关系，预测待分析外显子区域拷贝数变化。

2.根据权利要求1所述的基于目标区域测序的单个外显子拷贝数变异预测方法，其特征在于：

所述步骤(2)c步骤中标准化单个样本的待分析外显子区域的覆盖度的具体方法如下：

对于特定样本，其对照外显子区域为m个，待分析外显子区域为n个，第j个对照外显子区域定义为ControlRegion_j,j∈{1,2,3…m}，其覆盖的reads数目定义为region_reads_j；第i个待分析外显子区域定义为TestRegion_i,i∈{1,2,3…n}，其覆盖的reads数目定义为region_reads_i；

该样本覆盖到对照外显子区域中的总reads数为：Total_reads_control＝∑region_readsj,j∈{1,2,3...m}，该样本待分析外显子区域的覆盖度标准化基数数值等于对照外显子区域的总reads数Total_readscontrol除以对照区域总个数m：该样本的待分析外显子区域TestRegioni的标准化之后的覆盖度为：

对照样本数量为p个，第c个实验样本定义smpc,c∈{1,2,3,…p}，根据上述方法标准化其每个待分析外显子区域的覆盖度，记为smp_{c_i}；

实验样本数量为q个，第t个实验样本定义为smpt,t∈{1,2,3,…q}，根据上述方法标准化其每个待分析外显子区域的覆盖度，记为smp_{t_i}；

所述步骤(2)d步骤中计算对照样本中待分析外显子区域标准化后的覆盖度的平均值、标准差和变异系数的具体方法如下：

对于p个对照样本，待分析外显子区域TestRegion_i,i∈{1,2,3...n}的标准化之后的覆盖度的平均值为：

标准差为：

变异系数为：

所述步骤(2)e步骤中根据变异系数、实验样本待分析外显子区域标准化后的覆盖度偏离对照样本标准化后的覆盖度平均值与标准差的关系，预测待分析外显子区域拷贝数变化的具体方法如下：

如果CV_i＞0.2，这提示该区域变异较大，预测不可靠；如果CV_i≤0.2，假设外显子区域的覆盖度在样本之间的符合正态分布，则正态分布的均值为标准差为σ_i；则对于特定实验样本smpt,t∈{1,2,3,...q}，其待分析外显子区域数值偏离对照样本覆盖度平均值为：

如果|diff smp_{t_i}|≤2σ_i，则认为没有发生拷贝数变化；

如果|diff smp_{t_i}|＞2σ_i，diff smp_{t_i}分大于0和小于0两种情况：

第一种情况diff smp_{t_i}＞0，可能为拷贝数增加情况，则预测为重复变异，并根据正态分布的密度函数计算样本落入smp_{t_i}至+∞区间内的概率；第二种情况diff smp_{t_i}＜0,可能为拷贝数减少情况，则预测为缺失变异，并根据正态分布的密度函数计算样本落入-∞至区smp_{t_i}间内的概率。

3.根据权利要求1所述的基于目标区域测序的单个外显子拷贝数变异预测方法，其特征在于：所述测序原始数据默认为通用的FASTQ格式，FASTQ文件中的原始序列经过去接头、去除低质量碱基的质控方法，得到用于比对的高质量序列数据。

4.根据权利要求1所述的基于目标区域测序的单个外显子拷贝数变异预测方法，其特征在于：所述目标区域为BED文件格式，至少记录目标区域所在的染色体、所在染色体的目标区域起始位置和目标区域染色体终止位置、单个区域的长度及区域编号或者所在基因等信息。

5.根据权利要求1所述的基于目标区域测序的单个外显子拷贝数变异预测方法，其特征在于：所述开源比对工具为BWA，bowtie2。