一种面向多信号的快速断点检测方法、系统、设备及存储介质与流程

1.本发明属于信号断点检测技术领域，具体涉及一种面向多信号的快速断点检测方法、系统、设备及存储介质。


背景技术：

2.传统的断点检测方法主要包含循环二分分割(circular binary segmentation，cbs)，最优分割(optimal partitioning，op)，删减的精确线性时间(pruned exact linear time，pelt)等方法。
3.cbs是针对基因芯片数据开发的拷贝数变异(copy number variation，cnv)检测方法，目前也广泛应用于基于高通量测序(high
‑
throughput sequencing，hts)数据的拷贝数变异检测，r语言中dnacopy就是基于cbs算法实现的。假设染色体上只有一个断点，且数据服从正态分布，则对断点两侧的数据进行双样本的t检验，就可以判断断点的有无。通过遍历检验该染色体所有可能的位点，就可以确定这个断点的位置。如果染色体上存在多个断点，则首先确定第一个断点的位置，之后分别对断点两侧的染色体再做同样的处理，就可以确定第二及第三个断点的位置；如此反复处理，即实现所谓的二分分割。若信号不服从正态分布，则可以使用其他合适的检验代替t检验。因此，cbs通过反复计算最大对数似然比，使用假设检验的方法将数值相近的区域连接起来，进而完成了信号二分分割。
4.cbs方法不能从全局角度证明分割的最优性，而op方法则可以证明其全局最优性。op方法利用动态规划(dynamic programming，dp)思想，将长信号的分割问题分解为若干短信号分割的子问题，并通过合并子问题的解，得到原问题的解。op从长度为1信号的分割问题开始，通过归纳法及分解子问题的方法，可以求解任意长度信号的分割问题。
5.虽然op方法具有全局最优性，但由于op在应用动态规划的过程中，需要遍历测试所有可能情况，其计算量随信号长度呈平方增长。对于大规模问题(如上万点长度的信号)其计算量惊人，很难满足实际应用需求。然而可以证明，在某些特定条件下，一些情况不具备可能性，因此在遍历过程中可删除这些情况，进而避免不必要的计算。因此pelt方法突破了op在计算量方面的瓶颈，其计算量降低到与信号长度呈线性增长关系，使得应用范围极大得到扩展，但只能检测单一信号的断点。
6.许多的科学研究、工程应用中需要检测信号中的断点。这里所谓的断点指的是信号中的一个位置，信号在其两边呈现不同的模式(数学上严格地讲是分布)。例如图1中包含两个高低台阶的信号，信号在台阶内保持同一模式(均值不变)，而在台阶间会有模式的改变。如果能检测到断点(粗线)，则断点间信号的均值容易求得。
7.上述cbs、op、pelt等现有方法均只能检测单一信号中的断点，而实际应用中往往需要检测多个信号中共有的断点(如图2所示)。如基于高通量测序技术的拷贝数变异检测问题，由于读深信号包含噪声，因此单个信号的真阳性率(true positive rate，tpr)低，假阳性率(false positive rate，fpr)高。提高检测性能的直接方法是增加测序覆盖度，然而这会导致实验成本的提高。替代方案是对样品进行多次中等或低覆盖度的测序，或使用多
个平台测序，即多重测序。多重测序可以减少由单个样品或平台引入的系统误差，能够提高检测性能，但需要能够检测多个信号断点的技术。另外，一个群体中的多个个体会共享cnv，而许多复杂疾病也可能共享cnv，因此需要从多样本的角度出发，在多个信号中检测出共有的cnv。
8.因此，上述cbs、op、pelt等现有方法均只能检测单一信号中的断点，不能检测多个信号中共有的断点，而且检测计算量大导致实验成本高。


技术实现要素：

9.为了克服上述现有技术的缺点，本发明的目的在于提供一种面向多信号的快速断点检测方法、系统、设备及存储介质，旨在解决现有技术中断点检测方法只能检测单一信号中的断点，检测效率低的缺陷性技术问题。
10.本发明提出的一种面向多信号的快速断点检测方法，包括如下步骤：
11.s1、对原始信号进行预处理，得到预处理后的n
×
m尺寸的预处理信号矩阵y；
12.s2、根据实际情况确定断点数k或惩罚参数λ；
13.s3、当选择惩罚参数λ作为输入参数时，求解最小化优化问题，获取处理后的n
×
m尺寸的信号矩阵x；当选择断点数k作为输入参数时，则根据预处理信号矩阵y计算最大可能λ
max
，在区间[0,λ
max
]内搜索估计惩罚参数λ，并求解给定λ下的最小化优化问题，使得断点数为k，获取处理后的n
×
m尺寸的信号矩阵x；
[0014]
s4、根据求得的断点位置计算信号矩阵x；对信号矩阵x进行去归一化处理，得到处理后的信号x0，实现分段信号的断点快速检测。
[0015]
优选地，在s1具体包括如下步骤：
[0016]
s1.1、将采集得到的原始信号存储入n
×
m的原始矩阵y0；
[0017]
s1.2、计算原始矩阵y0的最大绝对值c；
[0018]
s1.3、采用最大绝对值c对原始信号进行预处理操作。
[0019]
优选地，在s1.3中，采用归一化法对原始信号进行预处理，预处理结果如公式(1)所示：
[0020]
y＝y0/c
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)。
[0021]
优选地，在s4中，处理后的信号x0的计算公式如公式(2)所示：
[0022]
x0＝cx
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)。
[0023]
优选地，在s3中，采用二分法搜寻区间[0,λ
max
]内的惩罚参数λ。
[0024]
优选地，在s3中，在给定预处理信号矩阵y及惩罚参数λ的情况下，求解如公式(3)所示的最小化优化问题，即求解信号矩阵x：
[0025][0026]
其中，y是n
×
m的预处理信号矩阵，x是处理后的n
×
m尺寸的信号矩阵，n是每个信号的采样点数，m是信号的个数，λ是对每个断点的惩罚参数，p(x)是x中的断点个数。
[0027]
优选地，s3的具体操作步骤如下：
[0028]
当选择断点数k作为输入参数时：
[0029]
1)、根据预处理信号矩阵y计算最大可能λ
max
，并令最小可能λ
min
＝0计算惩罚参数
[0030]
2)、在给定预处理信号矩阵y及惩罚参数λ的情况下，求解如公式(3)所示的最小化优化问题，即求解信号矩阵x和断点个数p(x)；
[0031]
3)若断点数k小于断点个数p(x)时，令λ
min
＝λ并重复以上步骤；
[0032]
若断点数k大于断点个数p(x)时，令λ
max
＝λ并重复以上步骤；
[0033]
若断点数k等于断点个数p(x)时，输出信号矩阵x；
[0034]
当选择惩罚参数λ作为输入参数时：求解如公式(3)所示的最小化优化问题，即求解信号矩阵x。
[0035]
本发明还公开了一种面向多信号的快速断点检测方法的系统，包括：
[0036]
信号预处理模块，用于对获取的原始信号进行预处理；
[0037]
惩罚参数估计模块，用于在给定断点数或惩罚参数的情况下获取处理后的信号矩阵；
[0038]
信号处理模块，用于确定信号矩阵的断点位置和处理后的信号，实现分段信号的断点快速检测。
[0039]
一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行计算机程序时实现面向多信号的快速断点检测方法的步骤。
[0040]
一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现面向多信号的快速断点检测方法的步骤。
[0041]
与现有技术相比，本发明具有以下有益效果：
[0042]
本发明提出了一种面向多信号的处理方法，通过对采集的信号进行预处理操作，根据实际情况确定断点数和惩罚参数，就可以结合预处理结果和断点数及惩罚参数确定最终的断点位置和处理后的信号。该方法能够快速、准确地检测出信号中相同位置共享的断点，从而为信号进一步的分割、拟合、参数估计提供可靠的起始及终止位置信息。
[0043]
进一步地，采用二分法查找惩罚参数，首先与序列中间的元素进行比较，如果大于这个元素，就在当前序列的后半部分继续查找；如果小于这个元素，就在当前序列的前半部分继续查找，直到找到相同的元素或者所查找的序列范围为空为止，采用二分法优点是比较次数少，查找速度快，平均性能好。
[0044]
进一步地，采用归一化法对原始信号进行预处理可以提升收敛速度，同时为了保证预处理的信号的准确性，因此对原始信号进行标准化处理。
[0045]
本技术提出的一种面向多信号的快速断点检测方法的系统，将断点检测根据内容的关联性分解为不同且相互独立的模块进行实现，采用模块化思想实现多信号断点检测，当哪个模块出现问题时能单独管理，模块之间相互独立且互不影响。
附图说明
[0046]
图1为包含两个高低台阶的信号断点检测示意图；
[0047]
图2为多个信号中共有断点检测示意图；
[0048]
图3为本发明的快速断点检测方法流程图；
[0049]
图4为本发明应用于家庭共有拷贝数变异检测((a)子代继承了父母双亲共有的变
异；(b)子代继承了来自父亲的变异)；
[0050]
图5为本发明应用于可穿戴式设备的步态分析((a)使用x、y及z三个方向加速度传感器检测的步态；(b)只使用z方向加速度传感器检测的步态)；
[0051]
图6为本发明相比于传统方法对计算时间的改善((a)计算时间随信号长度n的改善；(b)计算时间随信号维数m的改善)。
具体实施方式
[0052]
为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。
[0053]
需要说明的是，本发明的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
[0054]
下面结合附图对本发明做进一步详细描述：
[0055]
本发明的提出的一种面向多信号的快速断点检测方法，具体信号处理流程如图3所示，快速断点检测方法包含以下步骤：
[0056]
s1、对原始信号进行归一化预处理，得到预处理后的n
×
m尺寸的预处理信号矩阵y；
[0057]
s1.1、将采集得到的原始信号存储入n
×
m的原始矩阵y0；
[0058]
s1.2、计算原始矩阵y0的最大绝对值c；
[0059]
s1.3、使用最大绝对值c对原始信号进行预处理操作。
[0060]
s2、根据实际情况确定断点数k或惩罚参数λ；
[0061]
s3、若选择输入参数为λ，则直接求解最小化优化问题，得到处理后的n
×
m尺寸的信号矩阵x；若选择输入参数为k，则首先根据预处理信号矩阵y计算最大可能λ
max
，之后在区间[0,λ
max
]内搜索估计惩罚参数λ，并求解给定λ下的最小化优化问题，使得断点数为k，得到处理后的n
×
m尺寸的信号矩阵x；
[0062]
以二分法搜寻区间[0,λ
max
]内的惩罚参数λ，使得检测断点数p(x)为k；
[0063]
s4、根据求得的断点位置计算信号矩阵x；对信号矩阵x乘以最大绝对值c进行去归一化后处理，得到处理后信号x0，从而实现分段信号的断点快速检测。
[0064]
其中，预处理结果如公式(1)所示：
[0065]
y＝y0/c
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0066]
处理后的信号x0的计算公式如公式(2)所示：
[0067]
x0＝cx
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0068]
在给定预处理信号矩阵y及惩罚参数λ的情况下，求解如公式(3)所示的最小化优化问题，即求解信号矩阵x：
[0069][0070]
其中，y是n
×
m的待处理信号矩阵，x是同尺寸的处理后信号矩阵，n是每个信号的采样点数(信号长度)，m是信号的个数，λ是对每个断点的惩罚参数，p(x)是x中的断点个数。当x求得，则断点位置及个数易知。
[0071]
当选择断点数k作为输入参数时：
[0072]
1)、根据预处理信号矩阵y计算最大可能λ
max
，并令最小可能λ
min
＝0计算惩罚参数
[0073]
2)、在给定预处理信号矩阵y及惩罚参数λ的情况下，求解如公式(3)所示的最小化优化问题，即求解信号矩阵x和断点个数p(x)；
[0074]
3)若断点数k小于断点个数p(x)时，令λ
min
＝λ并重复以上步骤；
[0075]
若断点数k大于断点个数p(x)时，令λ
max
＝λ并重复以上步骤；
[0076]
若断点数k等于断点个数p(x)时，输出信号矩阵x；
[0077]
当选择惩罚参数λ作为输入参数时：求解如公式(3)所示的最小化优化问题，即求解信号矩阵x。
[0078]
上述技术方案的核心步骤是s3，s3的具体操作步骤如下：
[0079]
a)输入：预处理信号矩阵y及惩罚参数λ；
[0080]
b)初始化：目标函数存储向量f为n 1长的全零向量，断点存储数组bp为元胞数组，其第一个元胞为空向量，有效索引列表r＝1，分段能量e＝0，平均值z等于预处理信号矩阵y的第一列；
[0081]
c)进入i＝1的循环，直到步骤x)；
[0082]
d)将目标函数存储向量f的有效索引列表的第r个元素与分段能量e相加，存入临时向量v；
[0083]
e)搜寻临时向量v的最小值及位置，并将临时向量v的最小值及位置分别存入a及i1；
[0084]
f)计算a λ，存入目标函数存储向量f的第i 1个存储单元中；
[0085]
g)读取有效索引列表r的第i1个元素，存入最大绝对值c中；
[0086]
h)读取断点存储数组bp的第c个元胞，并与最大绝对值c拼接，存入断点存储数组bp的第i 1个元胞；
[0087]
i)寻找临时向量v中小于目标函数存储向量f中第i 1个元素的位置，将其存入i2；
[0088]
j)如果i小于采样点数n，则逐条执行到步骤x)，否则直接跳到步骤x)；
[0089]
k)将预处理信号矩阵y的第i 1列存入y中；
[0090]
l)在有效索引列表r中只保留第i2个元素，删除 其余元素；
[0091]
m)计算i减去有效索引列表r并加1，并存入长度向量l；
[0092]
n)将长度向量l以列方式复制m次，存入矩阵l中；
[0093]
o)将y以行方式复制n次(n为长度向量l的长度)，存入矩阵b中；
[0094]
p)读取平均值z的i2行，存入矩阵t中；
[0095]
q)计算矩阵b减矩阵t，对每一行计算所有元素的平方和，存入向量e中；
[0096]
r)计算向量e点乘长度向量l，再点除l 1，存入向量e中；
[0097]
s)计算向量e与分段能量e中第i2位置元素之和，存入分段能量e中；
[0098]
t)在分段能量e末尾增添0；
[0099]
u)计算矩阵t点乘矩阵l，再加矩阵b，再点除(l 1)，存入平均值z中；
[0100]
v)在平均值z末尾列增添y；
[0101]
w)在有效索引列表r末尾增添i 1；
[0102]
x)i加1，返回步骤c)；
[0103]
y)输出：输出断点位置为断点存储数组bp最后一个元胞中第2到最后一个元素；输出处理后数据矩阵x以列为单位进行处理，对每列信号，用相邻两个断点间数据矩阵y中该列的平均值作为x该列该段的信号值。
[0104]
本发明还公开了一种面向多信号的快速断点检测方法的系统，包括：
[0105]
信号预处理模块，用于对获取的原始信号进行预处理；
[0106]
惩罚参数估计模块，用于在给定断点数或惩罚参数的情况下获取处理后的信号矩阵；
[0107]
信号处理模块，用于确定信号矩阵的断点位置和处理后的信号，实现分段信号的断点快速检测。
[0108]
一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行计算机程序时实现面向多信号的快速断点检测方法的步骤。
[0109]
一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现面向多信号的快速断点检测方法的步骤。
[0110]
如图4所示，展示了本发明应用于高通量测序技术的群体共有拷贝数变异检测有益效果。该测序数据源于一个三口之家(父、母、子，m＝3)，数据的前处理包括(i)使用映射工具bowtie将hts下机数据文件fastq对比对到参考基因组hg19；(ii)计算读深(read depth，rd)信号，其为基因组上固定宽度窗口内每个碱基位点的测序覆盖深度；(iii)对rd信号的gc及可映射性矫正。图4(a)展示了使用本发明在22号染色体32.8～33.4mb区间内检测到的拷贝数变异。可以看到，子代继承了父母双亲共有的变异。更进一步，图4(b)展示了使用本发明在22号染色体39.2～39.5mb区间内检测到的拷贝数变异。可以看到，这里子代继承了来自父亲的变异，母亲在该处无变异。说明使用多个信号检测断点具有重要的应用价值。
[0111]
如图5所示，展示了本发明应用于可穿戴式设备的健康数据分析。该数据源于一名被试12秒的跑步运动，传感器安装于被试的右踝关节，检测到x、y及z三个方向(m＝3))的加速度。图5(a)展示了使用本发明对三个加速度信号的检测效果，可以看到，步与步之间被很好检测。作为对比，图5(b)展示了只使用一个加速度信号(z方向)的检测效果，可以看到，在8秒左右的步态不能被很好地区分。说明使用多个信号可以提升断点的检测精度。
[0112]
如图6所示，展示了本发明相比于传统方法对计算时间的改善。这里使用了仿真数据。图6(a)展示了信号个数m＝10时，计算时间与信号长度n的关系。可以看到，随着信号长度n的增加，计算时间随之增加，且本发明只使用传统方法约百分之一的计算时间。对于n＝
100000点数的极长信号，本发明只使用约10秒，而传统方法则需使用约1000秒。图6(b)展示了信号长度n＝3000时，计算时间与信号维数m的关系。可以看到，随着信号维数m的增加，本发明所需的计算时间几乎保持不变，而传统方法的计算时间随之呈线性增加。对于m＝1000的众多信号，本发明计算时间少于1秒，而传统方法则需使用约400秒。说明本发明可大幅降低计算时间，尤其适用于大量的信号。
[0113]
本发明提出了一种快速信号处理方法。该方法可以快速、准确地检测出多维信号中共有的断点位置，进而对多维信号进行分割、拟合、参数估计提供可靠的起始及终止位置信息。该发明技术具有广泛的实用特点，该方法可以应用于生物、医学、工程等领域，如基于高通量测序技术的群体拷贝数变异检测、基于可穿戴设备的运动状态检测等。
[0114]
本发明出现和使用的缩略语和关键术语定义如下：
[0115]
cnv
ꢀꢀꢀꢀ
copy number variation
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
拷贝数变异
[0116]
hts
ꢀꢀꢀꢀ
high
‑
throughput sequencing
ꢀꢀꢀꢀꢀ
高通量测序
[0117]
dp
ꢀꢀꢀꢀꢀ
dynamic programming
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
动态规划
[0118]
rd
ꢀꢀꢀꢀꢀ
read depth
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
读深
[0119]
cbs
ꢀꢀꢀꢀ
circular binary segmentation
ꢀꢀꢀ
循环二分分割
[0120]
op
ꢀꢀꢀꢀꢀ
optimal partitioning
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
最优分割
[0121]
pelt
ꢀꢀꢀ
pruned exact linear time
ꢀꢀꢀꢀꢀꢀꢀ
删减的精确线性时间
[0122]
tpr
ꢀꢀꢀꢀ
true positive rate
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
真阳性率
[0123]
fpr
ꢀꢀꢀꢀ
false positive rate
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
假阳性率
[0124]
gc
ꢀꢀꢀꢀꢀ
guanine
‑
cytosine content
ꢀꢀꢀꢀꢀꢀꢀ
鸟嘌呤
‑
胞嘧啶含量
[0125]
以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。