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

Abstract

The present invention relates to a kind of peptide mass spectra peak characteristic parameter extraction method.When existing method exists relatively large deviation for its distribution of each sampling point forming spectral peak in peptide fragment mass spectrogram, there is the deficiency being difficult to ensure that extracted mass spectra peak characteristic parameter precision.The present invention proposes peptide mass spectra peak characteristic parameter extraction method based on nonlinear fitting mode, utilize multiple sampling point data, with the minimum guiding of difference between real data and fitting result, alternative manner is used to constantly update time parameters estimation value, until meeting the condition of convergence, thus obtain final characteristic parameter valuation.The method is effectively reduced sampling point distribution bias and Gaussian curve characteristic parameter solves the adverse effect brought, and improves characteristic parameter numerical value precision, and then beneficially peptide fragment identifies the improvement of precision.

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

It is widely used in current proteome research field that peptide based on tandem mass spectrum is identified Technology.Peptide to be identified is fractured as fragment ion in a mass spectrometer, thus generates tandem mass spectrum Data, and compare with theoretical tandem mass spectra storehouse or the peptide fragment mass spectral database identified and analyze, Finally complete the qualification to unknown peptide fragment.

Certain ion carries out Mass Spectrometer Method under normal circumstances, and detected mass-to-charge ratio data are not Being single numerical point, but there is some sampling points, on mass spectrogram, it fits to Gaussian curve, I.e. Gaussian peak.For determining the charge-mass ratio of this ion, these sampling points need to be pre-processed, calculate Go out the barycenter (Centroid) in its X direction, i.e. the actual measurement mass-to-charge ratio of this ion.According to institute Seek barycenter, other characteristic parameters such as this ion maximum Abundances and then can be extrapolated.Barycenter at present Method for solving has multiple, and relatively common thinking is: assuming that constitute each of Gaussian peak on mass spectrogram Individual sampling point is the most strictly distributed on certain Gaussian curve, utilize each sampling point numerical value (mass-to-charge ratio and Abundances), it is updated in the common Gaussian curvilinear function expression formula of unknown parameters, constructs simultaneous Equation group, thus solve the characteristic parameter of corresponding Gaussian peak, including barycenter, maximum Abundances etc.. The extremely wide a proteomics data of current application is analyzed software MAXQUANT and is adopted Be i.e. this method.But in actually detected, by experiment condition, place environment and The impact of the factors such as instrument and equipment noise, on mass spectrogram, each sampling point is often and non-critical is distributed in On Gaussian curve, but there is certain deviation.When each sampling point amount of deflection is relatively big, the most above-mentioned Assumed condition in method is difficult to set up, thus the characteristic parameter solved certainly will be caused at numerical value The bigger error of upper existence, and then have influence on the precision that peptide fragment is identified.

Summary of the invention

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

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

A={ (m₁,d₁),(m₂,d₂),…(m_N,d_N)}

Wherein, m_iRepresent mass-to-charge ratio, d_iRepresent abundance, its value be more than 0, i ∈ 1,2 ..., N}.Need Its functional form of Gaussian curve 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) represents Abundances, and independent variable x represents mass-to-charge ratio, p₁、p₂And p₃For Gaussian curve characteristic parameter to be solved, characterizes zoom factor, barycenter, standard deviation, structure respectively Become characteristic parameter vector P=[p₁p₂p₃].Described characteristic parameter extraction method processes step such as Under:

Gaussian curve feature, according to 3 sampling point data of Abundances maximum, is joined by step (1) Number composes initial value.

$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 to initialize iteration step length parameter lambda, this parameter initialization The size of numerical value will affect iterations and the convergence rate for the treatment of method.

Step (3) digital simulation resultant error Err, it is determined that whether iterative process terminates.

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

Set decision threshold ε₁If, Err≤ε₁, then processing procedure terminates, the spy in current vector P Levy parameter value and be the final result solved.Whereas if Err ＞ ε₁, then step (4) is entered.

Thresholding ε₁Value determine the precision of extracted characteristic ginseng value, affect place simultaneously The iterations of reason process.ε₁Value is the least, and the precision of characteristic ginseng value is the highest, processes Required iterations is the most.It should be noted that if ε₁Value is too small, then this iteration Process may will be unable to finally restrain.Otherwise, ε₁Value is the biggest, the parameter precision extracted To reduce accordingly, 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, the renewal vector of characteristic parameter vector P H=[Δ p₁ Δp₂ Δp₃]^T, Δ p₁, Δ p₂With Δ p₃It is respectively characteristic parameter p₁, p₂And p₃Treat Determine updated value.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) updating vector H.

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

Step (7) updates characteristic parameter vector P and iteration step length parameter lambda.Set decision threshold ε₂If updating metric ρ (H) the ＞ ε of vector H₂, then current signature parameter vector P numerical value Substituted by P+H, i.e. P ← P+H, complete to update, current iteration step parameter λ numerical value simultaneously It is decreased to λ/K, i.e. λ ← λ/K.Whereas if ρ (H)≤ε₂, then current signature parameter to Amount P keeps constant, and Simultaneous Iteration step parameter λ numerical value increases K times, i.e. λ ← K × λ.K is Scale factor, general span is 5～20.Decision threshold ε₂Should according to sampling point extent of deviation, The specific targets such as convergence rate requirement arrange appropriate value.

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

Peptide mass spectra peak characteristic parameter extraction method in the present invention, uses many sampling points Nonlinear Quasi Conjunction mode solves characteristic parameter, decreases the adverse effect that sampling point distribution bias is brought, carries Rise parameter extraction precision, and then beneficially peptide fragment has identified the improvement of precision.

Detailed description of the invention

A kind of peptide mass spectra peak characteristic parameter extraction method based on nonlinear fitting mode, Specific as follows:

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

A={ (m₁,d₁),(m₂,d₂),…(m_N,d_N)}

Wherein, m_iRepresent mass-to-charge ratio, d_iRepresent abundance, its value be more than 0, i ∈ 1,2 ..., N}.Need Its functional form of Gaussian curve 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) represents Abundances, and independent variable x represents mass-to-charge ratio, p₁、p₂And p₃For Gaussian curve characteristic parameter to be solved, characterizes zoom factor, barycenter, standard deviation, structure respectively Become characteristic parameter vector P=[p₁p₂p₃].It is as follows that characteristic parameter extraction method processes step:

Gaussian curve feature, according to 3 sampling point data of Abundances maximum, is joined by step (1) Number composes initial value.

$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 to initialize iteration step length parameter lambda, this parameter initialization The size of numerical value will affect iterations and the convergence rate for the treatment of method.

Step (3) digital simulation resultant error Err, it is determined that whether iterative process terminates.

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

Set decision threshold ε₁If, Err≤ε₁, then processing procedure terminates, the spy in current vector P Levy parameter value and be the final result solved.Whereas if Err ＞ ε₁, then step (4) is entered.

Thresholding ε₁Value determine the precision of extracted characteristic ginseng value, affect place simultaneously The iterations of reason process.ε₁Value is the least, and the precision of characteristic ginseng value is the highest, processes Required iterations is the most.It should be noted that if ε₁Value is too small, then this iteration Process may will be unable to finally restrain.Otherwise, ε₁Value is the biggest, the parameter precision extracted To reduce accordingly, 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, the renewal vector of characteristic parameter vector P H=[Δ p₁ Δp₂ Δp₃]^T, Δ p₁, Δ p₂With Δ p₃It is respectively characteristic parameter p₁, p₂And p₃Treat Determine updated value.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) updating vector H.

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

Step (7) updates characteristic parameter vector P and iteration step length parameter lambda.Set decision threshold ε₂If updating metric ρ (H) the ＞ ε of vector H₂, then current signature parameter vector P numerical value Substituted by P+H, i.e. P ← P+H, complete to update, current iteration step parameter λ numerical value simultaneously It is decreased to λ/K, i.e. λ ← λ/K.Whereas if ρ (H)≤ε₂, then current signature parameter to Amount P keeps constant, and Simultaneous Iteration step parameter λ numerical value increases K times, i.e. λ ← K × λ.K is Scale factor, general span is 5～20.Decision threshold ε₂Should according to sampling point extent of deviation, The specific targets such as convergence rate requirement arrange appropriate value.

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

Claims (1)

1. a peptide mass spectra peak characteristic parameter extraction method based on nonlinear fitting mode, its feature exists In:

If in mass spectrogram, the Gaussian peak of certain ion is made up of N number of sampling point, N >=3；To sampling point by its abundance After value sequence from big to small, its coordinate constitutes set A；

A={ (m₁,d₁),(m₂,d₂),…(m_N,d_N)}

Wherein, m_iRepresent mass-to-charge ratio, d_iExpression Abundances, i ∈ 1,2 ..., N}；It is ready to pass through sampling point to intend Its functional form of the Gaussian curve closed out 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₂With p₃For Gaussian curve characteristic parameter to be solved, characterize zoom factor, barycenter, standard deviation, structure respectively Become characteristic parameter vector P=[p₁p₂p₃]；

Specifically comprise the following steps that

Step (1) is according to 3 sampling point data of Abundances maximum, at the beginning of composing Gaussian curve characteristic parameter Value；

$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{{\left(m_2-p_2\right)}^2-{\left(m_1-p_2\right)}^2}{\ln \left(d_1\right)-\ln \left(d_2\right)}}+\sqrt{\frac{{\left(m_3-p_2\right)}^2-{\left(m_2-p_2\right)}^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 to initialize iteration step length parameter lambda, this parameter initialization numerical value Size will affect iterations and convergence rate；

Step (3) digital simulation resultant error Err, it is determined that whether iterative process terminates；

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

Set decision threshold ε₁If, Err≤ε₁, then processing procedure terminates, the feature in current vector P Parameter value is the final result solved；Whereas if Err ＞ ε₁, 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\left(m_1,P\right)}{\∂p_1}& \frac{\∂f\left(m_1,P\right)}{\∂p_2}& \frac{\∂f\left(m_1,P\right)}{\∂p_3}\\ \frac{\∂f\left(m_2,P\right)}{\∂p_1}& \frac{\∂f\left(m_2,P\right)}{\∂p_2}& \frac{\∂f\left(m_2,P\right)}{\∂p_3}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ \frac{\∂f\left(m_N,P\right)}{\∂p_1}& \frac{\∂f\left(m_N,P\right)}{\∂p_2}& \frac{\∂f\left(m_N,P\right)}{\∂p_3}\end{array}\right]$

Step (5) calculates in each iterative process, the renewal vector of characteristic parameter vector P H=[Δ p₁ Δp₂ Δp₃]^T, Δ p₁, Δ p₂With Δ p₃It is respectively characteristic parameter p₁, p₂And p₃Renewal undetermined Value；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) updating vector H；

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

Step (7) updates characteristic parameter vector P and iteration step length parameter lambda；Set decision threshold ε₂, If updating metric ρ (H) the ＞ ε of vector H₂, then current signature parameter vector P numerical value is replaced by P+H Generation, i.e. P ← P+H, completing to update, current iteration step parameter λ numerical value is decreased to λ/K, i.e. simultaneously λ←λ/K；Whereas if ρ (H)≤ε₂, then current signature parameter vector P keeps constant, changes simultaneously Long parameter lambda numerical value of riding instead of walk increases K times, 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, it is back to step (3), next round iteration is carried out.