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

本发明专利技术涉及一种任意阶分段多项式信号的处理方法，以解决现有分段多项式信号分割效果差的问题。本发明专利技术提供一种任意阶分段多项式信号的处理方法：首先，从原子力显微镜或从高通量基因组测序仪中提取信号；然后按照函数模型对信号进行分割拟合，再通过计算机算法程序进行动态规划和矩阵分解，根据输出分割值位置，绘制分割拟合后的新信号；最后，根据新信号检测断点位置，分割出伸直长度及持续长度估计所需要的合适区域，完成蛋白质解折叠，或分割出杨氏模量估计所需要的合适区域，或检测拷贝数变异区间及变异类型。

【技术实现步骤摘要】
一种任意阶分段多项式信号的处理方法
本专利技术涉及一种任意阶分段多项式信号分割拟合的方法。
技术介绍
在工程实践与科学实验中，常常需要对试验数据进行拟合，其中多项式曲线拟合是一种较常用的数据拟合方法。当数据点较多时，多项式阶数太低，拟合精度和效果不太理想，要提高拟合精度和效果就需要提高曲线阶数，但阶数太高又带来计算上的复杂性及其他方面的不利。现有高阶分段多项式信号分析出现在许多科学领域，如原子力显微镜(atomicforcemicroscopy，AFM)对蛋白质折叠的数据分析、对杨氏模量的测量，以及新一代测序技术(Nextgenerationsequencing，NGS)对拷贝数变异(copynumbervariation，CNV)的检测等，而信号分析关键则是对信号进行分割，即断点检测，一旦确定了分割，就可以分别对各个片段进行分析(曲线拟合、参数估计等)。现有信号分割及拟合常常采用经典的分割算法循环二分分割(circularbinarysegmentation，CBS)，但该分割算法分割效果差、不能快速找到断点且对噪声敏感。
技术实现思路
本专利技术目的在于提供了一种任意阶分段多项式信号的处理方法，用以克服现有分段多项式信号分割效果差的问题。本专利技术所采用的技术方案为：一种任意阶分段多项式信号的处理方法，包括以下步骤：1)信号采集获取原子力显微镜输出的原始信号y0，并进行归一化处理，得到信号y；或者获取高通量基因组测序仪输出的短读的fasta格式文件，并从中提取信号y；2)对步骤1)得到的信号y，按如下函数模型进行分割拟合：所述的函数模型为其中：v为信号分割位置；ε(vk-1+1，vk)是信号y中第k段的拟合误差；λ是给定参数的惩罚值，取任意正数，用于调节分割拟合的质量；K为信号的分割段数；3)向量和矩阵的初始化优化值向量Φ(1)＝[0]；分割值位置向量q{1}＝[]；向量S＝[1]；矩阵矩阵Ω(1)＝[G0]；向量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＝[α1，α2，...，αs]其中：li＝i-S(j)+2+P；A为P+1行s列的矩阵，每列均具有α的形式；l是子信号的长度，每列的l的计算公式为i-S(j)+2+P；A6)根据A、Г和ρ三个变量，按下式分别更新矩阵B、矩阵Ω和向量e；B＝B+yv+1Ae＝e+ρ⊙(yv+11s-(ΓT⊙BT)1P+1)2A7)得到更新后的矩阵B为矩阵Ω为[Ω，G0]以及向量e为[e，O]，并将循环标志i插入到步骤A4)中保留的向量S，使得向量S为[S，i+1]；5)Φ向量中第(N+1-P)位作为输出优化值，q向量中第(N+1-P)位再加上数值(P-1)后作为输出分割值位置，所述输出分割值位置即为信号分割位置υ的实际分割位置；根据输出分割值位置，绘制分割拟合后的新信号y1；6)根据新信号y1找到两个断点位置，分割出伸直长度及持续长度估计所需要的合适区域；或者，根据新信号y1找到两个断点位置，分割出杨氏模量估计所需要的合适区域；或者，根据新信号y1，检测拷贝数变异区间及变异类型。进一步地，步骤1)中原子力显微镜输出的原始信号y0为一组力曲线y0(z)，其中y0是探针与样品之间的相互作用力，z是偏移距离；步骤(1)中信号y是通过对原子力显微镜输出的原始信号y0归一化值进行线性内插得到的。进一步地，步骤1)中，从短读的fasta格式文件中提取信号y的具体方式是：首先通过比对软件从短读的fasta格式文件得到SAM格式文件或者压缩的BAM格式文件，再利用计算程序从SAM格式文件或者压缩的BAM格式文件得到读深信号y。所述比对软件为MAQ软件或bowtie软件，所述计算程序为samtools程序；进一步地，在对力曲线y0(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测量杨氏模量的酵母菌压力测试信号、分割拟合处理结果以及信号拟合误差(放大十倍)结果的对比分析图，其中，横坐标为探针偏移距离，纵坐标为探针与样品本文档来自技高网...

【技术保护点】
1.一种任意阶分段多项式信号的处理方法，其特征在于，包括以下步骤：/n1)信号采集/n获取原子力显微镜输出的原始信号y

【技术特征摘要】
1.一种任意阶分段多项式信号的处理方法，其特征在于，包括以下步骤：
1)信号采集
获取原子力显微镜输出的原始信号y0，并进行归一化处理，得到信号y；或者获取高通量基因组测序仪输出的短读的fasta格式文件，并从中提取信号y；
2)对步骤1)得到的信号y，按如下函数模型进行分割拟合：
所述的函数模型为



其中：
v为信号分割位置；
ε(vk-1+1，vk)是信号y中第k段的拟合误差；
λ是给定参数的惩罚值，取任意正数，用于调节分割拟合的质量；
K为信号的分割段数；
3)向量和矩阵的初始化
优化值向量Φ(1)＝[0]；
分割值位置向量q{1}＝[]；
向量S＝[1]；
矩阵
矩阵Ω(1)＝[G0]；
向量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＝[α1，α2，...，αs]






其中：






A为P+1行s列的矩阵，每列均具有α的形式；
l是子信号的长度，每列的l的计算公式为i-S(j)+2+P；
A6)根据A、Г和ρ三个变量，按下式分别更新矩阵B、矩阵Ω和向量e；
B＝B+yυ+1A



e＝e+ρ⊙(yv+11s-(Γ...

【专利技术属性】
技术研发人员：段君博，王青，王玉平，
申请(专利权)人：西安交通大学，
类型：发明
国别省市：陕西;61