一种基于非线性拟合方式的肽质谱峰特征参数提取方法
【专利摘要】本发明涉及一种肽质谱峰特征参数提取方法。现有方法针对在肽段质谱图中形成谱峰的各样点其分布存在较大偏差时，存在难以保证所提取出的质谱峰特征参数精准度的不足。本发明提出基于非线性拟合方式的肽质谱峰特征参数提取方法，利用多个样点数据，以实际数据与拟合结果之间差值最小为导向，采用迭代方法不断更新特征参数估计值，直至满足收敛条件，从而获得最终特征参数估值。该方法可有效减少样点分布偏差对高斯曲线特征参数求解带来的不利影响，提高特征参数数值精准度，进而有利于肽段鉴定精度的改善。
【专利说明】一种基于非线性拟合方式的肽质谱峰特征参数提取方法
【技术领域】
[0001]本发明属于生物质谱数据预处理及信息提取【技术领域】，具体涉及一种基于非线性拟合方式的肽质谱峰特征参数提取方法。
【背景技术】
[0002]基于串联质谱的肽鉴定是目前蛋白质组研究领域中广泛使用的技术。待鉴定的肽在质谱仪中被碎裂为碎片离子，从而生成串联质谱数据，并与理论串联质谱库或已鉴定的肽段质谱库进行比对及分析，最后完成对未知肽段的鉴定。
[0003]通常情况下对某种离子进行质谱检测，所检测到的质荷比数据不是单一数值点，而是存在若干样点，在质谱图上其拟合为高斯曲线，即高斯峰。为确定该离子的荷质比，需对这些样点进行预处理，计算出其横轴方向上的质心(Centroid)，即该离子的实测质荷比。根据所求质心，可进而推算出该离子最大丰度值等其他特征参数。目前质心求解方法有多种，比较常见的思路是:假定质谱图上构成高斯峰的各个样点均严格分布在某条高斯曲线上，利用各样点的数值(质荷比和丰度值)，代入到参数未知的通用高斯曲线函数表达式中，构造联立方程组，从而解出相应高斯峰的特征参数，包括质心，最大丰度值等。当前应用极为广泛的一款蛋白质组学数据分析软件MAXQUANT采用的即是这一方法。然而在实际检测中，受实验条件、所在环境以及仪器设备噪声等因素的影响，质谱图上各个样点往往并非严格分布在高斯曲线上，而是存在一定偏差。当各个样点偏差数值较大，则上述方法中的假设条件难以成立，因而势必造成求解出的特征参数在数值上存在较大误差，进而影响到肽段鉴定的精度。
【发明内容】
[0004]本发明的目的在于解决上述方法的缺点和不足，提出一种基于非线性拟合方式的肽质谱峰特征参数提取方法。
[0005]设质谱图中某离子的高斯峰由N个样点组成，通常情况下N彡3。对样点按其丰度值从大到小排序后，其坐标构成集合A。
[0006]A = {(Hi1, (I1)，(m2, d2)，…(mN, dN)}
[0007]其中，Hii表示质荷比，(Ii表示丰度,其值大于O, i e {1,2，3，...，N}。需要通过样点拟合出的高斯曲线其函数形式设为:
[0008]^ -(T)
/(.V, P) = p{ Xe
[0009]其中，函数f(x，p)代表丰度值，自变量X代表质荷比，Pl、PjP P3S待求解的高斯曲线特征参数，分别表征缩放因子、质心、标准差，构成特征参数向量P = [P1 P2 P3I °所述的特征参数提取方法处理步骤如下:
[0010]步骤(I)根据丰度值最大的3个样点数据，对高斯曲线特征参数赋初值。
(呷-]).|~(所:-匕 |~(衍3 ~^2 |~
[0011] ^^ ^ 十32 X 4 ^」+ ^3 ^6 ^
1 卜⑷)—)] X ！11\ + [111(4 )—卜1“,! )] 乂 1111 + [1^(^, ) — 1*1(4 )] X
〔 〕 巧 2 卜⑷)—!!!“/;)^"、+(10(4)-(门⑷)…"。+^1^(^/,)-^]^/^；
[0013]卜=1?|』卜—巧)2+1(^ — ^ )2-(^ — ^ )2
[0014]其中，IV 0表示取自然对数操作。
[0015]步骤(2)选择合适数值初始化迭代步长参数入，该参数初始化数值的大小将影响本处理方法的迭代次数和收敛速度。
[0016]步骤(3)计算拟合结果误差此!"，判定迭代过程是否结束。
[0017]￡^ = ^/(^,., ^)-^,-]2
卜1
[0018]设定判决门限、，如果此!"彡8工，则处理过程结束，当前向量？中的特征参数值即为求解的最终结果。反之，如果此1~? 8工，则进入步骤(幻。
[0019]门限8 ,的取值决定了所提取的特征参数值的精准度，同时影响处理过程的迭代次数。8工取值越小，特征参数值的精准度就越高，处理所需迭代次数越多。需要注意的是，如果8工取值过小，则该迭代过程可能将无法最终收敛。反之，8工取值越大，所提取的参数精准度将相应降低，而迭代次数将会减少。
[0020]步骤⑷根据当前特征参数向量？，构造矩阵了。
0” 卩)(:/.(坩丨,？?0丨，？?
8(), €/),
0/ (1112 , 11)0/ (1112 ,
[0021]1= 01^0^25/?
￠7/.(川、,卩)句'11.1、, 0) ￡/ (111.,，卩)
— 5/^,5厂2(? 一
[0022]步骤(5)计算每次迭代过程中，特征参数向量？的更新向量!I =[八巧八％八巧]7，八巧，八巧和八分别为特征参数巧，和的待定更新值。构造误差矢量2。
[0023]2 = (叫，户)，(叫，？；?，…0%，户)11
[0024]则:
[0025]11= [了?了+入父乜叫(了?了)]—1※】1※^
[0026]其中，(11^( ? )表示矩阵对角元素提取和创建对角阵操作。
[0027]步骤(6)计算更新向量！！的度量值9⑶。
「 ^ [[/(坩” ？卜式]2-[[,(岬’ ^ + 0)-^-]2〔0028〕 /7(?) =^~-~,-
[0029]步骤(7)更新特征参数向量P和迭代步长参数λ。设定判决门限ε2，如果更新向量H的度量值P⑶> ε2，则当前特征参数向量P数值由Ρ+Η替代，即P —Ρ+Η，完成更新，同时当前迭代步长参数λ数值减小至λ/K，S卩λ — λ/Κ。反之，如果P (H) ( ε2，则当前特征参数向量P保持不变，同时迭代步长参数λ数值增加K倍，即λ — ΚΧλ。K为比例因子，一般取值范围为5?20。判决门限ε2应根据样点偏差程度、收敛速度要求等具体指标设置合适数值。
[0030]在完成特征参数向量P和迭代步长参数λ更新后，返回至步骤(3)，进入行下一轮迭代。
[0031]本发明中肽质谱峰特征参数提取方法，采用多样点非线性拟合方式求解特征参数，减少了样点分布偏差带来的不利影响，提升了参数提取精准度，进而有利于肽段鉴定精度的改善。
【具体实施方式】
[0032]步骤(I)根据丰度值最大的3个样点数据，对高斯曲线特征参数赋初值。
II rn「p2 j~I m1-Pi Γ? ?H-Pi
[0033]P1 = Jx xe" ρ> J +d2xe- Pl ^ +d^xe" Pi」
I [?η(?/9) - ln(c/,)] X mf + [Ιη(?/,) - \χ\{?{)] χ ηι?, + [lr^i/,) - Ιη(?/,,)] χ m\
[ ]2 [?η(?/7)- Ιη(?/,)] χ + [ln(t/,) - ln(f/,)] χ nu + [?η(?/,) -\n(cL)] χ m,
I Ι( nu — ) " - (ηι — ) " I(/η, - ρΛ - ? nu - Γ
[0035]P3 =-x Jv 二 ' J~J +J J~h
3 2 ]j ln(i/,)-ln(r/:) ]j In(V2)-1n(Ji)
[0036]其中，ln(.)表示取自然对数操作。
[0037]步骤(2)选择合适数值初始化迭代步长参数λ，该参数初始化数值的大小将影响本处理方法的迭代次数和收敛速度。
[0038]步骤(3)计算拟合结果误差Err，判定迭代过程是否结束。
N?
[0039]Err = Yj [/ (m,, P) - cl ]"
/-1
[0040]设定判决门限S1，如果Err彡ε i，则处理过程结束，当前向量P中的特征参数值即为求解的最终结果。反之，如果Err > S1，则进入步骤(4)。
[0041]门限ε ,的取值决定了所提取的特征参数值的精准度，同时影响处理过程的迭代次数。^取值越小，特征参数值的精准度就越高，处理所需迭代次数越多。需要注意的是，如果^取值过小，则该迭代过程可能将无法最终收敛。反之，^取值越大，所提取的参数精准度将相应降低，而迭代次数将会减少。
[0042]步骤(4)根据当前特征参数向量P，构造矩阵J。
￠/.(/0丨,？? ￠/ (/^7,,卩)？?
(7/^5/7,
￠7/ (//7,,卩)0/ (1)1,,0/()11-,,卩)
[0043]I = 0^01？25/7,
0)趴坩',飞)
— 印丨印20/1, —
[0044]步骤(5)计算每次迭代过程中，特征参数向量？的更新向量!I =[八巧八％八巧]7，八巧，八巧和八分别为特征参数巧，？2和的待定更新值。构造误差矢量2。
[0045]2 = (叫，户)，(叫，？；?，…0%，户)11
[0046]则:
[0047]！！ = [了?了+入父乜叫(了1)^)]—1※/※^
[0048]其中，(11^( ? )表示矩阵对角元素提取和创建对角阵操作。
[0049]步骤(6)计算更新向量！！的度量值9⑶。
[[/(历,.，^)-4] -[[/(?，妒十0)-式]
[0050]即=^~-~,-
[0051]步骤(7)更新特征参数向量？和迭代步长参数\。设定判决门限、，如果更新向量0的度量值0⑶? 8 2，则当前特征参数向量？数值由？+?替代，即？ 一 9+?，完成更新，同时当前迭代步长参数\数值减小至^/1(,即\ \ II反之，如果0 (?) ( 82，则当前特征参数向量？保持不变，同时迭代步长参数、数值增加X倍，即、—XXX。X为比例因子，一般取值范围为5?20。判决门限8 2应根据样点偏差程度、收敛速度要求等具体指标设置合适数值。
[0052]在完成特征参数向量？和迭代步长参数\更新后，返回至步骤(3),进入行下一轮迭代。
【权利要求】
1.一种基于非线性拟合方式的肽质谱峰特征参数提取方法，其特征在于: 设质谱图中某离子的高斯峰由N个样点组成，N ^ 3 ;对样点按其丰度值从大到小排序后，其坐标构成集合A ;
A = {(HI1, (I1)，(m2, d2)，…(mN, dN)} 其中，Hii表示质荷比，Cli表示丰度值，i e {1，2，3，…，N};准备通过样点拟合出的高斯曲线其函数形式设为:
Jfzft V
f(x, P) = p, xe ' Pi J 其中，函数f(x，P)代表理论丰度值，自变量X代表质荷比，PpP2和P3S待求解的高斯曲线特征参数，分别表征缩放因子、质心、标准差，构成特征参数向量P = [P1 P2 P3I ；具体步骤如下: 步骤(I)根据丰度值最大的3个样点数据，对高斯曲线特征参数赋初值；
-f -P1 X? η>2~Ρι[ ^%-Ρι Γ
ρ, = —χ d, xe^ Pi J -\-d0xe^ Pi」+d,xeK ih」.1 323
_ I [ln(i/2)- ln(i/.)] Xinf + [ln(i/.)- ln(i/,)] X m; +[lnii/,) - ln(i/.)] x m:
P1 2 [ln(i1.)-ln(i/.)]x//z, + [ln(i/,) -1n(i/,)] x nu + ) - ln(i/,)] x ///,
=Ix ?θ?2-Ρ:? -(ln1-PzY.Pl)2 -('nZ- PlY . 3 ~2 ]j ln(i/,)-ln(i/2) ^ In(J2)-ln(i/;) 其中，ln(.)表示取自然对数操作； 步骤(2)选择合适数值初始化迭代步长参数λ，该参数初始化数值的大小将影响迭代次数和收敛速度； 步骤(3)计算拟合结果误差Err，判定迭代过程是否结束； ["=Σ[/(3，p)-4]2
/=1 设定判决门限S1，如果ErrS S1，则处理过程结束，当前向量P中的特征参数值即为求解的最终结果；反之，如果Err > ε i，则进入步骤(4)； 步骤(4)根据当前特征参数向量P，构造矩阵J ;
CfQni, P) Cfjniy, P) Cfjmi, P)
Φ,op2cpy
Cfim1, P) cf (nu, P) cf(nu, P) J = φ,Gp2dp.,
cf{ms, P) cf(m,, P) cf(m,, Ρ)
Gplδρ--ρ.步骤(5)计算每次迭代过程中，特征参数向量P的更新向量H= [ΛΡι Ap2 Δρ3]τ,ΔΡι，ΔΡ2和Δ P3分别为特征参数Pi，P2和Ρ3的待定更新值；构造误差矢量E ;
E = [dff (Hi1, P)，d2-f (m2, P)，…dN-f (mN, P) ]τ
则: H= [JTXJ+A Xdiag(JTXJ) F1XJtXE 其中，diag(.)表示矩阵对角元素提取和创建对角阵操作； 步骤(6)计算更新向量H的度量值P⑶；
Σ[/(";,, P)-c-/.]:-1[/(m, P + H)-^.]2 p(H) = ^?~?ξΙ--
2HTx(/lxH + JTxE) 步骤(7)更新特征参数向量P和迭代步长参数λ ;设定判决门限ε2，如果更新向量H的度量值P⑶> ε2，则当前特征参数向量P数值由Ρ+Η替代，即P —Ρ+Η，完成更新，同时当前迭代步长参数λ数值减小至λ/K，S卩λ — λ/K ;反之，如果P (H) ( ε2，则当前特征参数向量P保持不变，同时迭代步长参数λ数值增加K倍，S卩λ — ΚΧλ ;Κ为比例因子，取值范围为5?20 ;在完成特征参数向量P和迭代步长参数λ更新后，返回至步骤(3)，进行下一轮迭代。
【文档编号】G01N27/62GK104316591SQ201410498854
【公开日】2015年1月28日 申请日期:2014年9月25日 优先权日:2014年9月25日
【发明者】易志强, 李芸, 章剑秋, 曾嵘, 姚英彪, 张福洪, 李希元 申请人:杭州电子科技大学