CN114487072A - 一种飞行时间质谱峰拟合方法 - Google Patents

Abstract

本发明公开了一种飞行时间质谱峰拟合方法，包括：将整条飞行时间核酸质谱图按照一定间隔长度进行分段，并计算出各分段强度阈值和二阶导数阈值：在每个分段内找到离子强度过零点的一阶导数，确定可能的高斯峰位置；从可能的高斯峰位置中，选取符合条件的离子强度，得到确定的高斯峰位置，确定出所有的高斯峰个数l，及高斯峰的质荷比m/z_i和峰高H_i；根据半经验公式，计算得到高斯峰的峰宽W_i；将高斯峰的质荷比m/z_i和峰宽W_i，取4W_i作为高斯拟合区域，作为高斯峰最小二乘法迭代初始值；将高斯拟合区域4W_i内的所有采样点的核酸质谱数据带入高斯拟合函数f(x_i)中进行迭代求解，得到精确的拟合高斯峰的质荷比m/z_fit、峰宽W_fit、峰高H_fit。

Description

一种飞行时间质谱峰拟合方法

技术领域

本发明属于质谱分析数据处理技术领域，具体涉及一种飞行时间质谱峰拟合方法。

背景技术

飞行时间质谱是一种高通量、高精度、检测质量范围宽等特性的基因检测技术。样本在激光等外在条件的激发下离子化，通过加速电场加速，然后自由飞行一定距离到达离子检测器；经过相关转换生成质谱数据，随后再对质谱数据信息进行处理完成基因位点诊断。其中对质谱数据的处理主要是高斯峰的确定拟合及相关参数的提取。

质谱峰拟合一般分为两步：第一步质谱峰的检测及重叠峰数量的判断；第二步检测到的质谱峰的拟合及特性参数提取。目前，第一步质谱峰的检测及重叠峰数量的判断方法主要包括傅里叶变化的去卷积技术和基于傅里叶变换技术扩展的小波变化技术，但是它们在处理过程中都会增加噪音降低信噪比同时小波技术还严重依赖于小波函数的选取、小波分解层数确定等问题。第二步检测到的质谱峰的拟合及特性参数提取方法主要包括最小二乘法高斯拟合，但是目前基于最小二乘法的高斯拟合结果严重依赖于原始数据点的，在原始数据点存在分布偏差较大噪音较明显的情况下很难获取准确度拟合结果。

发明内容

发明目的：为解决现有质谱峰拟合存在拟合结果准确度低的问题，本发明提出了一种飞行时间质谱峰拟合方法。

技术方案：一种飞行时间质谱峰拟合方法，包括以下步骤：

步骤1：获取整条飞行时间核酸质谱图，该飞行时间核酸质谱图，以采样点的质荷比为横坐标，以采样点的离子强度为纵坐标；对于第i个采样点，i∈{1,2,3,…，n}，其核酸质谱数据表示为(I_i,m/z_i)，其中，I_i表示离子强度，m/z_i表示质荷比；

步骤2：将整条飞行时间核酸质谱图按照一定间隔长度进行分段，并计算出各分段强度阈值和二阶导数阈值，对于第j分段，其强度阈值表示为I_{threshold j}和二阶导数阈值表示为I″_{threshold j}：

步骤3：在每个分段内找到离子强度过零点的一阶导数，该一阶导数对应的原始的质荷比，为可能的高斯峰位置；

步骤4：从可能的高斯峰位置中，选取符合以下条件的离子强度，该离子强度对应的原始的质荷比，为确定的高斯峰位置：

且

式中，I′_(i+1)j为第j分段内第i+1采样点的离子强度I(_i+1)_j的一阶导数I′_(i+1)j；

确定出所有的高斯峰个数l，及高斯峰的质荷比

和峰高H_i，峰高H_i由以下公式计算得到：

对于第i个采样点，当质谱基线c_i为零，峰高H_i＝I_i；当质谱基线c_i不为零，峰高H_i＝I_i-c_i；

步骤5：根据半经验公式，计算得到高斯峰的峰宽W_i＝2×0.001×m/z_i；

步骤6：将高斯峰的质荷比m/z_i和峰宽W_i，取4W_i作为高斯拟合区域，作为高斯峰最小二乘法迭代初始值；

步骤7：将高斯拟合区域4W_i内的所有采样点的核酸质谱数据带入高斯拟合函数f(x_i)中进行迭代求解：

其中，H_fit、m/z_fit、W_fit分别为拟合高斯峰的峰高、质荷比和峰宽，α为特征峰宽系数；

质谱基线c_fit根据下式进行求解：

c_fit＝A\I (9)

其中，I为由高斯拟合区域内所有采样点的离子强度I_i组成的向量，表示为I＝(I₁,I₂,I₃,…,I_N)，N为高斯拟合区域内所有采样点的总数，A为向量内部所有数值均为1且长度与向量I相同，表示为A＝(1₁,1₂,1₃,…,1_N)；

将求解得到的拟合高斯峰的峰高H_fit和质荷比m/z_fit与作为迭代初始值的质荷比m/z_i和峰高H_i进行匹配，判断以下条件是否成立：

|m/z_i-m/z_fit|≤Δ₁和|H_i-H_fit|≤Δ₂或者

式中，

为第k次和第k+1次迭代拟合求解得到的质荷比，k为大于1的正整数，Δ₁、Δ₂和Δ₃为微小量常数；

若不成立，则将拟合高斯峰的峰高H_fit和质荷比m/z_fit作为下一次迭代初始值，进行新一轮的求解直到满足式(10)的条件；

若成立，则拟合求解结束，得到精确的拟合高斯峰的质荷比m/z_fit、峰宽W_fit、峰高H_fit。

进一步的，步骤3中，所述的在每个分段内找到离子强度过零点的一阶导数，具体步骤为：

对于第j分段内第个采样点i，其一阶导数表示为I′_ij，数学表达式为：

式中，I_i+1j表示第j分段内第i+1个采样点的离子强度，I_i-1j表示第j分段内第i-1个采样点的离子强度；

对于第1个采样点，其一阶导数表示为：

I′₁＝I₂-I₁ (4)

式中，I₂表示第2个采样点的离子强度，I₁表示第1个采样点的离子强度；

对于第n个采样点，其一阶导数表示为：

I′_n＝I_n-I_n-1 (5)

式中，I_n表示第n个采样点的离子强度，I_n-1表示第n-1个采样点的离子强度；

对每个采样点进行一阶导数计算；

基于一阶导数计算，在每个分段内找到离子强度过零点的一阶导数。

进一步的，步骤2中，对于第j分段，按照下式计算强度阈值I_{threshold j}和二阶导数阈值I″_{threshold j}：

式中，

为第j分段内离子强度的平均值，σ_j为第j分段内离子强度的标准差；a为信号强度阈值系数；d为特定系数，H_j和W_j为第j分段内峰高最小的特征峰的峰高和半高半宽。

有益效果：本发明方法，克服了在寻峰及重叠峰判断时存在噪音增大，影响信噪比的问题，以及在根据原始质谱数据进行最小二乘法拟合高斯峰情况下，本发明引入在寻找一阶导数过零点的质谱点确定质谱峰位置，准确记录下原始数据中质谱峰高和位置作为高斯拟合结果的判定条件，确保拟合结果的准确性。

附图说明

图1为本发明方法的流程图；

图2为采用本发明方法估算出峰的位置质谱图；

图3为采用本发明方法得到的高斯拟合结果质谱图。

具体实施方式

现结合附图和实施例对本发明的技术方案做进一步说明。

如图1所示，本发明的一种飞行时间质谱峰拟合方法，主要包括以下步骤：

步骤1：获取整条飞行时间核酸质谱图信号，对于第i个采样点，i∈{1,2,3,…，n}，其核酸质谱数据表示为(I_i,m/z_i)，其中，I_i表示离子强度，也称为离子丰度，m/z_i表示离子质荷比，与离子强度I_i一一对应；

步骤2：将整条飞行时间核酸质谱图信号按照一定间隔长度进行分段，并根据下式，计算得到各分段强度阈值和二阶导数阈值，对于第j分段，其强度阈值表示为I_{threshold j}和二阶导数阈值表示为I″_{threshold j}：

式中，

为第j分段内离子强度的平均值，σ_j为第j分段内离子强度的标准差；a为信号强度阈值系数，a∈{1,2,3}，优选为2；d为特定系数，取值为2，H_j和W_j为第j分段内峰高最小的特征峰的峰高和半高半宽，当很难评估出H_j和W_j可取值为0；

步骤3：在每个分段内找到信号强度过零点的一阶导数，其对应的原始离子质荷比，即为可能的高斯峰位置；在本步骤中，采用以下公式求解各采样点的一阶导数：

对于第j分段内第个采样点i，其一阶导数表示为I′_ij，数学表达式为：

式中，I_i+1j表示第j分段内第i+1个采样点的离子强度，I_i-1j表示第j分段内第i-1个采样点的离子强度；

对于第1个采样点，其一阶导数表示为：

I′₁＝I₂-I₁ (4)

式中，I₂表示第2个采样点的离子强度，I₁表示第1个采样点的离子强度；

对于第n个采样点，其一阶导数表示为：

I′_n＝I_n-I_n-1 (5)

式中，I_n表示第n个采样点的离子强度，I_n-1表示第n-1个采样点的离子强度。

步骤4：从可能的高斯峰位置中，选取符合以下条件的质谱峰(离子强度)，其对应的原始离子质荷比，即为确定的高斯峰位置：

式中，I′_(i+1)j为第j分段内第i+1采样点的离子强度I_(i+1)j的一阶导数I′_(i+1)j。

基于确定的高斯峰位置，确定出所有可能的高斯峰个数l，位置m/z_i和峰高H_i，峰高H_i由以下公式计算得到：

对于第i个采样点，当质谱基线c_i为零，峰高H_i＝I_i；当质谱基线c_i不为零，峰高H_i＝I_i-c_i，该质谱基线c_i为质谱峰(离子强度)I_i对应的基线。

步骤5：根据半经验公式，高斯峰半峰高的半宽是其峰位置(离子质荷比)m/z_i的0.1％，也就是说，峰宽W_i＝2×0.001×m/z_i；在高斯拟合过程中一般设定高斯拟合区域为4倍的峰宽W_i，当两个峰位置(离子质荷比)m/z_i之差小于4W_i则认为是重叠峰，否则认为是单峰；初步判定是否存在重叠峰及重叠峰中谱峰的数量d；

步骤6：将估算出的高斯峰位置m/z_i和峰宽W_i，作为高斯峰最小二乘法迭代初始值；

利用单纯性、算法Levenberg-Marquardt算法或者信赖与反射算法将高斯拟合区域4W_i内的(I_i,m/z_i)数据点带入高斯拟合函数f(x_i)中进行迭代求解；

高斯拟合函数f(x_i)的表达式为：

其中，H_fit、m/z_fit、W_fit分别为拟合高斯峰高度、位置和宽度，α为特征峰宽系数，取值为0.60056120439323。

最小二乘法迭代最优收敛公式为：

λ²＝∑(f(x_i)-I_i)² (8)

基线强度c_fit根据Flat Baseline Correction方法进行求解：

c_fit＝A\I (9)

其中，I为由高斯拟合区域内所有采样点的离子强度I_i组成的向量，表示为I＝(I₁,I₂,I₃,…,I_N)，N为高斯拟合区域内所有采样点的总数，A为向量内部所有数值均为1且长度与向量I相同，表示为A＝(1₁,1₂,1₃,…,1_N)。

将获取高斯峰特征参数(拟合质荷比m/z_fit和拟合峰高H_fit)与估算出的质荷比m/z_i和峰高H_i进行匹配，确认以下条件是否成立：

|m/z_i-m/z_fit|≤Δ₁和|H_i-H_fit|≤Δ₂或者

式中，

为第k次和第k+1次迭代拟合求解的质荷比，k为大于1的正整数，Δ₁、Δ₂和Δ₃为微小量常数。

当不成立，将拟合结果m/z_fit和W_fit作为下一次迭代初始值，进行新一轮的求解直到满足式(10)的条件，则拟合求解结束，从而获得精确的高斯峰特征参数：m/z_fit、W_fit、H_fit。

现通过一具体实施例说明本发明方法。

本实施例中采用真实飞行时间图谱数据如图2所示，横坐标为质荷比，纵坐标为离子强度，图2中存在由两个高斯峰叠加形成的重叠峰。

峰1和峰2是由本发明方法寻找到高斯峰的位置和个数，位置和峰高具体参数如下表1；

表1原始估算质谱峰参数

	峰位置	峰高
			峰1	9239.4	17
峰2	9256.3	17.4

根据半经验公式，峰1的峰宽为18.48，峰2的峰宽为18.51。

峰1和峰2位置之差为16.9小于4倍的峰1宽度，所以峰1和峰2为重叠峰；

对峰1和峰2进行高斯拟合，拟合结果如表2所示：

表2拟合质谱峰参数

	峰位置	峰高	峰宽
				峰1	9240.5	16.028	12.372
峰2	9256.1	16.366	12.256

Claims (3)

1.一种飞行时间质谱峰拟合方法，其特征在于：包括以下步骤：

步骤1：获取整条飞行时间核酸质谱图，该飞行时间核酸质谱图，以采样点的质荷比为横坐标，以采样点的离子强度为纵坐标；对于第i个采样点，i∈{1,2,3,…，n}，其核酸质谱数据表示为(I_i,m/z_i)，其中，I_i表示离子强度，m/z_i表示质荷比；

步骤2：将整条飞行时间核酸质谱图按照一定间隔长度进行分段，并计算出各分段强度阈值和二阶导数阈值，对于第j分段，其强度阈值表示为I_thresholdj和二阶导数阈值表示为I″_thresholdj：

步骤3：在每个分段内找到离子强度过零点的一阶导数，该一阶导数对应的原始的质荷比，为可能的高斯峰位置；

步骤4：从可能的高斯峰位置中，选取符合以下条件的离子强度，该离子强度对应的原始的质荷比，为确定的高斯峰位置：

式中，I′_(i+1)j为第j分段内第i+1采样点的离子强度I_(i+1)j的一阶导数I′_(i+1)j；

确定出所有的高斯峰个数l，及高斯峰的质荷比m/z_i和峰高H_i，峰高H_i由以下公式计算得到：

对于第i个采样点，当质谱基线c_i为零，峰高H_i＝I_i；当质谱基线c_i不为零，峰高H_i＝I_i-c_i；

步骤5：根据半经验公式，计算得到高斯峰的峰宽W_i＝2×0.001×m/z_i；

步骤6：将高斯峰的质荷比m/z_i和峰宽W_i，取4W_i作为高斯拟合区域，作为高斯峰最小二乘法迭代初始值；

步骤7：将高斯拟合区域4W_i内的所有采样点的核酸质谱数据带入高斯拟合函数f(x_i)中进行迭代求解：

其中，H_fit、m/z_fit、W_fit分别为拟合高斯峰的峰高、质荷比和峰宽，α为特征峰宽系数；

质谱基线c_fit根据下式进行求解：

c_fit＝A\I (9)

其中，I为由高斯拟合区域内所有采样点的离子强度I_i组成的向量，表示为I＝(I₁,I₂,I₃,…,I_N)，N为高斯拟合区域内所有采样点的总数，A为向量内部所有数值均为1且长度与向量I相同，表示为A＝(1₁,1₂,1₃,…,1_N)；

将求解得到的拟合高斯峰的峰高H_fit和质荷比m/z_fit与作为迭代初始值的质荷比m/z_i和峰高H_i进行匹配，判断以下条件是否成立：

式中，

为第k次和第k+1次迭代拟合求解得到的质荷比，k为大于1的正整数，Δ₁、Δ₂和Δ₃为微小量常数；

若不成立，则将拟合高斯峰的峰高H_fit和质荷比m/z_fit作为下一次迭代初始值，进行新一轮的求解直到满足式(10)的条件；

若成立，则拟合求解结束，得到精确的拟合高斯峰的质荷比m/z_fit、峰宽W_fit、峰高H_fit。

2.根据权利要求1所述的一种飞行时间质谱峰拟合方法，其特征在于：步骤3中，所述的在每个分段内找到离子强度过零点的一阶导数，具体步骤为：

对于第j分段内第个采样点i，其一阶导数表示为I′_ij，数学表达式为：

式中，I_i+1j表示第j分段内第i+1个采样点的离子强度，I_i-1j表示第j分段内第i-1个采样点的离子强度；

对于第1个采样点，其一阶导数表示为：

I′₁＝I₂-I₁ (4)

式中，I₂表示第2个采样点的离子强度，I₁表示第1个采样点的离子强度；

对于第n个采样点，其一阶导数表示为：

I′_n＝I_n-I_n-1 (5)

式中，I_n表示第n个采样点的离子强度，I_n-1表示第n-1个采样点的离子强度；

对每个采样点进行一阶导数计算；

基于一阶导数计算，在每个分段内找到离子强度过零点的一阶导数。

3.根据权利要求1所述的一种飞行时间质谱峰拟合方法，其特征在于：步骤2中，对于第j分段，按照下式计算强度阈值I_thresholdj和二阶导数阈值I″_thresholdj：

式中，

为第j分段内离子强度的平均值，σ_j为第j分段内离子强度的标准差；a为信号强度阈值系数；d为特定系数，H_j和W_j为第j分段内峰高最小的特征峰的峰高和半高半宽。

Images