CN113436678A - 一种基于滤波降噪的基因组结构变异检测方法 - Google Patents

Abstract

本发明提出了在一种基于滤波降噪的基因组拷贝数变异检测方法，该方法考虑了染色体序列本身的特性以及序列中碱基GC含量对读深的影响，在更好的对数据进行预处理的基础上，将读深数据与高斯核函数进行卷积运算得到尺度空间图像函数，并对尺度空间图像进行边缘检测及基准检测，确定候选拷贝数区域，及检测拷贝数变异类型和位置，提高了拷贝数变异检测的精度。

Description

一种基于滤波降噪的基因组结构变异检测方法

技术领域

本发明涉及生物信息学领域，具体涉及一种基于滤波降噪的基因组拷贝数变异检测方法。

背景技术

新一代测序(New generation sequencing，NGS)技术的发展越来越成熟，各测序平台层出不穷，基因序列的测序成本大幅度地下降，测序的速度越来越高，这使得测序产生的DNA序列数据非常庞大，数据处理的准确程度成为当务之急。

小波变换是信号处理当中去除数据噪声的有利利器，而尺度空间滤波是形象直观的描述信号的重要手段。生物信息学领域的测序数据，由于生物序列本身的高重复性以及测序过程中不可避免的误差，往往在对数据建模的时候，需要考虑噪声对模型带来的影响。在拷贝数变异的检测中，由于人类染色体序列的重复序列片段的存在会导致测序发生不可避免的错误，为了用合适的尺度描述拷贝数变异，采用了多尺度的尺度空间滤波，找出读深信号中的零交叉点，最终确定拷贝数变异区域。

另外，伴随着人类基因组计划及1000genomes project的实施与发展，蛋白质、DNA、RNA的序列数据的规模日趋增加，仅仅依靠生物实验来研究生物基因变异及疾病产生早已不能满足现实需要，因此必须借助计算机、数学等学科的理论及思想方法从海量数据中来研究和阐明生物学问题。拷贝数变异检测是生物信息学中研究生物基因结构改变的有效方法之一。

目前应用于拷贝数变异检测的技术主要有：

1.比较基因组杂交(CGH)：该技术发展至今，已与芯片技术(Microarray)结合后衍生为芯片比较基因组杂交技术(Array-CGH)。该技术可以在全部染色体或染色体亚带水平上，对不同基因组之间DNA序列的拷贝数进行检测，从而发现拷贝数变异。然而该技术分辨率在Mb水平，更小片段的拷贝数片段则不易检出。同时该技术操作繁琐，通量低、耗时长且成本昂贵，需要较为大量的模板DNA，不利于大范围的推广。

2.MLPA：全称为多重连接探针扩增技术，是2002年发展起来的一种拷贝数检测方法。目前已有相应的试剂盒检测如SMA、唐氏综合征等疾病。该技术具有较准确的相对定量功能。但是该方法探针制备较为复杂，同时操作步骤繁琐，耗时长。并且采用毛细管电泳作为分析手段，通量较低、成本较高且属于开放式操作，易于造成PCR产物的污染。

发明内容

在本发明中，提出了一种基于滤波降噪的基因组拷贝数变异检测方法，该方法考虑了染色体序列本身的特性以及序列中碱基GC含量对读深的影响，在更好的对数据进行预处理的基础上，将读深数据与高斯核函数进行卷积运算得到尺度空间图像函数，并对尺度空间图像进行边缘检测及基准检测，确定候选拷贝数区域，及检测拷贝数变异类型和位置，提高了拷贝数变异检测的精度。

具体包括如下步骤：

S1.数据预处理；

利用SAMtools工具从bam文件中提取出读深信号，读深信号由以下两部分信号组成的:R_m＝r_m+E_m，R_m代表观察到的读深信号的实际值，r_m代表在染色体序列期望得到的读深信号，E_m代表噪声信号；

采用haar函数

进行噪声的去除:

采用GC校正去除碱基GC含量对读深信号的影响；

S2.获得尺度空间图像；

将读深数据r[i-j]与高斯核函数K(x，y)进行卷积运算得到尺度空间图像函数I_SS[i，l]:

其中，

σ_l代表第l层的尺度参数，m代表高斯核函数K(i，j)的窗口值大小；

S3.尺度空间图像边缘检测；

将尺度空间图像函数I_SS[i，l]在不同的尺度x,y下求得其分量，分量为

和

在每一尺度下，求每个像素点的模值MI_SS[i，l]和相角AI_SS[i，l]：

MI_SS[i，l]在相角AI_SS[i，l]上取得极大值的点对应着尺度空间图像的突变点，由MI_SS[i，l和AI_SS[i，l]可以求得极值点，从而对尺度空间图像进行边缘检测；

S4.尺度空间图像基准检测；

设置三个基准标准m_t(l)、m_t(l)+λδ_t(l)和m_t(l)-λδ_t(l)，其中m_t(l)和δ_t(l)是尺度空间图像函数y_SS[i，l]在有两个非零的零交叉点函数Z_SS[i，l]在第l层的均值和标准差，λ是基准校验系数，尺度空间值的正常范围为m(k)±2δ(k)，在所述正常范围之外的尺度空间函数值将被滤除；

S5.确定候选拷贝数区域；

若Z_SS[s_m，l,l]·Z_SS[e_m，l,l]<0，第l层中间区{is_m，l≤i≤e_m，l}中所有的点满足Z_SS[i，l]＝0，且区间{i|s_m，l≤i≤e_m，l}上尺度空间图像函数I_SS[i，l]的均值

在m_t(l)+λδ_t(l)和m_t(l)-λδ_t(l)之间，则[s_m，l,e_m，l]是一个候选拷贝数变异的区域；其中i在零交叉点函数Z_SS[i，l]在第l层上的对应位置区间内；

S6.拷贝数变异类型和位置检测；

尺度空间图像函数均值

在m_t(l)+λδ_t(l)之上，则拷贝数变异增加；尺度空间图像函数均值

在m_t(l)-λδ_t(l)下，则拷贝数变异缺失。

进一步地，步骤S1中，噪声信号E_m的数学期望为E(X_K)＝μ,(k＝1,2,……)，方差为D(X_K)＝σ²≠0，(k＝1,2,……)；其中μ是随机变量的期望值，σ²是方差。

进一步地，步骤S1中，

向右或向左平移单位k后的函数变为

其中k为整数。

进一步地，步骤S3中，对各尺度的边缘图像设置阈值TH，求取边缘点；大于等于阈值TH的点作为边缘点保留，小于阈值TH的点置零；求取尺度空间图像边缘，将不同尺度下的边缘点链接起来，得到不同尺度下的尺度空间图像边缘。

进一步地，步骤S6中，当在第l层的第m个区间[s_m，l,e_m，l]发现拷贝数变异时，从第l层开始往下搜索，找到l-1层的第m个区间[s_m，l-1,e_m，l-1]，如此循环，直到迭代到原始的读深信号为止。

进一步地，步骤S2中，σ_l代表第l层的尺度参数，随着σ的增大，尺度空间图像越来越平滑，并且在平滑的过程中保持其轮廓不变。

进一步地，步骤S2中，m代表高斯核函数K(i，j)的窗口值大小，m＝3σ_l。

本文使用的术语解释如下：

拷贝数变异(CNV)：是指与正常样品中的相应核酸序列相比，待测试样品中核酸分子的至少一部分的拷贝数变化，其中所述部分具有大于1kb的长度。拷贝数变异的情况和原因可包括：缺失，诸如微缺失；插入，例如微插入、微复制、复制；倒位、转座和复杂的多位点变异。

测序：是指获得样品的核酸序列的信息的过程。测序可以通过各种方法进行，包括但不限于双脱氧链终止；优选地，高通量测序方法，包括但不限于下一代测序技术或单分子测序技术。测序深度越高，检测的灵敏度越高，即可以检测的缺失片段和重复片段的长度越小。

读深信号：是指具有一定长度(通常长于20bp)的核酸序列，例如由测序仪产生的序列的测序结果，其可以通过序列比对方法与参考序列的特定区域或位置比对。

索引：指具有特定长度并发挥标记功能的核酸序列。当待测试的DNA分子衍生自多个待测试的样品时，多个样品中的每一个可以添加有不同的索引，用于在测序期间区分多个样品。

GC含量偏差：批次之间或一个批次内存在一定的GC偏差，这可能导致拷贝数偏差呈现在基因组的具有高GC含量或低GC含量的区域中。用基于对照集的测序数据进行CG校正，以获得每个窗口中的校正的相对读段数，由此可以消除这种偏差，并且可以提高检测拷贝数变异的准确性。

附图说明

附图1为本发明的基于滤波降噪的基因组拷贝数变异检测方法的流程图。

附图2为变异检测方法进过数据预处理后的信号对比示例图。

附图3为变异检测方法中零交叉点和平滑后信号的拐点示例图。

具体实施方式

以下实施例结合附图对本发明作进一步的说明，所给出的是本发明的一些具体实施例，这些实施例只是说明而不表示本发明所有的可能性，本发明并不局限于这些实施例中提到的材料、反应条件或参数，任何在相关领域具备经验的人，都可以按照本发明的原理，利用其它类似材料或反应条件实现本发明所描述的基因拷贝数变异检测。这些并不脱离本发明描述的基本概念。

本文使用的词语“包括”、“包含”、“具有”或其任何其他变体意欲涵盖非排它性的包括。例如，包括列出要素的工艺、方法、物品或设备不必受限于那些要素，而是可以包括其他没有明确列出或属于这种工艺、方法、物品或设备固有的要素。

参考附图1为本发明的基于滤波降噪的基因组拷贝数变异检测方法的流程图：该变异检测方法具体包括以下步骤：

S1、数据预处理

读深信号是从bam文件中通过工具SAMtools提取得到的，bam文件中存储的是测序读段匹配到参考序列中的信息。利用SAMtools工具从bam文件中提取出read count文件，文件包括read counts值和对应的位置信息，我们将读深信号看作是由以下两部分信号组成的:

R_m＝r_m+E_m (1)

在式(1)中，其中R_m代表观察到的读深信号的实际值，r_m代表在染色体序列期望得到的读深信号，E_m代表噪声信号，一般都被认为是高斯白噪声信号。白噪声是在无限宽的频率范围内，功率分布均匀的噪声，只是一种理想化的噪声模型。高斯白噪声其幅度的统计规律服从高斯分布，而定义中的“白”是指它的功率谱在整个频域内为常数。噪声信号E_m的数学期望为E(X_K)＝μ,(k＝1,2,……)，方差为D(X_K)＝σ²≠0，(k＝1,2,……)。其中μ是随机变量的期望(或均值)，σ²是它的方差。

我们可以使用基于信号的方法找出读深信号当中的断点。小波理论在去噪和检测信号的断点信息方面都很有用处。小波去噪是从含噪信号中寻找到小波函数空间的最佳映射对信号进行滤波处理分析能够将高频干扰信号滤掉同时能成功地保留信号的原始特征将得到的特征信号与低通滤波后的信号进行组合重构。

我们采用haar函数

进行噪声的去除:

R_m代表观察到的读深信号的实际值:

对于整数k，

的图形为

向右或向左平移单位k后的结果。

参考图2为进过数据预处理后的信号对比示例图。我们之所以采用haar函数进行噪声去除主要是因为它与读深数据的结构吻合的很自然，沿着染色体的拷贝数变异都是以块存在的，而且被标记的相邻的基因座具有相同的拷贝数增加或减少。

在进行完小波去噪后，由于碱基GC含量对读深信号的影响，我们需要采用GC校正，GC校正采用现有技术中的GC校正方法，这里就不再详细的说明了。

S2、获得尺度空间图像

将读深数据r[i-j]与高斯核函数K(i，j)进行卷积运算得到尺度空间图像函数I_SS[i，l]:

σ_l代表第l层的尺度参数，随着σ的增大，尺度空间图像将变得越来越平滑，并且在平滑的过程中保持其轮廓不变。m代表高斯核函数K(i，j)的窗口值大小，默认为m＝3σ_l，σ_l的变化范围决定了能检测到的拷贝数变异区间长度的大小。两个相邻尺度参数之间的比例决定了时间复杂度和能检测到的拷贝数变异区间的精度。若采用较小的比例，则是以时间为代价换取高的检测精度，反之，若采用较大的比例，则是牺牲检测精度来换取较低的时间复杂度。所以选择一个合适的比例既能获得较高的检测精度又能获得较小的时间复杂度。

S3、尺度空间图像边缘检测

尺度空间函数滤波过程是利用一个平滑函数，在不同的尺度下平滑所要检测的图像信号，根据平滑后信号的小波变换系数模的一阶或二阶导数找出信号的突变点。一阶导数的极值点对应二阶导数的零交叉点和平滑后信号的拐点，参考图2。因此可由小波变换模局部极大值检测图像边缘。

将尺度空间图像函数I_SS[i，l]在不同的尺度x,y下求得其分量，分量为

和

在每一尺度下，求每个像素点的模值MI_SS[i，l]和相角AI_SS[i，l]：

MI_SS[i，l]在相角AI_SS[i，l]上取得极大值的点对应着尺度空间图像的突变点，由MI_SS[i，l和AI_SS[i，l]可以求得极值点，从而对尺度空间图像进行边缘检测。

优选地，可以对各尺度的边缘图像设置阈值TH，求取边缘点，大于等于TH的点作为边缘点保留，小于TH的点置零；求取图像边缘，将不同尺度下的边缘点链接起来，得到了不同尺度下的图像边缘。

S4、尺度空间图形基准检测

为了进一步去除离群值对检测精度的影响，很有必要对尺度空间图形进行基准检测。这里我们设置三个基准标准，分别为m_t(l)、m_t(l)+λδ_t(l)和m_t(l)-λδ_t(l)。其中m_t(l)和δ_t(l是尺度空间图像函数I_SS[i，l]在有两个非零的零交叉点函数Z_SS[i，l]在第l层的均值和标准差,λ是基准校验系数，优选为3。为了滤除离群点，尺度空间值的正常范围为m(k)±2δ(k)，在所述正常范围之外的尺度空间函数值将被滤除；

S5、确定候选拷贝数区域

当将每层的零交叉点都找到后，对于拷贝数变异的候选区域就可以从有两个非零值的每层的零交叉点着手，对于区间[s_m，l,e_m，l]表示在第l层的第m个区间，在区域{i|s_m，l≤i≤e_m，l}，其中i在零交叉点函数Z_SS[i，l]在第l层上的对应位置区间内。

若Z_SS[s_m，l,l]·Z_SS[e_m，l,l]<0，第l层中间区{i|s_m，l≤i≤e_m，l}中所有的点满足Z_SS[i，l]＝0，且区间{i|s_m，l≤i≤e_m，l}上尺度空间图像函数I_SS[i，l]的均值

在m_t(l)+λδ_t(l)和m_t(l)-λδ_t(l)之间，则[s_m，l,e_m，l]是一个候选拷贝数变异的区域；其中i在零交叉点函数Z_SS[i，l]在第l层上的对应位置区间内；

S6、拷贝数变异检测

拷贝数变异检测包括检测拷贝数变异的类型(增加和缺失)以及检测拷贝数变异区域的精准位置。拷贝数变异增加定义为：尺度空间函数的均值

在m_t(l)+λδ_t(l)之上，拷贝数变异缺失(LOSS)定义为：尺度空间函数均值

在m_t(l)-λδ_t(l)之下。

当在第l层的第m个区间[s_m，l,e_m，l]发现拷贝数变异时，从第l层开始往下搜索，找到l-1层的第m个区间[s_m，l-1,e_m，l-1]，如此循环，直到迭代到原始的读深信号为止。

本发明另一方面，可提供一种计算机可读存储介质，所述可读存储介质上存储有机器可执行指令，所述机器可执行指令在被执行时使机器执行根据本发明所述的基于滤波降噪的基因组结构变异检测方法的步骤。

在本发明中，计算机可读存储介质可以是可以保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如包括但不限于，电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：便携式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。这里所使用的计算机可读存储介质不被解释为瞬时信号本身，诸如无线电波或者其他自由传播的电磁波、通过波导或其他传输媒介传播的电磁波(例如，通过光纤电缆的光脉冲)、或者通过电线传输的电信号。

这里所描述的机器可执行指令可以从计算机可读存储介质下载到各个计算/处理设备，或者通过网络、例如因特网、局域网、广域网和/或无线网下载到外部计算机或外部存储设备。网络可以包括铜传输电缆、光纤传输、无线传输、路由器、防火墙、交换机、网关计算机和/或边缘服务器。每个计算/处理设备中的网络适配卡或者网络接口从网络接收机器可执行指令，并转发该机器可执行指令，以供存储在各个计算/处理设备中的计算机可读存储介质中。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (7)

1.一种基于滤波降噪的基因组拷贝数变异检测方法，其特征在于，包括如下步骤：

S1.数据预处理；

利用SAMtools工具从bam文件中提取出读深信号，读深信号由以下两部分信号组成的：R_m＝r_m+E_m，R_m代表观察到的读深信号的实际值，r_m代表在染色体序列期望得到的读深信号，E_m代表噪声信号；

采用haar函数

进行噪声的去除：

采用GC校正去除碱基GC含量对读深信号的影响；

S2.获得尺度空间图像；

将读深数据r[i-j]与高斯核函数K(i，j)进行卷积运算得到尺度空间图像函数I_SS[i，l]：

其中，

σ_l代表第l层的尺度参数，m代表高斯核函数K(i，j)的窗口值大小；

S3.尺度空间图像边缘检测；

将尺度空间图像函数I_SS[i，l]在不同的尺度x，y下求得其分量，分量为

和

在每一尺度下，求每个像素点的模值MI_SS[i，l]和相角AI_SS[i，l]：

MI_SS[i，l]在相角AI_SS[i，l]上取得极大值的点对应着尺度空间图像的突变点，由MI_SS[i，l]和AI_SS[i，l]可以求得极值点，从而对尺度空间图像进行边缘检测；

S4.尺度空间图形基准检测；

设置三个基准标准m_t(l)、m_t(l)+λδ_t(l)和m_t(l)-λδ_t(l)，其中m_t(l)和δ_t(l)是尺度空间图像函数I_SS[i，l]在有两个非零的零交叉点函数Z_SS[i，l]在第l层的均值和标准差，λ是基准校验系数，尺度空间函数值的正常范围为m(k)±2δ(k)，在所述正常范围之外的尺度空间函数值将被滤除；

S5.确定候选拷贝数区域；

若Z_SS[s_m，l，l]·Z_SS[e_m，l，l]＜0，第l层中间区{i|s_m，l≤i≤e_m，l}中所有的点满足Z_SS[i，l]＝0，且区间{i|s_m，l≤i≤e_m，l}上尺度空间图像函数I_SS[i，l]的均值

在m_t(l)+λδ_t(l)和m_t(l)-λδ_t(l)之间，则[s_m，l，e_m，l]是一个候选拷贝数变异的区域；其中i在零交叉点函数Z_SS[i，l]在第l层上的对应位置区间内；

S6.拷贝数变异检测；

尺度空间图像函数均值

在m_t(l)+λδ_t(l)之上，则拷贝数变异增加；尺度空间图像函数均值

在m_t(l)-λδ_t(l)之下，则拷贝数变异缺失。

2.根据权利要求1所述的基于滤波降噪的基因组拷贝数变异检测方法，其特征在于：步骤S1中，噪声信号E_m的数学期望为E(X_K)＝μ，(k＝1，2，……)，方差为D(X_K)＝σ²≠0，(k＝1，2，……)；其中μ是随机变量的期望值，σ²是方差。

3.根据权利要求1所述的基于滤波降噪的基因组拷贝数变异检测方法，其特征在于：步骤S1中，

向右或向左平移单位k后的函数变为

其中k为整数。

4.根据权利要求1所述的基于滤波降噪的基因组拷贝数变异检测方法，其特征在于：步骤S3中，对各尺度的边缘图像设置阈值TH，求取边缘点；大于等于阈值TH的点作为边缘点保留，小于阈值TH的点置零；求取尺度空间图像边缘，将不同尺度下的边缘点链接起来，得到不同尺度下的尺度空间图像边缘。

5.根据权利要求1所述的基于滤波降噪的基因组拷贝数变异检测方法，其特征在于：步骤S6中，当在第l层的第m个区间[s_m，l，e_m，1]发现拷贝数变异时，从第l层开始往下搜索，找到l-1层的第m个区间[s_m，l-1，e_m，l-1]，如此循环，直到迭代到原始的读深信号为止。

6.根据权利要求1所述的基于滤波降噪的基因组拷贝数变异检测方法，其特征在于：步骤S2中，σ_l代表第l层的尺度参数，随着σ的增大，尺度空间图像越来越平滑，并且在平滑的过程中保持其轮廓不变。

7.根据权利要求1所述的基于滤波降噪的基因组拷贝数变异检测方法，其特征在于：步骤S2中，m代表高斯核函数K(i，j)的窗口值大小，m＝3σ_l。

Images