GB2405991A - Method of calibrating a mass spectrometer - Google Patents

Abstract

A method of calibrating a mass spectrometer is disclosed. The times of detection of a plurality of ions having known mass to charge ratios are measured. A data set comprising known mass to charge ratios and times of ion detection is generated. A prior probability distribution function is assigned for a flexibility parameter D , wherein the flexibility parameter D is a measure of curvature of the data set. The most probable flexibility parameter D consistent with the data set and with the prior probability distribution function is determined. Then, the most probable curve consistent with the data set and the most probable flexibility parameter D is determined.

Description

240599 1

MASS SPECTROMETER

The present invention relates to a method of calibrating a mass spectrometer.

It is known to calibrate a mass spectrometer by analysing ions having accurately known mass to charge ratios M1,M2... etc. and then experimentally measuring the mass to charge ratios X1,X2... etc. of the ions.

Ideally, each experimentally determined measurement X should be identical to the truth M. However, instrumentation is never perfect and in practice there will be a response function X=f(M) that will be different from the ideal X=M. It is known to model the response using a low-order polynomial or spline. These are smooth mathematical functions that can be fitted using some simple equations.

Calibrating the response of a mass spectrometer using a low-order polynomial or a spline requires a decision to be made concerning the number of coefficients to be used. If too few coefficients are used then the measurements cannot be fitted properly.

However, if too many coefficients are used then the measurements become over-fitted with small deviations due to noise being interpreted as truth. The number of coefficients used is decided upon using arbitrary rules of thumb and this can therefore introduce error into the calibration process.

Conventional calibration approaches also require the mass to charge ratios of calibration masses to be spaced reasonably uniformly apart. If two calibration masses have similar mass to charge ratios then there is a risk that their measurement errors will conspire to make the response function locally too steep or too - 2 - shallow. Once that has happened in one place then the error may then be propagated throughout the entire mass range.

By their nature, low-order polynomial or spline functions cannot properly describe the calibration accuracy. Clearly, the accuracy should be greatest at and near where precise measurements have been made and should deteriorate elsewhere. However, no matter how the measurement uncertainties are factored into the coefficients, a low-order function simply does not have the freedom to do this.

Low-order polynomials are also poor at mo.deling functions that are not low-order polynomials. For example, the characteristics of an instrument may well change such that the response function is different above and below a particular mass to charge ratio. A low-order polynomial is unable to model this very well.

It is therefore desired to provide an improved method of calibrating a mass spectrometer.

According to the present invention there is provided a method of calibrating a mass spectrometer comprising: measuring the times of detection of a plurality of ions having known mass to charge ratios; generating a data set comprising the known mass to charge ratios and the times of ion detection; assigning a prior probability distribution function for a flexibility parameter A, wherein the flexibility parameter is a measure of curvature of the data set; determining the most probable flexibility parameter consistent with the data set and with the prior probability distribution function; and - 3 - determining the most probable curve consistent with the data set and the most probable flexibility parameter A. The most probable curve is preferably used as a calibration of mass to charge ratio versus time of ion detection. An uncertainty in the calibration of mass to charge ratio versus time of ion detection is preferably determined.

The uncertainty in the calibration of mass to charge ratio versus time of ion detection is preferably determined from a measure of a spread of plausible calibration curves about the calibration curve.

The probability distribution function for the flexibility parameter is preferably proportional to ant wherein n is a constant. Preferably, 1 < n < 1.5, further preferably 1.125 < n < 1.375. According to a particularly preferred embodiment n is 1.25.

The most probable flexibility parameter is preferably that for which a joint probability of the data set and of the probability distribution function for the flexibility parameter is maximized.

The most probable curve is preferably that for which a joint probability of the data set and of a calibration of mass to charge ratio versus time of ion detection is maximized.

The joint probability of the data set and of a calibration of mass to charge ratio versus time of ion detection is preferably that for which the probability of the data set is that given by the most probable flexibility parameter A. In order to calibrate mass to charge ratio measurements from a sequence of calibration peaks at known "lock masses" the preferred embodiment uses 1 1 - 4 - probabilistic (Bayesian) analysis to assign a measure of calibration curvature and to average over the (extremely large) number of possible calibration functions having a similar degree of curvature. This provides a smooth but otherwise free-form calibration curve augmented with a known degree of reliability at every point on the axis.

The preferred embodiment therefore relates to a new method of performing mass calibration. Instead of restricting the response function f(M) to a specific low-order form, the preferred embodiment allows it to be a freeform curve chosen to be as straight as reasonably possible i.e. having plausibly low overall curvature.

Specifically, this is done by assigning to any candidate responsefunction an intrinsic or prior probability which decreases with its curvature according to Gaussian statistics. Each calibration measurement, with its own uncertainty if any, introduces a likelihood factor which pulls the family of plausible curves towards it according to the standard rules of probability. The preferred method yields a central average calibration curve augmented by the local uncertainty everywhere.

The preferred embodiment exhibits a number of advantages compared with conventional methods of mass calibration. Firstly and importantly, a user is not required to assign or decide on a particular number of coefficients to be used to produce the calibration function. Instead, there is just a single degree of curvature parameter and that is given its most probable value. Accordingly, there is no possibility of user- error being introduced at this point. Secondly, all calibration information is used and calibration masses can be arbitrarily close together. In the limit - 5 - calibration masses may merge seamlessly into a single measurement. Thirdly, the treatment of accuracy is mathematically defensible. It is computed as a standard-deviation uncertainty at any required mass so that calibration uncertainties can be factored into any list of observed peaks. In particular, the accuracy may deteriorate away from the calibration points. Fourthly, a free-form curve is likely to be better at treating response functions that happen not to be low-order polynomials or splines.

Various embodiments of the present invention will now be described, by way of example only, and with reference to the accompanying drawings in which: Fig. 1 shows a graph of measured value y against true value xi Fig. 2 shows the offset between the measured value y and the true value x as a function of true value x; Fig. 3 shows a data set which was subsequently analyzed to provide comparative data; Fig. 4 shows a graph of the x,y data pairs from the data set presented in Fig. 3; Fig. 5 shows the residual error as a function of true value x after conventional least squares straight line fitting to the data presented in Fig. 3; and Fig. 6 shows the residual error as a function of true value x based on probabilistic calibration of the data presented in Fig. 3 according to a preferred embodiment.

The preferred method of calibrating a mass spectrometer will now be described in more detail.

Calibration involves transforming a measured coordinate towards true values. The starting point may be considered to comprise a set of measurements Y1Y2--- Yn giving apparent locations that ought to have been xx2...x. In mass spectrometry, the x values represent the calculated mass to charge ratios of accurately known reference masses and the y values represent the mass to charge ratio as experimentally measured for these reference masses using a mass analyzer.

In a Time of Flight (TOF) mass spectrometer the mass to charge ratio of an ion is proportional to the square of the time of flight of the ion through a field free flight region. By comparing the square root of the reference mass to charge ratio value to the measured or recorded flight time, a calibration expression relating these two values may then be generated.

Graphically, a suitable curve y(x) can be interpolated through the given points. Calibration of a measured value y may then be achieved by reading off the corresponding x value. However, accurate calibration will tend to fail due to ambiguity if the chosen curve turns downwards anywhere. In that case the data may be essentially unusable and will presumably be in serious error.

Fig. 1 shows a graph of measured y values against calculated or true values x. The vertical bars represent the precision of the measurement of the y values as standard deviation around the measured value.

Fig. 2 illustrates the same data wherein the y axis now represents the offset (or difference) between the measured y value and the true x value.

It is apparent from Fig. 2 that there are numerous different curves which could be fitted to the data.

Accordingly, the calibration procedures should allow for uncertainty even though the measured y data is itself l l completely accurate. For the following discussion the data points are referred to as knots.

Considering a calibration curve y(x) interpolated through a set of measured data values D, a natural measure of smoothness or lack of it is curvature y''.

The total curvature Q is given by: Q(y) = Jyr,(x)2 Ax ( 1) oo The total curvature Q may be promoted into a prior probability varying inversely with Q in order to discriminate against curvature. Specifically, a Gaussian prior may be assigned: Pr(y|)xexp(- Q/2O ( 2) Here Pr(y|)xexp(/20 is a flexibility parameter which defines the amount of curvature that may be reasonably expected. The notation Pra|b) denotes the conditional probability of "a given b is true". Estimation of from the data will be discussed later - for the following it will be treated as a given parameter.

The recorded data values D enter through a likelihood factor based on the chi squared misfit: Pr(Dly)xexp(-/ /2) ( 3) (y(Xk)-Dk) z2 = k 2 (4) ark 8 wherein: Y(Xk) Dk -Ark ( 5) for k = 1,2,,n (6) and wherein Ark iS the uncertainty on the data Dk Multiplying these two quantities yields the joint distribution: Pr(y,D|)xexp-(Q/2+12/2) (7) from which the best (i.e most probable) interpolant y, its uncertainty By and the likelihood Pr(D|) for any given value of can all be recovered.

To define the quantities more carefully the x axis may be digitised into a large number of grid points separated by arbitrarily small intervals. Positions thereby become large integers which are scaled back to ordinary values at the end. Local curvature at grid point i may be defined as: Hi = Yi+1 - 2Yi + Yi--l ( 8) Accordingly, the global curvature measure is in matrix formulation and may be given by: Q (Y) (y) = y Ay (9) - 9 - where Bold face is used for x-dimensional vectors and matrices. A is the 4th-derivative matrix of arbitrarily high dimension: 1+ -2 1 -2 5 -4 1 1 -4 6 -4 1 1 -4 6 -4 1 1 -4 6 A=I. . . 1 (10) 6 -4 1 1 -4 6 1 1 1 -4 6 -4 1 l -4 5 -2 1 -2 l+ The "1, -4, 6, -4, 1" internal rows and columns represent the 4th-derivative product of the "1, -2, 1" rows and columns of y''. Were to be 0, the corners of A would ensure that any straight line (y = a + b x) would be exactly annihilated. However, this choice makes A singular with two null eigenvalues. The correction E is imposed to bring det A back up to 1 for later convenience. Its value and the other changes it induces are trivially small if the range of the matrix is allowed to extend out towards x = m. With Q now defined in full, the properly normalized prior on interpolating functions y can be assigned as: Pr(y|)= det(A/2)exp(-y T Ay/2O (11) - 10 The corresponding normalised likelihood is: Pr(D|y)=Ziexp(-72/2) (12) where: / Z= (13) k X (Xk) k D) Aid =yT-2 y ( 14 k and wherein -2 iS the diagonal matrix which is zero everywhere except at the data points where its components are l/2k. In fact it does not matter whether the data and uncertainty D and are used in n dimensional or in x-dimensional context. Using the latter formulation throughout, the required joint probability is: Pr(y,DI() = Z-; det(A/2)exp- (yTAy/ +(y D)Ta-2(y - D)) (15) This is an important quantity in the preferred method of probabilistic analysis.

Most probable interpolating function The most probable interpolant minimises the quadratic expression given above in (15) above. The full optimising spline function y is given by: - 11 - y = D _ 2lJT (lJa. U + jV) UD (16) where U and V are tridiagonal matrices defined as: lJy=Vp Ok Y(Xk) Pk y (Xk) ( 17) With Pl and Pr, both zero, V is (n-2) x (n-2) symmetric, with indicies ranging from 2 to n-1. U is a (n-2) x n matrix operating on the e-dimensional y.

With regard (16) above, the knot values y are seen to converge upon the data D as the uncertainty shrinks and the computation proceeds to that limit without relying upon cancellation of large terms.

Likelihood of flexibility To estimate reference is made to the joint distribution (15) which is integrated over y to obtain the likelihood of : Z-ld7(Al / _)/ )exp-(YTBY/+(Y-D)Ta-2(Y-D)) Pr(DI) = idy Pr(y,DI+) (18) Z \|det(A+á-2) exp-2(D U (Up UT+ jV) IUD) - 12 The exponential is already calculable, in the order of n operations (where n is the number of calibration points) because the matrix to be inverted is pentadiagonal. Z may be computed and the determinant of A has earlier been set to 1. It remains to evaluate the m-dimensional determinant: -2 5 -4 1 1 -4 6-41 1 -46+a-4 1 1-46 -4 1 det(A + fa-2) = 1-4 6 -4 1 1 -4 6+b -4 1 1 -4 G -4 1 1 -4 6 -4 1 1 -4 5 -2 1 -2 1 ( 1 9) where a,b,... are the individual knot components 2,22, appearing sporadically down the leading diagonal. Determinants that terminate with 1, -4, 6, 4, 1 corners obey a recurrence relation. Let: ::: ::: ::: ::: P::: 1 ::: ::: ::: - 4 1 1 -4 6.

Ar = 1. . . 1 6 -4 1 1 -4 fi -4 1 1 - 4 G -4 1 -4 6 (20) - 13 be of dimension r. Direct expansion by the last row or column of Ar and its associated determinants yields the recurrence relation: Air -5Ar 1 + 106r 3 + 56r-4-r--S = 0 for which the general solution is quartic in r. Special cases, proved by induction from low dimension, are: 6 -4 1 -4 6 -4 1 1 -4 6 -4 1 1 - 4 6. = 132 (r + l)(r + 2)2(r + 3) 6 -41 1-4 6-4 1 1 -4 6 -4 1 -4 6 1 - 2 1 -2 6 -4 1 (21) 1 -4 fi -4 1-4 6 fi -4 1 1-4 6 -4 1 1 -4 6 -4 1 -4 6 (22) The latter shows why it was natural to assign det A = 1 earlier. More generally, Ar takes a quartic form in r outside and between the knots. The effect of a knot is therefore to change the quartic coefficients from one side to the other. If c' at position x is the rightmost knot element contained in Ar, then: ::: ::: ;:: ::: ::: ::: :: P;:: 1::: P;:: 1, ::: ::: ::: -4 1::: ;:: ::: -4 1... ...

1 -4 6+c -4 1 = 1 -4 6 -4 1 +c::: ::: :: l 1 -4 1 Q ^ Q (23) where Q is a 1, -4, 6, -4, 1 pentadiagonal extension.

Without the extra contribution from the knot, the initial quartic, which can be written as: Ar =a4(r-x)4 +a3(r-x)3 +a2(r-x)2 +a(r-x)+aO (24) would have extended from r x before the knot, straight through to r > x beyond the knot, as in the first determinant on the right. The second determinant on the right gives the increment to the quartic: p c 1... ... =c(detP detQ-detP_detQ-) ::: Q::: :: ::; ::: (25) where P_ is P without its last row and column, and Q is Q without its first row and column. However: detP=6x-l = a4 - a3 + a2 - al + aO ( 2 6) detP_=Ax 2 = 16a4 - 8a3 + 4a2 gal +aO (27) detQ=i (r -x+l)(r -x+2)2(r-x+3) (28) detQ=I (r - x)r - x + 1)2 (r - x + 2) (29) - 15 Hence, beyond a knot, the quartic form of fir is augmented by the corrective quartic: 12 c((a4 - a3 + a2 a,+aO)(r-x+l)(r-x+2)2(r-x+3) (16a4-8a3+4a2-2a+aO)(r-x)(r-x+l)2(r-x+2)) (30) Some of this detail drops out as the integer labels become indefinitely large on a fine grid. The partial determinants AL can be stepped past all n knots in the order of only n operations. The final knot AL will take its rightmost quartic form: /\r = a4(r - )4 + a3(r-X)3 + a2(r - [)2 + a-)+ aO (31) Direct expansion of more determinants shows that this terminates at the right-hand boundary as: d:: ::: 1 A -4 1 -4 6 -4 1 = 12a4.

1 -4 5 -2 1 -2 1 (32) So each term in the likelihood (18) of can be evaluated in just n operations.

Most probable flexibility As a consequence of the singularity or near singularity of A, Det(A+2) behaves almost like ? as - 16 0, so that the likelihood diverges as -1. Hence naive maximization of the likelihood would make extremely small, which amounts to forcing the interpolant to be the least squares straight line.

Accordingly, a prior can be set that favours flexibility by incorporating more than one power of which will overcome the awkward singularity in the evidence. On the other hand, at large the likelihood only goes down as "/2 in accord with the determinant. It is not desired for the prior to overcome this, lest maximization yield = which would require the interpolant to fit all the data exactly.

Since usable curvature information is available with as few as three data, a prior is desired that preferably increases with more slowly than 13/2. To obey both requirements an intermediate prior may be set: Pr(O =5'4 (33) can be determined by maximizing: i E()= f5/4 Pr(DIf)=5/iZ-! dt(A)-2) exp-(D'U7(u2uT+> (34) Uncertainty Reverting to the joint probability (15) written as: Pr(y,D|)= exp--(y-y)7 (A/±2)(y-y) (35) 17 it can be seen that the uncertainty by around the optimal interpolant y has covariance: < {yyT > = (A/ + a-27-1 (36) From this, the point-wise variance of the interpolant at any desired position x may be selected: <(:yX)2≥(AI±2)-i (37) Using the standard adjoint/determinant formula for matrix inverse leads to: < (y)2 > = jadj(A + ) (38) It has already been demonstrated how to evaluate the determinant in the denominator. Deleting the xth row and column of the matrix gives the adjoins component as: ::: P::: adj(A+á-2) 2z= ::: ::: ::. 1=detPdetQ-detP_detQ (39) where P is the entire matrix before x, including all leftward knots and Q is the entire matrix after x, including all rightward knots. However: detP=Ax-l detP_=Ax2 (40) which can be read off from the local quartic which was constructed as part of the calculation of the overall determinant.

A similar calculation, starting from +x and stepping backwards gives V, being the partial determinants calculated from the right, giving: detQ = Vr+l7 detQ- = V+2 (41) Thus the variance (or uncertainty) of the interpolant can be evaluated at any point by looking up the relevant interval and evaluating a couple of pre calculated quartic polynomials. The result is piecewise 7th order in x with continuous 3rd derivative.

Example data

An example data set is shown in Fig. 3 and will be used to provide a simple comparison between calibration using a conventional least squares linear regression approach and calibration according to a preferred embodiment of the present invention. Fig. 3 shows X (true) data, corresponding Y (measured) data and the standard deviation of the Y (measured) data. The data presented in Fig. 3 is shown plotted in Fig. 4. As can be seen from Fig. 4, data from the middle region of X values has been deliberately left absent from the data set in order to illustrate the differences between using a conventional calibration approach and the method of calibration according to the preferred embodiment. - 19

Using a simple weighted linear regression, a straight line of the form You = a+bX was fitted to the data. The calculated coefficients together with standard deviation were determined to be a = 0.69 _ 0.051 and b = 1.005 + 0.0041.

The 959; confidence interval (predicted mean error on an average of repeated new measurements) may be calculated from: in Sit where n is the number of data points, SO is the standard error of estimate, Sx is the sum of the squared errors in x from the mean value of the x data and t/2 is the probability from the student T distribution. For a 95'-;, confidence interval a/2 = 0.025. This interval allows an estimate of the predicted error in the measured value of Y for a new value of X. This error may be used to predict the error in a calculated value of X from an observed value of Y. Fig. 5 shows the residual error (Ycalc-Y) after linear (i.e. conventional) calibration for each of the X values. The vertical error bars represent the error (+ 2 standard deviations) of the original Y measured data values. The 95% confidence interval for X given a new observation of Y is shown in Fig. 5. As will be seen from Fig. 5, the 95% confidence interval predicts a relatively small uncertainty in the central region of the data where in fact little or no data is actually available.

In contrast, Fig. 6 shows the residual error (Ycalc Y) for each X value after application of the - 20 probabilistic calibration method according to the preferred embodiment. The vertical error bars represent the error (+ 2 standard deviations) of the original Y data values. The predicted uncertainty (+ 2 standard deviations) in X for new observations of Y is shown. As can be seen from Fig. 6 in contrast to the result shown in Fig. 5 the errors predicted are less in the areas where multiple data points are available and are greater where data points are sparse or absent. This reflects the probability that given the standard deviations in the initial Y data, the calibration curve fitted is non- linear. A simple least squares linear regression does not adequately describe the potential errors in this data set.

Although the present invention has been described with reference to preferred embodiments, it will be understood by those skilled in the art that various changes in form and detail may be made without departing from the scope of the invention as set forth in the accompanying claims.

Claims (11)

1. Claims 1. A method of calibrating a mass spectrometer comprising:

   measuring the times of detection of a plurality of ions having known mass to charge ratios; generating a data set comprising said known mass to charge ratios and said times of ion detection) assigning a prior probability distribution function for a flexibility parameter 1, wherein said flexibility parameter is a measure of curvature of said data set; determining the most probable flexibility parameter consistent with said data set and with said prior probability distribution function; and determining the most probable curve consistent with said data set and said most probable flexibility parameter A.
2. 2. A method as claimed in claim 1, wherein said most probable curve is used as a calibration of mass to charge ratio versus time of ion detection.
3. 3. A method as claimed in claim 2, wherein an uncertainty in said calibration of mass to charge ratio versus time of ion detection is determined.
4. 4. A method as claimed in claim 3, wherein said uncertainty in said calibration of mass to charge ratio versus time of ion detection is determined from a measure of a spread of plausible calibration curves about said calibration curve. - 22
5. 5. A method as claimed in any preceding claim, wherein said probability distribution function for said flexibility parameter is proportional to i', wherein n is a constant.
6. 6. A method as claimed in claim 5, wherein 1 < n < 1.5.
7. 7. A method as claimed in claim 6, wherein 1.125 < n < 1.375.
8. 8. A method as claimed in claim 7, wherein n is 1.25.
9. 9. A method as claimed in any preceding claim, wherein said most probable flexibility parameter is that for which a joint probability of said data set and of said probability distribution function for said flexibility parameter is maximized.
10. 10. A method as claimed in any preceding claim, wherein said most probable curve is that for which a joint probability of said data set and of a calibration of mass to charge ratio versus time of ion detection is maximised.
11. 11. A method as claimed in claim 10, wherein said joint probability of said data set and of a calibration of mass to charge ratio versus time of ion detection is that for which the probability of said data set is that given by said most probable flexibility parameter a.