阅读说明：本技术 一种任意阶分段多项式信号的处理方法 (Method for processing arbitrary-order segmented polynomial signals ) 是由 段君博 王青 王玉平 于 2020-06-22 设计创作，主要内容包括：本发明涉及一种任意阶分段多项式信号的处理方法,以解决现有分段多项式信号分割效果差的问题。本发明提供一种任意阶分段多项式信号的处理方法：首先,从原子力显微镜或从高通量基因组测序仪中提取信号；然后按照函数模型对信号进行分割拟合,再通过计算机算法程序进行动态规划和矩阵分解,根据输出分割值位置,绘制分割拟合后的新信号；最后,根据新信号检测断点位置,分割出伸直长度及持续长度估计所需要的合适区域,完成蛋白质解折叠,或分割出杨氏模量估计所需要的合适区域,或检测拷贝数变异区间及变异类型。(The invention relates to a method for processing a piecewise polynomial signal of any order, which aims to solve the problem of poor segmentation effect of the traditional piecewise polynomial signal. The invention provides a method for processing a piecewise polynomial signal of any order, which comprises the following steps: firstly, extracting signals from an atomic force microscope or a high-throughput genome sequencer; then, carrying out segmentation fitting on the signals according to a function model, carrying out dynamic planning and matrix decomposition through a computer algorithm program, and drawing new signals after segmentation fitting according to the positions of output segmentation values; and finally, according to the position of the breakpoint of the new signal detection, segmenting a proper region required by extension length and continuous length estimation to finish protein unfolding, or segmenting a proper region required by Young modulus estimation, or detecting a copy number variation interval and a variation type.)

一种任意阶分段多项式信号的处理方法

技术领域

本发明涉及一种任意阶分段多项式信号分割拟合的方法。

背景技术

在工程实践与科学实验中，常常需要对试验数据进行拟合，其中多项式曲线拟合是一种较常用的数据拟合方法。当数据点较多时，多项式阶数太低，拟合精度和效果不太理想，要提高拟合精度和效果就需要提高曲线阶数，但阶数太高又带来计算上的复杂性及其他方面的不利。

现有高阶分段多项式信号分析出现在许多科学领域，如原子力显微镜(atomicforce microscopy，AFM)对蛋白质折叠的数据分析、对杨氏模量的测量，以及新一代测序技术(Next generation sequencing，NGS)对拷贝数变异(copy number variation，CNV)的检测等，而信号分析关键则是对信号进行分割，即断点检测，一旦确定了分割，就可以分别对各个片段进行分析(曲线拟合、参数估计等)。现有信号分割及拟合常常采用经典的分割算法循环二分分割(circular binary segmentation，CBS)，但该分割算法分割效果差、不能快速找到断点且对噪声敏感。

发明内容

本发明目的在于提供了一种任意阶分段多项式信号的处理方法，用以克服现有分段多项式信号分割效果差的问题。

本发明所采用的技术方案为：一种任意阶分段多项式信号的处理方法，包括以下步骤：

1)信号采集

获取原子力显微镜输出的原始信号y₀，并进行归一化处理，得到信号y；或者获取高通量基因组测序仪输出的短读的fasta格式文件，并从中提取信号y；

2)对步骤1)得到的信号y，按如下函数模型进行分割拟合：

所述的函数模型为

其中：

v为信号分割位置；

ε(v_k-1+1，v_k)是信号y中第k段的拟合误差；

λ是给定参数的惩罚值，取任意正数，用于调节分割拟合的质量；

K为信号的分割段数；

3)向量和矩阵的初始化

优化值向量Φ(1)＝[0]；

分割值位置向量q{1}＝[]；

向量S＝[1]；

矩阵

矩阵Ω(1)＝[G₀]；

向量e(1)＝[0]；

4)循环执行步骤A1至步骤A7，循环标志i＝1:(N-P)，共循环(N-P)次，其中N为信号y的长度，P为多项式阶数，

A1)按下式对向量S中的每个元素t进行重新赋值，得到向量中的元素：

A2)根据重新赋值后的向量

求矩阵

的最小值，其中j为向量S的标志位置；

A3)根据矩阵

的最小值，对信号y重新分割，并将新的分割位置存入分割值位置向量q{i+1}；

A4)选出同时满足条件①和条件②的标志位置j，保留满足标志位置j条件的矩阵Ω、矩阵e、矩阵B和矩阵S中的相应位置元素，其余元素删除；其中：

条件①为：1≤j≤s，s为向量S中的元素个数

条件②为：

A5)根据所保留的元素，按下式分别计算A、Г和ρ三个变量；

A＝[α₁，α₂，...，α_s]

其中：

li＝i-S(j)+2+P；

A为P+1行s列的矩阵，每列均具有α的形式；

l是子信号的长度，每列的l的计算公式为i-S(j)+2+P；

A6)根据A、Г和ρ三个变量，按下式分别更新矩阵B、矩阵Ω和向量e；

B＝B+y_v+1A

e＝e+ρ⊙(y_v+11_s-(Γ^T⊙B^T)1_P+1)²

A7)得到更新后的矩阵B为矩阵Ω为[Ω，G₀]以及向量e为[e，O]，并将循环标志i***到步骤A4)中保留的向量S，使得向量S为[S，i+1]；

5)Φ向量中第(N+1-P)位作为输出优化值，q向量中第(N+1-P)位再加上数值(P-1)后作为输出分割值位置，所述输出分割值位置即为信号分割位置υ的实际分割位置；根据输出分割值位置，绘制分割拟合后的新信号y₁；

6)根据新信号y₁找到两个断点位置，分割出伸直长度及持续长度估计所需要的合适区域；

或者，根据新信号y₁找到两个断点位置，分割出杨氏模量估计所需要的合适区域；

或者，根据新信号y₁，检测拷贝数变异区间及变异类型。

进一步地，步骤1)中原子力显微镜输出的原始信号y₀为一组力曲线y₀(z)，其中y₀是探针与样品之间的相互作用力，z是偏移距离；步骤(1)中信号y是通过对原子力显微镜输出的原始信号y₀归一化值进行线性内插得到的。

进一步地，步骤1)中，从短读的fasta格式文件中提取信号y的具体方式是：首先通过比对软件从短读的fasta格式文件得到SAM格式文件或者压缩的BAM格式文件，再利用计算程序从SAM格式文件或者压缩的BAM格式文件得到读深信号y。所述比对软件为MAQ软件或bowtie软件，所述计算程序为samtools程序；

进一步地，在对力曲线y₀(z)进行处理时，步骤4)中的多项式阶数P取2或者3。

进一步地，在对高通量基因组测序仪的数据检测拷贝数变异时，步骤4)中的多项式阶数P取0。

进一步地，本发明可处理0到任意阶分段多项式信号，对高阶分段多项式信号处理更佳。同时，本发明又可降低原子力显微镜及高通量测序数据分析的成本。

本发明与现有技术相比具有以下有益效果。

一、本发明采用的一种任意阶分段多项式信号的处理方法，基于L-0范数惩罚的最小二乘稀疏模型，最小二乘法提高了拟合的保真度，同时L-0范数惩罚优选了分割的个数，通过动态规划及矩阵分解的方法，得到了最优分割(即最优解)，加快了运算速度，缩短了运算时间。

二、本发明采用的一种任意阶分段多项式信号的处理方法，可以实现任意阶的多项式拟合，即P可选任意非负整数，例如在原子力显微镜的力曲线分析中P的取值为2或3，在分析新一代测序技术数据检测拷贝数变异中P的取值为0，应用范围广泛。

三、本发明采用的一种任意阶分段多项式信号的处理方法，在原子力显微镜蛋白质解折叠的数据分析和杨氏模量的测量中，分割效果好，提高了检测断点的精确性，降低了数据分析的成本；在使用基于新一代测序技术的检测技术，提升了拷贝数变异的检测准确度。

附图说明

图1为本发明动态规划和矩阵分解的算法程序图。

图2为本发明实施例1中蛋白质解折叠的信号图。

图3为本发明实施例1中蛋白质解折叠的分割拟合处理结果图，其中多个与纵轴平行的线段为蛋白质解折叠的分割位置。

图4为本发明实施例1中蛋白质解折叠的信号与分割拟合处理结果的对比分析图，其中，横坐标为探针偏移距离，纵坐标为探针与样品之间的相互作用力。

图5为实施例1现有软件Fordis的数据分析图。

图6本发明为实施例2中测量杨氏模量的酵母菌压力测试信号图。

图7为本发明实施例2中测量杨氏模量的分割拟合处理结果图，其中两个星号位置即为分割位置。

图8为实施例2测量杨氏模量中信号拟合误差结果图。

图9为实施例2测量杨氏模量的酵母菌压力测试信号、分割拟合处理结果以及信号拟合误差(放大十倍)结果的对比分析图，其中，横坐标为探针偏移距离，纵坐标为探针与样品之间的相互作用力。

图10为本发明实施例3中基于新一代测序技术的拷贝数变异检测的信号图。

图11为实施例3基于循环二分分割(CBS)算法的拷贝数变异的检测结果图。

图12为本发明实施例3中基于新一代测序技术的拷贝数变异的检测结果图。

图13为实施例3检测拷贝数变异中图10、图11和图12的对比分析图，其中横坐标为基因座，纵坐标为读深。

具体实施方式

下面将结合本发明的实施例和附图，对本发明的技术方案进行清楚、完整地描述。显然，所描述的实施例并非对本发明的限制。

实施例1

本发明可对蛋白质解折叠的数据进行分析，使用原子力显微镜采集的非洲爪蟾***的环状核苷酸负通道亚单位α1(CNGA1)，原子力显微镜的力曲线包含多个断点，相邻的两个断点间可以用虫状连接链(worm like chain，WLC)模型或自由连接链(freelyjoint chain，FJC)模型拟合，通过拟合这些模型可以估计持续长度(persistence length)及伸直长度(contour length)等。

本实施例为二阶分段多项式信号的处理方法，包括以下步骤：

(1)信号采集

获取原子力显微镜输出的原始信号y₀，并进行归一化处理，如图2所示，得到信号y；

(2)对步骤1)得到的信号y，按如下函数模型进行分割拟合：

所述的函数模型为

其中：

v为信号分割位置；

ε(v_k-1+1，v_k)是信号y中第k段的拟合误差；

λ是给定参数的惩罚值，取任意正数，用于调节分割拟合的质量；

K为信号的分割段数；

(3)调用上述函数模型，通过图1所示的计算机算法程序进行动态规划和矩阵分解，输出Φ向量优化值，以及输出q向量分割值位置即为实际分割位置υ；并根据输出实际分割位置υ绘制分割拟合后的新信号y₁，对新信号y₁进行分析。

在计算机算法程序中：

(8)为：

(9)为：B＝B+y_v+1A

(10)为：

(11)为：

(12)为：

(13)为：e＝e+ρ⊙(y_v+11_s-(Γ^T⊙B^T)1_P+1)²

(14)为：

如图3和图4所示，使用了二阶(P＝2)分段多项式信号进行分割拟合，可对蛋白质解折叠的数据分析，可以检测到力曲线y₁在探针偏移距离30nm处的断点，本发明分割效果更好。

如图5所示的现有软件Fordis的数据分析结果，则无法检测到探针偏移距离30nm处的断点。

实施例2

使用原子力显微镜的力曲线也可以对杨氏模量进行测量。图6为测量杨氏模量的酵母菌压力测试信号图，图7为测量杨氏模量的分割拟合处理结果图。其中两个星号位置即为分割位置，这段力曲线在探针偏移距离600nm之前的星号左边，由于探针未和样品接触，探针与样品之间的相互作用力为0；而在力曲线在探针偏移距离600nm之后的第二个星号右边，由于压力过大，为负载饱和的线性区域。只有在力较小的范围内(两个星号之间)为杨氏模型合适的拟合区域。使用传统的方法很难检测到第二个星号的位置。

使用本发明可以自动、快速地找到两个星号断点，从而分割出杨氏模量估计所需要的合适区域。该力曲线测量质量较高，噪声很小，进而该力曲线拟合误差很小。图8为实施例2测量杨氏模量中信号拟合误差结果图，即图6中测量杨氏模量的酵母菌压力测试信号图和图7中测量杨氏模量的分割拟合处理结果图之间拟合误差的结果图。

图9为实施例2测量杨氏模量的酵母菌压力测试信号、分割拟合处理结果以及信号拟合误差(放大十倍)结果的对比分析图，其中，横坐标为探针偏移距离，纵坐标为探针与样品之间的相互作用力，可以看出使用本发明的信号处理方法对分段多项式新信号分割和拟合的效果最佳。

实施例3

本实施例为零阶分段多项式信号的处理方法，应用在新一代测序技术中，包括以下步骤：

1)信号采集

获取高通量基因组测序仪输出的短读的fasta格式文件，并从中提取信号y，如图10所示。

从短读的fasta格式文件中提取信号y的具体方式为：首先通过比对软件从短读的fasta格式文件得到SAM格式文件或者压缩的BAM格式文件，再利用计算程序从SAM格式文件或者压缩的BAM格式文件得到读深信号y。

2)对步骤1)得到的信号y，按如下函数模型进行分割拟合：

所述的函数模型为

其中：

v为信号分割位置；

ε(v_k-1+1，v_k)是信号y中第k段的拟合误差；

λ是给定参数的惩罚值，取任意正数，用于调节分割拟合的质量；

K为信号的分割段数；

2)调用上述函数模型，通过图1所示的计算机算法程序进行动态规划和矩阵分解，输出Φ向量优化值，以及输出q向量分割值位置即为实际分割位置υ；并根据输出实际分割位置υ绘制分割拟合后的新信号y₁，对新信号y₁进行分析。

在计算机算法程序中：

(8)为：

(9)为：B＝B+y_v+1A

(10)为：

(11)为：

(12)为：

(13)为：e＝e+ρ⊙(y_v+11_s-(Γ^T⊙B^T)1_P+1)²

(14)为：

本实施例使用了零阶(P＝0)分段多项式信号进行分割拟合，从新一代测序技术数据中检测拷贝数变异。作为对照，如图11所示，使用了经典的分割算法循环二分分割，可以看到，循环二分分割在基因座约920处有虚假的检测结果。并且在同一台计算机上，循环二分分割是35秒；如图12和图13所示，本发明所有的拷贝数变异都可以检测到，本发明的计算时间是0.5秒，本发明的计算时间高效快速。

以上所述仅为本发明的实施例，并非对本发明保护范围的限制，凡是利用本发明说明书及附图内容所作的等效结构变换，或直接或间接运用在其他相关的技术领域，均包括在本发明的专利保护范围内。