一种基于非线性拟合方式的肽质谱峰特征参数提取方法

摘要

本发明涉及一种肽质谱峰特征参数提取方法。现有方法针对在肽段质谱图中形成谱峰的各样点其分布存在较大偏差时，存在难以保证所提取出的质谱峰特征参数精准度的不足。本发明提出基于非线性拟合方式的肽质谱峰特征参数提取方法，利用多个样点数据，以实际数据与拟合结果之间差值最小为导向，采用迭代方法不断更新特征参数估计值，直至满足收敛条件，从而获得最终特征参数估值。该方法可有效减少样点分布偏差对高斯曲线特征参数求解带来的不利影响，提高特征参数数值精准度，进而有利于肽段鉴定精度的改善。

主权项

一种基于非线性拟合方式的肽质谱峰特征参数提取方法，其特征在于：设质谱图中某离子的高斯峰由N个样点组成，N≥3；对样点按其丰度值从大到小排序后，其坐标构成集合A；A＝{(m₁,d₁),(m₂,d₂),…(m_N,d_N)}其中，m_i表示质荷比，d_i表示丰度值，i∈{1,2,3,…,N}；准备通过样点拟合出的高斯曲线其函数形式设为：$f(x,P)=p_1\times e^{-{\left(\frac{x-p_2}{p_3}\right)}^2}$其中，函数f(x，P)代表理论丰度值，自变量x代表质荷比，p₁、p₂和p₃为待求解的高斯曲线特征参数，分别表征缩放因子、质心、标准差，构成特征参数向量P＝[p₁ p₂ p₃]；具体步骤如下：步骤(1)根据丰度值最大的3个样点数据，对高斯曲线特征参数赋初值；$p_1=\frac{1}{3}\times [d_1\times e^{{\left(\frac{m_1-p_2}{p_3}\right)}^2}+d_2\times e^{{\left(\frac{m_2-p_2}{p_3}\right)}^2}+d_3\times e^{{\left(\frac{m_3-p_2}{p_3}\right)}^2}]$$p_2=\frac{1}{2}\times \frac{[\ln \left(d_2\right)-\ln \left(d_3\right)]\times m_1^2+[\ln \left(d_3\right)-\ln \left(d_1\right)]\times m_2^2+[\ln \left(d_1\right)-\ln \left(d_2\right)]\times m_3^2}{[\ln \left(d_2\right)-\ln \left(d_3\right)]\times m_1+[\ln \left(d_3\right)-\ln \left(d_1\right)]\times m_2+[\ln \left(d_1\right)-\ln \left(d_2\right)]\times m_3}$$p_3=\frac{1}{2}\times [\sqrt{\frac{{(m_2-p_2)}^2-{(m_1-p_2)}^2}{\ln \left(d_1\right)-\ln \left(d_2\right)}}+\sqrt{\frac{{(m_3-p_2)}^2-{(m_2-p_2)}^2}{\ln \left(d_2\right)-\ln \left(d_3\right)}}]$其中，ln(·)表示取自然对数操作；步骤(2)选择合适数值初始化迭代步长参数λ，该参数初始化数值的大小将影响迭代次数和收敛速度；步骤(3)计算拟合结果误差Err，判定迭代过程是否结束；$\mathrm{Err}=\underset{i=1}{\overset{N}{\Sigma }}{[f(m_i,P)-d_i]}^2$设定判决门限ε₁，如果Err≤ε₁，则处理过程结束，当前向量P中的特征参数值即为求解的最终结果；反之，如果Err＞ε₁，则进入步骤(4)；步骤(4)根据当前特征参数向量P，构造矩阵J；$J=\left[\begin{array}{ccc}\frac{\partial f(m_1,P)}{{\partial p}_1}& \frac{\partial f(m_1,P)}{{\partial p}_2}& \frac{\partial f(m_1,P)}{{\partial p}_3}\\ \frac{\partial f(m_2,P)}{{\partial p}_1}& \frac{\partial f(m_2,P)}{{\partial p}_2}& \frac{\partial f(m_2,P)}{{\partial p}_3}\\ .& .& .\\ .& .& .\\ .& .& .\\ \frac{\partial f(m_N,P)}{{\partial p}_1}& \frac{\partial f(m_N,P)}{{\partial p}_2}& \frac{\partial f(m_N,P)}{{\partial p}_3}\end{array}\right]$步骤(5)计算每次迭代过程中，特征参数向量P的更新向量H＝[Δp₁ Δp₂ Δp₃]^T，Δp₁，Δp₂和Δp₃分别为特征参数p₁，p₂和p₃的待定更新值；构造误差矢量E；E＝[d₁‑f(m₁,P),d₂‑f(m₂,P),…d_N‑f(m_N,P)]^T则：H＝[J^T×J+λ×diag(J^T×J)]^‑1×J^T×E其中，diag(·)表示矩阵对角元素提取和创建对角阵操作；步骤(6)计算更新向量H的度量值ρ(H)；$\rho \left(H\right)=\frac{\underset{i=1}{\overset{N}{\Sigma }}{[f(m_i,P)-d_i]}^2-\underset{i=1}{\overset{N}{\Sigma }}{[f(m_i,P+H)-d_i]}^2}{{2H}^T\times (\lambda \times H+J^T\times E)}$步骤(7)更新特征参数向量P和迭代步长参数λ；设定判决门限ε₂，如果更新向量H的度量值ρ(H)＞ε₂，则当前特征参数向量P数值由P+H替代，即P←P+H，完成更新，同时当前迭代步长参数λ数值减小至λ/K，即λ←λ/K；反之，如果ρ(H)≤ε₂，则当前特征参数向量P保持不变，同时迭代步长参数λ数值增加K倍，即λ←K×λ；K为比例因子，取值范围为5～20；在完成特征参数向量P和迭代步长参数λ更新后，返回至步骤(3)，进行下一轮迭代。