CN104316591A - Non-linear fitting mode-based peptide mass spectrum speak characteristic parameter extraction method - Google Patents

Abstract

The invention relates to a non-linear fitting mode-based peptide mass spectrum speak characteristic parameter extraction method. The prior art is difficult to guarantee a precision of the extracted mass spectrum speak characteristic parameter when distributing error of sample points formed in a peptide mass spectrum is large. The non-linear fitting mode-based peptide mass spectrum speak characteristic parameter extraction method is characterized in that data of multiple sample points is utilized, the difference minimum of actual data and a fitting result is used as a guide, and a characteristic parameter estimated value is continuously updated by an iteration method until convergence conditions are satisfies so that the final characteristic parameter estimated value is obtained. The method provided by the invention can effectively reduce the unfavorable influence produced by sample point distributing error on Gauss curve characteristic parameter solving, improve a characteristic parameter value precision and is conducive to peptide segment identification precision improvement.

Description

A kind of peptide mass spectra peak characteristic parameter extraction method based on nonlinear fitting mode

Technical field

The invention belongs to biological mass spectrometry data prediction and information extraction technology field, be specifically related to a kind of peptide mass spectra peak characteristic parameter extraction method based on nonlinear fitting mode.

Background technology

Peptide qualification based on tandem mass spectrum is widely used technology in current proteome research field.Peptide to be identified in a mass spectrometer by cracked be fragmention, thus generate tandem mass spectrum data, and compare with theoretical tandem mass spectra storehouse or the peptide section mass spectral database identified and analyze, finally complete the qualification to unknown peptide section.

Carry out Mass Spectrometer Method to certain ion under normal circumstances, detected mass-to-charge ratio data are not single numerical points, but there is some sampling points, and on mass spectrogram, it fits to Gaussian curve, i.e. Gaussian peak.For determining the specific charge of this ion, pre-service need be carried out to these sampling points, calculating the barycenter (Centroid) in its X direction, i.e. the actual measurement mass-to-charge ratio of this ion.According to required barycenter, other characteristic parameters such as the maximum Abundances of this ion and then can be extrapolated.Current barycenter method for solving has multiple, more common thinking is: assuming that each sampling point mass spectrogram being formed Gaussian peak is all strictly distributed on certain Gaussian curve, utilize the numerical value (mass-to-charge ratio and Abundances) of each sampling point, be updated in the common Gaussian curvilinear function expression formula of unknown parameters, structure Simultaneous Equations, thus solve the characteristic parameter of corresponding Gaussian peak, comprise barycenter, maximum Abundances etc.A proteomics data analysis software MAXQUANT adopts current application very widely is namely this method.But in reality detects, by the impact of the factors such as experiment condition, place environment and instrument and equipment noise, on mass spectrogram, each sampling point often and non-critical is distributed on Gaussian curve, but there is certain deviation.When each sampling point amount of deflection is comparatively large, then the assumed condition in said method is difficult to set up, and the characteristic parameter solved thus certainly will be caused numerically to there is comparatively big error, and then has influence on the precision of peptide section qualification.

Summary of the invention

The object of the invention is to the shortcoming and defect solving said method, propose a kind of peptide mass spectra peak characteristic parameter extraction method based on nonlinear fitting mode.

If in mass spectrogram, the Gaussian peak of certain ion is made up of N number of sampling point, N >=3 under normal circumstances.After sorting from big to small by its Abundances to sampling point, its coordinate forms set A.

A＝{(m ₁,d ₁),(m ₂,d ₂),…(m _N,d _N)}

Wherein, m _irepresent mass-to-charge ratio, d _irepresent abundance, its value be greater than 0, i ∈ 1,2,3 ..., N}.Its functional form of Gaussian curve gone out by spot fitting is needed to be set to:

$f(x,P)=p_1\×e^{-{\left(\frac{x-p_2}{p_3}\right)}^2}$

Wherein, function f (x, P) represents Abundances, and independent variable x represents mass-to-charge ratio, p ₁, p ₂and p ₃for Gaussian curve characteristic parameter to be solved, characterize zoom factor, barycenter, standard deviation respectively, constitutive characteristic parameter vector P=[p ₁p ₂p ₃].Described characteristic parameter extraction method treatment step is as follows:

Step (1) according to maximum 3 the sampling point data of Abundances, to Gaussian curve characteristic parameter initialize.

$p_1=\frac{1}{3}\×[d_1\×e^{{\left(\frac{m_1-p_2}{p_3}\right)}^2}+d_2\×e^{{\left(\frac{m_2-p_2}{p_3}\right)}^2}+d_3\×e^{{\left(\frac{m_3-p_2}{p_3}\right)}^2}]$

$p_2=\frac{1}{2}\×\frac{[\ln \left(d_2\right)-\ln \left(d_3\right)]\×m_1^2+[\ln \left(d_3\right)-\ln \left(d_1\right)]\×m_2^2+[\ln \left(d_1\right)-\ln \left(d_2\right)]\×m_3^2}{[\ln \left(d_2\right)-\ln \left(d_3\right)]\×m_1+[\ln \left(d_3\right)-\ln \left(d_1\right)]\×m_2+[\ln \left(d_1\right)-\ln \left(d_2\right)]\×m_3}$

$p_3=\frac{1}{2}\×[\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)}}]$

Wherein, right log operations is taken from ln () expression.

Step (2) selects appropriate value initialization iteration step length parameter lambda, and the large young pathbreaker of this parameter initialization numerical value affects iterations and the speed of convergence for the treatment of method.

Step (3) digital simulation resultant error Err, judges whether iterative process terminates.

$\mathrm{Err}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{[f(m_i,P)-d_i]}^2$

Setting decision threshold ε ₁if, Err≤ε ₁, then processing procedure terminates, and the characteristic ginseng value in current vectorial P is the net result solved.Otherwise, if Err > is ε ₁, then step (4) is entered.

Thresholding ε ₁value determine the precision of extracted characteristic ginseng value, affect the iterations of processing procedure simultaneously.ε ₁value is less, and the precision of characteristic ginseng value is higher, and the required iterations of process is more.It should be noted that if ε ₁value is too small, then this iterative process may will cannot finally restrain.Otherwise, ε ₁value is larger, and the parameter precision extracted is by corresponding reduction, and iterations will reduce.

Step (4) is according to current signature parameter vector P, structural matrix J.

$J=\left[\begin{array}{ccc}\frac{\∂f(m_1,P)}{{\∂p}_1}& \frac{\∂f(m_1,P)}{{\∂p}_2}& \frac{\∂f(m_1,P)}{{\∂p}_3}\\ \frac{\∂f(m_2,P)}{{\∂p}_1}& \frac{\∂f(m_2,P)}{{\∂p}_2}& \frac{\∂f(m_2,P)}{{\∂p}_3}\\ .& .& .\\ .& .& .\\ .& .& .\\ \frac{\∂f(m_N,P)}{{\∂p}_1}& \frac{\∂f(m_N,P)}{{\∂p}_2}& \frac{\∂f(m_N,P)}{{\∂p}_3}\end{array}\right]$

Step (5) calculates in each iterative process, renewal vector H=[the Δ p of characteristic parameter vector P ₁Δ p ₂Δ p ₃] ^t, Δ p ₁, Δ p ₂with Δ p ₃be respectively characteristic parameter p ₁, p ₂and p ₃updated value undetermined.Instrument error vector E.

E＝[d ₁-f(m ₁,P),d ₂-f(m ₂,P),…d _N-f(m _N,P)] ^T

Then:

H＝[J ^T×J+λ×diag(J ^T×J)] ^-1×J ^T×E

Wherein, diag () representing matrix diagonal element extracts and creates diagonal matrix operation.

Step (6) calculates the metric ρ (H) upgrading vectorial H.

$\mathrm{\ρ}\left(H\right)=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{[f(m_i,P)-d_i]}^2-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{[f(m_i,P+H)-d_i]}^2}{{2H}^T\×(\mathrm{\λ}\×H+J^T\×E)}$

Step (7) regeneration characteristics parameter vector P and iteration step length parameter lambda.Setting decision threshold ε ₂if upgrade metric ρ (H) the > ε of vectorial H ₂, then current signature parameter vector P numerical value is substituted by P+H, i.e. P ← P+H, completes renewal, and current iteration step parameter λ numerical value is decreased to λ/K, i.e. λ ← λ/K simultaneously.Otherwise, if ρ (H)≤ε ₂, then current signature parameter vector P remains unchanged, and Simultaneous Iteration step parameter λ numerical value increases K doubly, i.e. λ ← K × λ.K is scale factor, and general span is 5 ~ 20.Decision threshold ε ₂appropriate value should be set according to specific targets such as sampling point extent of deviation, speed of convergence requirements.

After completing characteristic parameter vector P and the renewal of iteration step length parameter lambda, be back to step (3), enter row next round iteration.

Peptide mass spectra peak characteristic parameter extraction method in the present invention, adopts various some nonlinear fitting mode to solve characteristic parameter, decreases the adverse effect that sampling point distribution bias is brought, improve parameter extraction precision, and then is conducive to the improvement of peptide section qualification precision.

Embodiment

Step (1) according to maximum 3 the sampling point data of Abundances, to Gaussian curve characteristic parameter initialize.

$p_1=\frac{1}{3}\×[d_1\×e^{{\left(\frac{m_1-p_2}{p_3}\right)}^2}+d_2\×e^{{\left(\frac{m_2-p_2}{p_3}\right)}^2}+d_3\×e^{{\left(\frac{m_3-p_2}{p_3}\right)}^2}]$

$p_2=\frac{1}{2}\×\frac{[\ln \left(d_2\right)-\ln \left(d_3\right)]\×m_1^2+[\ln \left(d_3\right)-\ln \left(d_1\right)]\×m_2^2+[\ln \left(d_1\right)-\ln \left(d_2\right)]\×m_3^2}{[\ln \left(d_2\right)-\ln \left(d_3\right)]\×m_1+[\ln \left(d_3\right)-\ln \left(d_1\right)]\×m_2+[\ln \left(d_1\right)-\ln \left(d_2\right)]\×m_3}$

$p_3=\frac{1}{2}\×[\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)}}]$

Wherein, right log operations is taken from ln () expression.

Step (2) selects appropriate value initialization iteration step length parameter lambda, and the large young pathbreaker of this parameter initialization numerical value affects iterations and the speed of convergence for the treatment of method.

Step (3) digital simulation resultant error Err, judges whether iterative process terminates.

$\mathrm{Err}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{[f(m_i,P)-d_i]}^2$

Setting decision threshold ε ₁if, Err≤ε ₁, then processing procedure terminates, and the characteristic ginseng value in current vectorial P is the net result solved.Otherwise, if Err > is ε ₁, then step (4) is entered.

Thresholding ε ₁value determine the precision of extracted characteristic ginseng value, affect the iterations of processing procedure simultaneously.ε ₁value is less, and the precision of characteristic ginseng value is higher, and the required iterations of process is more.It should be noted that if ε ₁value is too small, then this iterative process may will cannot finally restrain.Otherwise, ε ₁value is larger, and the parameter precision extracted is by corresponding reduction, and iterations will reduce.

Step (4) is according to current signature parameter vector P, structural matrix J.

$J=\left[\begin{array}{ccc}\frac{\∂f(m_1,P)}{{\∂p}_1}& \frac{\∂f(m_1,P)}{{\∂p}_2}& \frac{\∂f(m_1,P)}{{\∂p}_3}\\ \frac{\∂f(m_2,P)}{{\∂p}_1}& \frac{\∂f(m_2,P)}{{\∂p}_2}& \frac{\∂f(m_2,P)}{{\∂p}_3}\\ .& .& .\\ .& .& .\\ .& .& .\\ \frac{\∂f(m_N,P)}{{\∂p}_1}& \frac{\∂f(m_N,P)}{{\∂p}_2}& \frac{\∂f(m_N,P)}{{\∂p}_3}\end{array}\right]$

Step (5) calculates in each iterative process, renewal vector H=[the Δ p of characteristic parameter vector P ₁Δ p ₂Δ p ₃] ^t, Δ p ₁, Δ p ₂with Δ p ₃be respectively characteristic parameter p ₁, p ₂and p ₃updated value undetermined.Instrument error vector E.

E＝[d ₁-f(m ₁,P),d ₂-f(m ₂,P),…d _N-f(m _N,P)] ^T

Then:

H＝[J ^T×J+λ×diag(J ^T×J)] ^-1×J ^T×E

Wherein, diag () representing matrix diagonal element extracts and creates diagonal matrix operation.

Step (6) calculates the metric ρ (H) upgrading vectorial H.

$\mathrm{\ρ}\left(H\right)=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{[f(m_i,P)-d_i]}^2-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{[f(m_i,P+H)-d_i]}^2}{{2H}^T\×(\mathrm{\λ}\×H+J^T\×E)}$

Step (7) regeneration characteristics parameter vector P and iteration step length parameter lambda.Setting decision threshold ε ₂if upgrade metric ρ (H) the > ε of vectorial H ₂, then current signature parameter vector P numerical value is substituted by P+H, i.e. P ← P+H, completes renewal, and current iteration step parameter λ numerical value is decreased to λ/K, i.e. λ ← λ/K simultaneously.Otherwise, if ρ (H)≤ε ₂, then current signature parameter vector P remains unchanged, and Simultaneous Iteration step parameter λ numerical value increases K doubly, i.e. λ ← K × λ.K is scale factor, and general span is 5 ~ 20.Decision threshold ε ₂appropriate value should be set according to specific targets such as sampling point extent of deviation, speed of convergence requirements.

After completing characteristic parameter vector P and the renewal of iteration step length parameter lambda, be back to step (3), enter row next round iteration.

Claims (1)

1., based on a peptide mass spectra peak characteristic parameter extraction method for nonlinear fitting mode, it is characterized in that:

If in mass spectrogram, the Gaussian peak of certain ion is made up of N number of sampling point, N >=3; After sorting from big to small by its Abundances to sampling point, its coordinate forms set A;

A＝{(m ₁,d ₁),(m ₂,d ₂),…(m _N,d _N)}

Wherein, m _irepresent mass-to-charge ratio, d _irepresent Abundances, i ∈ 1,2,3 ..., N}; Its functional form of Gaussian curve preparing to be gone out by spot fitting is set to:

$f(x,P)=p_1\×e^{-{\left(\frac{x-p_2}{p_3}\right)}^2}$

Wherein, function f (x, P) representation theory Abundances, independent variable x represents mass-to-charge ratio, p ₁, p ₂and p ₃for Gaussian curve characteristic parameter to be solved, characterize zoom factor, barycenter, standard deviation respectively, constitutive characteristic parameter vector P=[p ₁p ₂p ₃];

Concrete steps are as follows:

Step (1) according to maximum 3 the sampling point data of Abundances, to Gaussian curve characteristic parameter initialize;

$p_1=\frac{1}{3}\×[d_1\×e^{{\left(\frac{m_1-p_2}{p_3}\right)}^2}+d_2\×e^{{\left(\frac{m_2-p_2}{p_3}\right)}^2}+d_3\×e^{{\left(\frac{m_3-p_2}{p_3}\right)}^2}]$

$p_2=\frac{1}{2}\×\frac{[\ln \left(d_2\right)-\ln \left(d_3\right)]\×m_1^2+[\ln \left(d_3\right)-\ln \left(d_1\right)]\×m_2^2+[\ln \left(d_1\right)-\ln \left(d_2\right)]\×m_3^2}{[\ln \left(d_2\right)-\ln \left(d_3\right)]\×m_1+[\ln \left(d_3\right)-\ln \left(d_1\right)]\×m_2+[\ln \left(d_1\right)-\ln \left(d_2\right)]\×m_3}$

$p_3=\frac{1}{2}\×[\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)}}]$

Wherein, right log operations is taken from ln () expression;

Step (2) selects appropriate value initialization iteration step length parameter lambda, and the large young pathbreaker of this parameter initialization numerical value affects iterations and speed of convergence;

Step (3) digital simulation resultant error Err, judges whether iterative process terminates;

$\mathrm{Err}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{[f(m_i,P)-d_i]}^2$

Setting decision threshold ε ₁if, Err≤ε ₁, then processing procedure terminates, and the characteristic ginseng value in current vectorial P is the net result solved; Otherwise, if Err > is ε ₁, then step (4) is entered;

Step (4) is according to current signature parameter vector P, structural matrix J;

$J=\left[\begin{array}{ccc}\frac{\∂f(m_1,P)}{{\∂p}_1}& \frac{\∂f(m_1,P)}{{\∂p}_2}& \frac{\∂f(m_1,P)}{{\∂p}_3}\\ \frac{\∂f(m_2,P)}{{\∂p}_1}& \frac{\∂f(m_2,P)}{{\∂p}_2}& \frac{\∂f(m_2,P)}{{\∂p}_3}\\ .& .& .\\ .& .& .\\ .& .& .\\ \frac{\∂f(m_N,P)}{{\∂p}_1}& \frac{\∂f(m_N,P)}{{\∂p}_2}& \frac{\∂f(m_N,P)}{{\∂p}_3}\end{array}\right]$

Step (5) calculates in each iterative process, renewal vector H=[the Δ p of characteristic parameter vector P ₁Δ p ₂Δ p ₃] ^t, Δ p ₁, Δ p ₂with Δ p ₃be respectively characteristic parameter p ₁, p ₂and p ₃updated value undetermined; Instrument error vector E;

E＝[d ₁-f(m ₁,P),d ₂-f(m ₂,P),…d _N-f(m _N,P)] ^T

Then:

H＝[J ^T×J+λ×diag(J ^T×J)] ^-1×J ^T×E

Wherein, diag () representing matrix diagonal element extracts and creates diagonal matrix operation;

Step (6) calculates the metric ρ (H) upgrading vectorial H;

$\mathrm{\ρ}\left(H\right)=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{[f(m_i,P)-d_i]}^2-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{[f(m_i,P+H)-d_i]}^2}{{2H}^T\×(\mathrm{\λ}\×H+J^T\×E)}$

Step (7) regeneration characteristics parameter vector P and iteration step length parameter lambda; Setting decision threshold ε ₂if upgrade metric ρ (H) the > ε of vectorial H ₂, then current signature parameter vector P numerical value is substituted by P+H, i.e. P ← P+H, completes renewal, and current iteration step parameter λ numerical value is decreased to λ/K, i.e. λ ← λ/K simultaneously; Otherwise, if ρ (H)≤ε ₂, then current signature parameter vector P remains unchanged, and Simultaneous Iteration step parameter λ numerical value increases K doubly, i.e. λ ← K × λ; K is scale factor, and span is 5 ~ 20; After completing characteristic parameter vector P and the renewal of iteration step length parameter lambda, be back to step (3), carry out next round iteration.