Automating LC-MS/MS mass chromatogram quantification. Wavelet transform based peak detection and automated estimation of peak boundaries and signal-to-noise ratio using signal processing methods
Florian Rupprecht, Sören Enge, Kornelius Schmidt, Wei Gao, Clemens Kirschbaum, Robert Miller
aa r X i v : . [ q - b i o . Q M ] J a n Automating LC–MS/MS mass chromatogram quantification:Wavelet transform based peak detection and automatedestimation of peak boundaries and signal-to-noise ratio usingsignal processing methods.
Florian Rupprecht a , S¨oren Enge a,b , Kornelius Schmidt a , Wei Gao a , ClemensKirschbaum a , Robert Miller c,a a Technische Universit¨at Dresden (TUD), Department of Psychology, Dresden, Germany b Medical School Berlin (MSB), Department of Psychology, Faculty of Science, Berlin, Germany c Health Technology Assessment & Outcomes Research, Pfizer Germany GmbH, Linkstr. 10, 10785Berlin, Germany
AbstractBackground and Objective:
While there are many different methods for peakdetection, no automatic methods for marking peak boundaries to calculate area un-der the curve (AUC) and signal-to-noise ratio (SNR) estimation exist. An algorithmfor the automation of liquid chromatography tandem mass spectrometry (LC–MS/MS)mass chromatogram quantification was developed and validated.
Methods:
Contin-uous wavelet transformation and other digital signal processing methods were used ina multi-step procedure to calculate concentrations of 6 different analytes. To evalu-ate the performance of the algorithm, the results of the manual quantification of 446hair samples with 6 different steroid hormones by two experts were compared to thealgorithm results.
Results:
The proposed approach of automating LC–MS/MS masschromatogram quantification is reliable and valid. The algorithm returns less non-detectables than human raters. Based on signal to noise ratio, human non-detectablescould be correctly classified with a diagnostic performance of AUC = 0 . Con-clusions:
The algorithm presented here allows fast, automated, reliable, and validcomputational peak detection and quantification in LC–MS/MS.
Keywords:
Chromatography, Peak quantification, Peak boundaries, Signalprocessing, Automation, LC–MS/MS
1. Introduction
A very time consuming part of LC–MS/MS quantification and related methods ismanual peak marking in the raw mass chromatograms. While many different methodsfor peak detection have been developed [1], to our knowledge, there is no open sourcemethod which also is able to mark peak boundaries for area under the curve (AUC)measurement and perform signal-to-noise ratio (SNR) estimation. This article describes
Preprint submitted to Biomedical Signal Processing and Control January 25, 2021 n automated procedure capable of performing these tasks, handles fluctuations inretention time (RT) and perform calibration to a known reference standard.There are three direct advantages of analyzing LC–MS/MS mass chromatograms byimplementing an automated process: Increased productivity and speed by decreasingmanual labour, increased objectivity and reliability resulting out of the absence ofhuman raters and, resulting from the non-proprietary nature of the described method,increased replicability [2].This article first describes underlying concepts and models applied by the algorithm.Then a test with real world data is performed and a comparison of the algorithm resultsto results manually quantified by human experts is made.Finally the following two hypotheses were evaluated:1. The algorithm returns analyte concentrations which correlate strongly with thosemanually quantified by experts.2. The algorithm returns analyte concentrations which are valid.The manuscript is structured as follows. In section 2 the algorithm and its theoret-ical background are described, in section 3 algorithm results are compared to resultsmanually compiled by human raters.
2. Theory and calculation
The following sections describe an algorithm for automatic LC–MS/MS mass chro-matogram quantification. The first section introduces a generic model of chromatogramcomposition which forms a theoretical basis for the subsequent sections about decom-position, transformation and processing of measured chromatograms. The final sectionsthen outline a multi-step procedure that is necessary to calculate the absolute concen-trations of an analyte.
A chromatogram can be described as a time series of LC–MS/MS detector intensities y ( t ) ∈ R , where t ∈ R + indicates time. The relevant sequence containing the peak canbe constructed from several component series using the following additive model. y ( t ) = B ( t ) + P ( t ) + N ( t ) , (1)where B ( t ) = b + ( t − t R ) · d is a linear trend component containing background b ∈ R + and linear drift d ∈ R of the chromatogram, P ( t ) is the peak distributioncomponent, and N ( t ) is an irregular component containing noise (besides low frequencybackground noise, this also can contain nearby peaks and low frequency noise or drift).These components are illustrated in Figure 1.Retention time (RT, t R ) is the time at which the peak component has its maximumintensity. Peak height ( h ) is background ( b = B ( t R )) subtracted by signal intensity atRT y ( t R ) − b . Slope ( s ) is the slope of the background component. All these measuresare annotated in Figure 1. 2o measure analyte concentration, the area of the peak needs to be calculated.Peak components are determined by performing several data transformations on thechromatogram time series. Time I n t en s i t y Peak HF noise BG drift Combined ht L t U t R Peak rangeb (cid:1) x (cid:0) ys = (cid:2) x (cid:3) y Figure 1: Generic model for chromatographic peaks. The actual LC–MS/MS sensor data ( y ( t )) iscombined background and drift ( B ( t )), peak distribution ( P ( t )) and noise ( N ( t )). Note that thisFigure does not include peak skew or nearby peaks. Peak range is the time interval between the lowerbound ( t L ) and upper bound ( t U ). Retention time (RT, t R ) indicates maximum sensor intensity inthe peak range. Peak height ( h ) is background ( b = B ( t R )) subtracted by sensor intensity of thecombined signal at RT y ( t R ) − b . Slope ( s ) is the slope of the background component. Note that inthe implementation background and slope are defined using the line between lower and upper boundsof the peak range to simplify the process. This is not necessarily equivalent to the definition in thismodel. The following number of transformations are calculated for a given chromatogram.Figure 2 illustrates the transformations and their order.
Smoothing / low-pass filtering.
The data are smoothed using a 1-D gaussian kernelsmoother [see 3, p. 40-42]. The following formula returns aggregation weights by whichchromatogram intensity values are smoothed. Smoothness can be adjusted using themanually set hyperparameter σ . G ( t ) = 1 √ πσ exp (cid:18) − t σ (cid:19) (2) High-pass filtering.
By differentiation of the original Y ( t ) and smoothed S ( t ) time se-ries, the high-pass filtered data H ( t ) = Y ( t ) − S ( t ) is received. This makes it dependenton the σ parameter of the gaussian kernel.3 igh-pass: H(t) Original: Y(t)
Gaussian smoothing
Smoothed: S(t)
Derivation
Wavelet transformation
CWT
Differentiation
Figure 2: Chromatogram time series transformations. The LC–MS/MS chromatographic time se-ries Y ( t ) is smoothed using a gaussian kernel. High-pass filtered data is obtained by differentiatingsmoothed and original time series. The second derivative of the smoothed data is approximated andthe continuous wavelet transformation (CWT) is calculated. Derivation.
The second derivative D ( t ) = S ′′ ( t ) of the smoothed data S ( t ) is approxi-mated by applying S ′ (cid:18) t i + t i − (cid:19) = S ( t i ) − S ( t i − ) t i − t i − (3)twice. Wavelet decomposition.
Continuous wavelet transformations are used to create time-frequency breakdowns of the chromatograms [see 4].Exponentially modified Gaussian distributions are commonly used to model chro-matographic peaks [5]. The second derivative of the gaussian distribution function isused to build a mother wavelet. Chromatographic peaks have a lower frequency thanthe interfering noise. This enables the separation and further analysis of the noise andpeak distributions by investigating CWT frequency breakdowns.Continuous wavelet transforms are used to divide the data into wavelets. Thecontinuous wavelet transform of a function x ( t ) with scale a ∈ R + and a translationvariable b ∈ R can be expressed by the following integral: X w ( a, b ) = 1 √ a Z ∞−∞ x ( t ) ψ (cid:18) t − ba (cid:19) dt (4) The original mass chromatogram is unevenly spaced. While the transformations up to this pointcan handle this type of data, the continuous wavelet transform implementation is constrained to theuse of evenly sampled data. The data are then resampled in even time steps using linear interpolation.If data loss is a concern the sampling resolution can be increased. mexican hat wavelet [6]: ψ ( t ) = 2 √ σπ / − (cid:18) tσ (cid:19) ! exp (cid:18) − t σ (cid:19) (5)Note that the processed sensor data also has to be differentiated so that the Gaussianderivative wavelet fits to the peak distribution. This furthermore has the advantage,that both background and linear drift are eliminated from the raw mass chromatogram.Then the continuous wavelet transform is used to construct a time-frequency repre-sentation of a signal. Figure 3 shows a chromatographic time series, its second deriva-tive, and the resulting CWT coefficient matrix. − Time I n t en s i t y S c a l e ● ● ● ● ● Coefficient
Figure 3: Chromatographic time series, second derivative (scaled for visibility), and the resultingcoefficient matrix after wavelet transformation. The time-intensity graph also contains smoothed andhigh-pass filtered (HPF) time series, as well as relevant points of the retention time fitness function(RT fit) which is referred to in section 2.3.4.
Each single chromatogram is analyzed in multiple steps. Possible peaks, their RTs,and area boundaries are identified first. Then signal-to-noise ratio (SNR), as wellas peak area, height, background and slope are calculated. Lastly a fitness value iscalculated by which the peak most likely associated with the analyte is selected.5 .3.1. Peak identification & boundary estimation
Peak RT ( t R ) and boundaries ( t L , t U ) (see Figure 1) are estimated in multiple steps.Figure 4 serves as a visual aid for picturing changes in the peak polygon in each step. Time I n t en s i t y Smoothed Original t ^ U (a) Peak area after initial boundaryestimation. Time I n t en s i t y Smoothed Original t ^ U^ (b) Peak area after friction bound-ary correction. Time I n t en s i t y Smoothed Original t ^ t U^ (c) Peak area after partial convexhull boundary correction. Figure 4: This Figure illustrates changes in the area polygon defined by estimations for lower ( b t L ) andupper ( c t U ) peak boundaries after initial approximation (4a) and the two subsequent correction steps.First a Newtonian friction inspired approach is used in conjunction with the smoothed time series towiden the bounds (4b). The resulting potentially self-overlapping area polygon is then modified byfirst calculating the lower part of the convex hull and removing separate polygons after (4c). Peak retention time (RT).
Non-maximum suppression (NMS) [7] is used to identifylocal maxima in the CWT coefficient matrix. NMS chooses points that have exclusivelygreater value neighbouring points. Local maxima in the CWT matrix correspond to lowfrequency components, peak components in this case, in the chromatogram data. Forthis reason NMS returns a list of peak RT with associated time, scale and coefficientvalues, each representing a potential peak in the data.
Peak boundary estimation.
To determine the peak area, lower and upper bounds of thepeak component are needed. These are approximated in two steps.While peak RT ( t R ) is projected to a local maxima in the CWT coefficient matrix,approximate inflection points of the peak distribution should correspond to minima.The vector of CWT coefficients with the same scale as the previously determined maxi-mum is inspected for the closest minima next to the maximum. This gives us estimatesfor lower and upper peak bounds. Friction boundary correction.
As the wavelet transformation models uniformly dis-tributed gaussian distributions and chromatographic peaks are often subject to skew,these approximate boundaries need to be corrected.A Newtonian friction inspired approach is performed. Similar to a physical bodygradually sliding down a hillside until a slope is reached at which friction stops it frommoving, the initially estimated boundaries of the peak component are increased untila specified slope threshold is reached. This is detailed in Algorithm 1. It is advised tonormalize the threshold to the intensity range of each analyzed chromatogram.6 lgorithm 1
Peak boundary correction. The size function returns the size of a list.
Input: Y : Smoothed chromatogram data Input: p : List of peak boundaries Input: t : Threshold t ← ( max ( Y ) − min ( Y )) · t for each ( u, l ) ∈ p do while u < size ( Y ) − do if Y u − Y u +1 > t then u ← u + 1 while l > do if Y l − Y l − > t then l ← l − Partial convex hull boundary correction.
At this point resulting peak polygons mayoverlap themselves at their lower and upper boundaries when there is significant noisein the data. To circumvent this and to mimic the behavior of a human rater the bottom of the convex hull of the peak data is calculated. A variation of Grahams scan [see 8]is implemented. It is described in Algorithm 2. Next, all areas of the polygon separatedfrom the main peak area are removed. Peak boundaries are then set to the minimumand maximum x (time) value of the polygon vertices and peak RT to the maximum y (intensity) value. Algorithm 2
Partial convex hull boundary correction. Modification of Grahams scan[see 8]. The size function return the current size of the stack, push and pop are standardstack operations of adding and removing an element and orient returns the orientation o (clockwise ( o > o <
0) or colinear ( o = 0)) of an ordered triplet. Input: W : Empty stack Input: p : Vector of points ordered by increasing x coordinate push ( W, p ) push ( W, p ) for i ← N do while size ( W ) ≥ orient ( p i , W top , W second ) ≤ do pop ( W ) push ( p ) Peak area ( a ∈ R + ), height ( h ∈ R + ), background ( b ∈ R + ) and slope ( s ∈ R )can be calculated after the peak boundaries are known (see Figure 1). Note that of When the points of the time series between temporal lower and upper boundaries are consideredpoints on a time-intensity plane, the convex hull is the smallest convex polygon that contains them all.
Area under the curve (AUC).
The AUC is defined as the area integral chromatographicdata constraint to the peak component subtracted by the area integral of the slopebetween lower and upper peak bounds. This is equivalent to the peak polygon area andcan easily be calculated using gauss’s area formula: a = 12 (cid:12)(cid:12)(cid:12) n − X i =1 x i y i +1 + x n y − n − X i =1 x i +1 y i − x y n (cid:12)(cid:12)(cid:12) , (6)where x ∈ R + are the temporal and y ∈ R + the intensity dimension coordinates ofall measurements between peak bounds. Peak slope.
The slope s ∈ R is defined through the time positions of lower and upperpeak boundaries: s = y ( t U ) − y ( t L ) t U − t L , (7)where y ( t ) ∈ R + is the chromatogram time series and t U ∈ R + is the time of thelower peak boundary. Peak background.
The background b ∈ R + is determined by the temporal position ofthe peak maximum on the peak slope with the lower peak bound as the line intercept: b = y ( t L ) + ( t R − t L ) · s , (8)where t L ∈ R + is the time of the lower peak boundary. Peak height.
The height h ∈ R + is defined as the difference of the peak maximum toits background: h = y ( t R ) − b , (9)where t R ∈ R + the time of the peak maximum. While there is a lack of definitive guidelines for SNR measurements, signal is oftendefined as maximum peak height above baseline and noise as standard deviation orroot mean squared of a manually selected part of the chromatogram where no peaksare present [9].We propose to approximate SNR ( r ∈ R + ), by dividing standard deviation of thehigh-pass data ( H ( x ) ∈ R + , mean H , number of data points N ) by doubled peak height( h ): 8 = q N − P Ni =1 ( H ( x ) − H ) h (10)This method can be performed with very little or no data points outside of thepeak bounds. Note that this method of estimating SNR is dependent on the smoothinghyperparameter σ . This means that SNR can only be compared across analytes whenthe same σ is used. The peak identification procedure performed previously returns multiple possiblepeaks.To select the appropriate one a peak fitness function f ( t R ) = c · g ( t R ) (11)is calculated for each peak depending on its RT ( t R ∈ R + ), its CWT coefficient( c ∈ R + ) and a RT fitness function ( g ( t R )). The peak with the highest fitness value isselected and peaks with negative signed fitness are rejected completely. Note that thiscan lead to cases where no peak is detected in a chromatogram. These are considerednon-detectables (ND).The RT fitness function g ( t R ) is used for the selection of peaks depending on availableinformation about the expected peak RT. When no information about the expected RTis known, the RT fitness function is equivalent to g ( t R ) = 1. This leads to selectionof the peak with the highest CWT coefficient and rejects none. In other cases wherean estimate of the RT ( b t R ∈ R + ) and a range hyperparameter ( a ∈ R + ) is available, aquadratic fitness function g ( t R ) = 1 − t R − b t R a ! , (12)is created to weigh peaks by proximity to expected RT and leads to rejected peaksoutside of the range. A whole dataset is comprised out of multiple samples. Each sample contains dataassociated with multiple measured analytes. Besides the chromatogram by which an(unknown) analyte concentration is measured directly, it also contains a chromatogramof an internal standard analyte (IS) with a known analyte concentration. The use ofstable isotopically labelled internal standards is a commonly used technique that helpsto control variability in quantitative assays [10, 11]. There are also several calibrationsamples which are used for concentration calibration. These samples contain a knownstaggered amount of the analyte in addition to the known amount of IS. The single chro-matogram processing described in the previous section, is modified by hyperparametersspecified for each analyte and IS. 9he overall process of whole dataset quantification consists out of three main parts,RT calibration, concentration calibration, and final quantification. These are describedin the following sections. Figure 5 contains a system diagram which gives an overviewover all parts of the routine.
Analyte + IS data
Samples with knownanalyte concentrationRelative known concentration ProcessingparametersIS ProcessingparametersAnalyteAll process
IS All process a nalyte
Sample 1
IS AnalyteKnown concentration Known concentration?
Sample 2
IS AnalyteKnown concentration Known concentration?
Sample 3
IS AnalyteKnown concentration Known concentration?
Sample n
IS AnalyteKnown concentration Known concentration? SNRRTAreaHeight, BG, etc. process
Processing parameter RT?IS s u b t r a c t
1: RT calibration
Processing parameter RT? Process performedfor each single sampleOutputISAnalyte ChromatogramSample dataRaw data
Figure 5: System diagram for a single analyte and its internal standard of a complete dataset. Samplescontain an IS chromatogram with a known concentration and an analyte chromatogram which mayhave a known concentration. There are processing parameters for both analyte and IS which modifysingle chromatogram processing. RT calibration, concentration calibration and the final quantificationare three separate steps which each process a different set of samples. Green boxes indicate outputsof processes.
While the RT of an analyte should be known approximately, it has some variance[12, 13]. The aim of this step is to analyze analyte chromatograms in which the rightpeaks are easily identifiable, to provide better estimates for the expected RT of a specificanalyte and subsequently make peak selection easier in all other samples of the analyte.Out of the calibration samples all are selected which contain a sufficient quantity ofthe calibrator analyte to be able to easily and distinctly detect the mass chromatogram10eaks. Selected samples have a known concentration ratio equal or greater than athreshold hyperparameter x ( C A /C IS ≥ x , where C A ∈ R + is the known concentrationof the analyte and C IS ∈ R + is the known concentration of the IS).Both IS and analyte chromatogram are processed as described in previous sections.When wrong peaks are present in the data a guessed rough estimate of the expectedRT and an acceptable range in the analyte parameters can be specified to create apreliminary RT fitness function (see 2.3.4).After obtaining all peak time positions of the selected samples, their arithmeticmean is used as an expected RT difference.Now a calibrated RT for both analyte and internal standard has been calculated.The data showed that the standard deviation of both of them individually is consis-tently greater than the standard deviation of their difference (see results Table 1). Insubsequent analyses the IS peak RT ( t ISR ∈ R + ) is retrieved first. IS peaks are generallyeasier to identify because of their constantly high concentration. Then the previouslycalculated mean RT difference ( t ∆ ∈ R ) is used to calculate the expected RT of theanalyte t AR = t ISR + t ∆ .In combination with a smaller acceptable RT range, this allows to build more accu-rate, dynamic RT fitness functions for the next chromatogram processing steps.Analyte t AR t ISR t ∆ β Cortisol 4.687 (0.018) 4.682 (0.019) 0.005 (0.007) 2.41Cortisol Q 4.688 (0.018) 4.682 (0.019) 0.006 (0.008) 4.35Cortisone 4.506 (0.014) 4.493 (0.013) 0.013 (0.007) 1.61Testosterone 5.841 (0.044) 5.810 (0.040) 0.036 (0.007) 0.31Progesterone 7.108 (0.055) 7.042 (0.054) 0.066 (0.006) 1.07DHEA 6.135 (0.046) 6.098 (0.045) 0.036 (0.007) 0.76
Table 1: Calibration data for RT and concentration calibration. Measured RT mean and standarddeviation (in parentheses) of the RT calibration samples (see section 2.4.1). Note that the standarddeviation of analyte RT ( t AR ) and IS RT ( t ISR ) is consistently greater then the standard deviation oftheir difference t ∆ . β is the slope of the concentration calibration regression model. Cortisol Q refersto the second most sensitive transition of cortisol, DHEA to dehydroepiandrosterone. The calibration samples, as previously described, contain known concentrations ofthe analyte as well as IS. While only calibration samples with a large known concen-tration ratio were used when calibrating the RT, now all of them are processed.First, the IS chromatogram is processed like before. Then the analyte chromatogramis processed using the dynamic RT fitness function which depends on the IS RT andthe mean RT difference calculated in the previous step. Then the analyte peak area( a A ∈ R +) is divided by IS peak area ( a IS ∈ R +) to obtain a relative measuredconcentration ( M = a A /a IS ). In combination with the relative known concentration11 C = C A /C IS ) available for each sample in this step, their linear relationship is modeledusing linear regression without an intercept. C = βM (13)The resulting regression model can be used to estimate absolute concentrations frommeasured relative concentrations whenever the linearity assumption holds. See Figure 6for a visual example of such a model. ●●●●●●●●●●● ●● ●● ●● ● ●●●●●●●●● ●●●● ●●●●● ●●●●●●●●● ●●●●●●●●● ●●●●●●●●● ●●●●●●●●● ●●●●●●●●● − − − − − − Measured K no w n Figure 6: Plot of the linear regression model for concentration calibration of cortisol using our data.
Absolute sample quantification can now be performed for all samples. The chro-matogram processing is performed just as in the previous, concentration calibrationstep. To obtain absolute measured concentrations the concentration calibration regres-sion function of the previous step is applied to each relative measured concentration.
We supply a open source reference application which implements all parts of thealgorithm. The software reads .mzML data, an open source standard for mass spec-trometry data [14].
3. Method
Specimens Collection and Preprocessing.
The raw mass chromatogram data were madeavailable from the Dresden longitudinal study of chronic stress and cognitive control(StressCog) [see 15]. Study approval was granted by the Dresden ethics committee(dossier EK23012016, IRB00001473, and IORG0001076). In the recruitment process8,400 eligible participants from the population registry of the City of Dresden between25 to 55 years old and within the metropolitan area, were selected and contacted byinvitation letter. The study Protocol encompassed the collection of hair samples with Code will be released upon acceptance of this paper.
12 min lenght of 3 cm at the posterior vertex region of the head which was required foranalyte extraction [see 16]. See Schmidt et al. [17] for more details of the sampling anddata collection process.The analyses in this paper were performed on a subset of 230 hair samples from 229different human subjects (137 females; age range 25 - 55; mean age = 38.05; age SD =8.51) of the StressCog study. Samples were processed following the protocol described inGao et al. [16] with slight adjustments. The non-pulverized samples were washed withisopropanol and steroid hormones were extracted from 7.5 mg hair by methanol incuba-tion. Column switching on-line solid phase extraction (SPE) was performed, followedby analyte detection on an AB Sciex API 5000 QTrap mass spectrometer. When enoughmaterial for two LC-MS/MS measurements was available the samples were split. Thisresulted in 446 samples which were processed by the LC-MS/MS (calibration samplesexcluded) to obtain selected-reaction monitoring (SRM) chromatograms.
Analyte Quantitation.
The performance of the algorithm was evaluated by comparingits results to the manual quantification by two trained experts, referred to as expert1 and expert 2. They are both trained in the manual marking of chromatographicpeaks. Expert 1 was more experienced, and has manually marked approximately 60,000samples with each 6 analytes before analyzing the data used in this paper, while expert 2has marked about a tenth of this sample count, 5,000 samples, with 5 different analytes.225 SRM chromatograms were analyzed by one expert, 221 by two experts. 216 ofthese LC-MS/MS samples were split from the same hair samples, which were thereforemeasured and analyzed twice. Multiplied by 6 different investigated analytes this resultsin 2676 individual measurements analyzed by one or both experts. All samples were alsoanalyzed using the algorithm. These numbers are further visualized in supplementaryFigure A.9.The measured analytes were cortisol, cortisone, testosterone, progesterone and dehy-droepiandrosterone (DHEA). The selection of analytes and the procedure follows [16].The second most sensitive transition of cortisol (indicated by the Q suffix) was alsoincluded in the selection of analytes.To be able to programmatically access the raw mass chromatograms, the data wasconverted from a proprietary file format to .mzML, using the ProteoWizard MSConvertsoftware [18, 19]. Then the supplied software implementation (see 2.5) was used toprocess the data.For the manual quantification data, the AUC of both analyte and IS were markedand extracted by the experts using the AB SCIEX Analyst software version 1.6.2.
Algorithm settings.
Algorithm hyperparameters were mostly the same across analyteswith the exception of fixed RT windows (expected IS RT and acceptable range) whichhad to be specified for the IS of progesterone, testosterone and DHEA, because of wrongIS peaks close to the correct one.The absolute measured concentration (analyte AUC divided by IS AUC) was usedfor all evaluative statistics. Manual compiled relative concentrations were also scaled13y the β of the algorithms concentration calibration regression. This increases compa-rability between analytes for calculations which include all analytes in comparison tousing the relative concentrations. And by not using the manually compiled values tocreate a separate calibration regression systematic errors are not overestimated. Statistical analyses.
Statistical analyses were performed with R version 3.6 [20], latentvariable models were estimated with lavaan [21].
4. Results
Table 2 contains descriptive statistics for measured relative concentrations as wellas SNR of each analyte.
Analyte N Concentration (pg/mg) SNR (Algorithm)Algorithm Expert 1 Analyte ISCortisol 389 8.07 (24.56) 8.10 (25.53) 15.12 (5.57) 25.23 (1.74)Cortisol Q 390 8.01 (23.73) 8.34 (24.62) 11.61 (5.06) 25.32 (1.64)Cortisone 446 16.11 (16.49) 15.91 (16.33) 20.77 (4.80) 24.19 (2.55)Testosterone 178 4.09 (26.43) 4.00 (25.46) 5.35 (6.51) 32.46 (8.29)Progesterone 415 8.09 (64.21) 7.73 (61.16) 13.39 (8.65) 30.81 (4.34)DHEA 414 1.89 (2.60) 1.85 (2.59) 13.74 (6.52) 34.52 (3.34)
Table 2: Descriptive statistics by analyte. Means followed by standard deviations in parentheses. Allcases in which both the expert and the algorithm could detect a peak are included.
Non-detectables (ND).
Chromatograms where an expert cannot find a peak are markedas non-detectable (ND). In these cases, the analyte concentration is presumed to be lessthan the measuring instruments lower limit of detection. These cases were comparedto cases when the algorithm cannot mark a peak with AUC > Y A+S = β + β X + β X , where Y is the probability of being an FP and X and X were AUC and SNR respectively, as well as models which just had AUC and SNRas predictors. As the number of FPs was greatly lower than the number of TPs, thesame number of TPs was randomly sampled to obtain equal group sizes. All threemodels were calculated for all analytes with at least 10 FPs as well as for the combineddata. Overall accuracy (AUROC) for peak AUC was 0 .
91, for SNR it was 0 .
95 and14 nalyte Expert 1 ( n = 446) Accuracy (FP class.)TP FN FP TN Sensitivity Specificity A+S AUC SNRCortisol 389 2 16 39 99.49 % 70.91 % 1.00 0.98 0.87Cortisol Q 390 0 45 11 100.00 % 19.64 % 0.91 0.92 0.75Cortisone 446 0 0 0 100.00 % - - - -Testosterone 178 7 121 140 96.22 % 53.64 % 0.95 0.88 0.95Progesterone 415 0 30 1 100.00 % 3.23 % 0.88 0.88 0.81DHEA 414 0 23 9 100.00 % 28.12 % 0.95 0.94 0.89All 2232 9 235 200 99.60 % 45.98 % 0.95 0.95 0.91Analyte Expert 2 ( n = 221) Accuracy (FP class.)TP FN FP TN Sensitivity Specificity A+S AUC SNRCortisol 199 4 4 14 98.03 % 77.78 % - - -Cortisol Q 205 2 10 4 99.03 % 28.57 % 1.00 1.00 0.93Cortisone 220 0 1 0 100.00 % 0.00 % - - -Testosterone 106 5 45 65 95.50 % 59.09 % 0.98 0.86 0.97Progesterone 217 0 4 0 100.00 % 0.00 % - - -DHEA 212 2 4 3 99.07 % 42.86 % - - -All 1159 13 68 86 98.89 % 55.84 % 0.96 0.96 0.94 Table 3: Contingency data for peak detection. The algorithms sensitivity of peak detection wasvery high. Specificity was very low, mainly due to false positives (FP). Binomial regression wasable to classify FPs (abbreviated FP class.) with high accuracy (area under the receiver operatingcharacteristic). Equal group sizes were randomly sampled for hormones with at least 10 false positives.A+S indicates the regression model with both AUC and SNR as predictors. the combined model had an accuracy of 0 .
95. Accuracies for all calculated models aredisplayed in Table 3. Figure 7a is a scatterplot of the resampled cases with peak AUCand SNR on the axes. Figure 7b shows the receiver operating characteristic (ROC) forboth dimensions.
Correlation & differences.
Pearson correlation coefficients of the log transformed con-centrations of all samples and analytes were all r = 0 .
98 between algorithm and expert1 and 2 as well as between the experts. Table 4 shows single correlation coefficients forall analytes and raters. The supplementary Figure A.10 shows scatter plots for all an-alytes and the supplementary Figure A.11 shows Bland-Altman plots of the differencesbetween raters.
Latent variable model.
The data contained several LC–MS/MS measurements of thesame original samples, rated by the two experts and the algorithm (see supplementaryFigure A.9). An analogous structural equation model was created. Expert and Algo-rithm ratings were considered manifest variables, LC–MS/MS data was considered firstorder latent variables and true hormone concentration as second order latent variable.15 − − A(cid:4)(cid:5)(cid:6)(cid:7) ithm: Peak AUC A l go r i t h m : S NR True positive False positive (a) Scatterplot of expert 1 true and false posi-tives.
S(cid:8)(cid:9)(cid:10)(cid:11)(cid:12)(cid:13)(cid:14)(cid:15)(cid:16)(cid:17)(cid:18)(cid:19)(cid:20)(cid:21)(cid:22)(cid:23)(cid:24)(cid:25)(cid:26)(cid:27)(cid:28) A U(cid:29) (cid:30) (cid:31) !
AUC SNR (b) ROC for classifying expert 1 FPs. Accuracy(AUROC) of AUC: 0.95, SNR: 0.89, combinedmodel: 0.95. The AUC model graph and thethe combined model graph are oversecting.
Figure 7: Binominal regression classifiers for false positive peak detection. AUC and SNR as measuredby the algorithm were investigated for their ability to predict false positives. These Figures show dataof all analytes for expert 1. True positives were randomly sampled to obtain equal group sizes.
Analyte Algorithm Algorithm Expert 1Expert 1 Expert 2 Expert 2Cortisol 0.981 0.981 0.985Cortisol Q 0.950 0.972 0.952Cortisone 0.998 0.983 0.988Testosterone 0.936 0.951 0.945Progesterone 0.970 0.964 0.970DHEA 0.953 0.898 0.952All 0.983 0.979 0.983
Table 4: Pearson correlation coefficients for log transformed hormone concentrations.
Based on the latent state-trait theory [22], the total variance of analyte concentra-tion per specimen replicate was decomposed into estimates of (1) reliability, (2) replicatespecificity and (3) replicate consistency for each quantitation method (i.e., algorithm,rater 1, and rater 2). The model’s structure and its parameters are shown in Figure 8.Additional method-specific variance components [see 23] were neither numerically iden-tifiable nor necessary to account for the covariance among cortisol cortisone, and DHEAquantitations. The fit between the observed and the model-implied covariance struc-ture was close for these analytes ( χ ≥ . p ≥ . ≤ . ≤ . χ ≥ . p ≥ . ≤ . ≤ . A E E A Sample 1 Sample 2Hormone ψ ψβ βλ λ λ λ λ ψ = 1 Figure 8: Path diagram of the measurement model. Rectangles represent manifest variables (Analyteconcentrations estimated by experts E and E as well as by the algorithm A ) for each LC–MS/MSsample. Rounded rectangles represent latent factors. First order latent factors are the measuredsamples, the second order latent factor is the hormone concentration. Residual variances are notshown in the diagram. the quantitation of progesterone by the algorithm differed slightly from the humanraters, whereas there were differences in the within-replicate variance of testosteroneirrespective of the quantatitation method. Reliability, specificity, and consistency ofeach quantitation method and analyte are reported in Table 5.Consistency Specificity Reliability A E E A E E A E E Cortisol 0.801 0.787 0.732 0.185 0.182 0.169 0.986 0.969 0.901Cortisol Q 0.760 0.762 0.708 0.206 0.207 0.192 0.966 0.968 0.900Cortisone 0.774 0.757 0.633 0.220 0.215 0.180 0.994 0.971 0.813Testosterone 0.869 0.896 0.699 0.102 0.105 0.082 0.971 1.002 0.782Progesterone 0.488 0.495 0.439 0.498 0.506 0.448 0.986 1.001 0.888DHEA 0.819 0.839 0.663 0.140 0.144 0.114 0.959 0.982 0.776
Table 5: Consistency, specificity, and reliability of analyte quantitation (2nd replicate) by the algorithm( A ), and both expert raters ( H , H ).
5. Discussion & Limitations
Measurements automatically calculated by the algorithm correlate very stronglywith those manually compiled by experts. The algorithm returns far less non-detectablesthan expert raters. It was demonstrated that these peak detection false positives canbe classified by the (algorithm measured) peak AUC and SNR. This leads us to rec-ommend setting a minimum threshold for one or both of these values for each analytewith high ND rates (e.g. testosterone).Latent variable models showed that the algorithm is able to measure the chro-matograms trait very similar to the experts, while being closer to the more experiencedexperts results. Consistency might be slightly overestimated due to list-wise exclusion.17ubsequent studies should experiment with modifications of algorithm parametersand the effect of disabling for example, the partial convex hull peak boundary correctionwhich was designed to imitate expert rating behavior or use the smoothed peak areafor calculating AUC. These modifications could potentially decrease model error.Furthermore, it should be investigated if the algorithm is reliable with other types ofdata. These could be other LC–MS/MS mass chromatogram types as well as other typesof chromatograms and spectrograms. If nonlinear drift would emerge as a problem whenanalyzing different types of data, we suggest investigating the addition of a high-passfilter.
In conclusion, the algorithm presented here allows fast, automated, reliable andvalid computational peak detection and quantification in LC–MS/MS and is also ableto quantify SNR automatically.
6. Acknowledgements
Special thanks to the employees of Dresden LabService GmbH and Christian Rup-precht (Department of Engineering Science, University of Oxford) for their continuoussupport and suggestions. This work was supported by the German Research Foundation(DFG, Grant No. SFB 940/2).
7. Conflict of interest
The authors have no conflict of interest to declare.
ReferencesReferences [1] C. Yang, Z. He, W. Yu, Comparison of public peak detection algorithms forMALDI mass spectrometry data analysis, BMC Bioinformatics 10 (1) (2009) 4,doi:10.1186/1471-2105-10-4.[2] H. E. Plesser, Reproducibility vs. Replicability: A Brief History of aConfused Terminology, Frontiers in Neuroinformatics 11 (2018) 76, doi:10.3389/fninf.2017.00076.[3] E. R. Davies, Computer and machine vision: theory, algorithms, practicalities,Academic Press, 2012.[4] I. Daubechies, Ten Lectures on Wavelets, vol. 61, Society for Industrial and AppliedMathematics, doi:10.1137/1.9781611970104, 1992.185] E. Grushka, Characterization of exponentially modified Gaussian peaksin chromatography, Analytical Chemistry 44 (11) (1972) 1733–1738, doi:10.1021/ac60319a011.[6] C. Torrence, G. P. Compo, A practical guide to wavelet analysis, Bulletin of theAmerican Meteorological society 79 (1) (1998) 61–78.[7] A. Neubeck, L. V. Gool, Efficient Non-Maximum Suppression, in: 18th Interna-tional Conference on Pattern Recognition (ICPR ' , 2013.[12] A. A. Klammer, X. Yi, M. J. MacCoss, W. S. Noble, Improving Tandem MassSpectrum Identification Using Peptide Retention Time Prediction across DiverseChromatography Conditions, Analytical Chemistry 79 (16) (2007) 6111–6118, doi:10.1021/ac070262k.[13] D. Guo, C. T. Mant, R. S. Hodges, Effects of ion-pairing reagents on the predictionof peptide retention in reversed-phase high-resolution liquid chromatography, Jour-nal of Chromatography A 386 (1987) 205–222, doi:10.1016/s0021-9673(01)94598-4.[14] L. Martens, M. Chambers, M. Sturm, D. Kessner, F. Levander, J. Shofstahl, W. H.Tang, A. R¨ompp, S. Neumann, A. D. Pizarro, L. Montecchi-Palazzi, N. Tasman,M. Coleman, F. Reisinger, P. Souda, H. Hermjakob, P.-A. Binz, E. W. Deutsch,mzML—a Community Standard for Mass Spectrometry Data, Molecular & Cellu-lar Proteomics 10 (1), doi:10.1074/mcp.r110.000133.[15] R. Miller, K. Schmidt, C. Kirschbaum, S. Enge, Chronic Stress and ExecutiveFunctioning: A longitudinal perspective, doi:10.17605/OSF.IO/KDMX5, URL osf.io/kdmx5 , 2019.[16] W. Gao, T. Stalder, P. Foley, M. Rauh, H. Deng, C. Kirschbaum, Quanti-tative analysis of steroid hormones in human hair using a column-switching19C–APCI–MS/MS assay, Journal of Chromatography B 928 (2013) 1–8, doi:10.1016/j.jchromb.2013.03.008.[17] K. Schmidt, S. Enge, R. Miller, Reconsidering the construct validity of self-reportedchronic stress: A multidimensional item response theory approach., PsychologicalAssessment 32 (11) (2020) 997–1014, doi:10.1037/pas0000829.[18] D. Kessner, M. Chambers, R. Burke, D. Agus, P. Mallick, ProteoWizard: opensource software for rapid proteomics tools development, Bioinformatics 24 (21)(2008) 2534–2536, doi:10.1093/bioinformatics/btn323.[19] M. C. Chambers, B. Maclean, R. Burke, D. Amodei, D. L. Ruderman, S. Neu-mann, L. Gatto, B. Fischer, B. Pratt, J. Egertson, K. Hoff, D. Kessner, N. Tas-man, N. Shulman, B. Frewen, T. A. Baker, M.-Y. Brusniak, C. Paulse, D. Creasy,L. Flashner, K. Kani, C. Moulding, S. L. Seymour, L. M. Nuwaysir, B. Lefebvre,F. Kuhlmann, J. Roark, P. Rainer, S. Detlev, T. Hemenway, A. Huhmer, J. Lan-gridge, B. Connolly, T. Chadick, K. Holly, J. Eckels, E. W. Deutsch, R. L. Moritz,J. E. Katz, D. B. Agus, M. MacCoss, D. L. Tabb, P. Mallick, A cross-platformtoolkit for mass spectrometry and proteomics, Nature Biotechnology 30 (10) (2012)918–920, doi:10.1038/nbt.2377.[20] R Core Team, R: A Language and Environment for Statistical Com-puting, R Foundation for Statistical Computing, Vienna, Austria, URL , 2020.[21] Y. Rosseel, lavaan: An R Package for Structural Equation Mod-eling, Journal of Statistical Software 48 (2) (2012) 1–36, URL .[22] R. Steyer, M. Schmitt, M. Eid, Latent state–trait theory and research in personalityand individual differences, European Journal of Personality 13 (5) (1999) 389–408,doi:10.1002/(SICI)1099-0984(199909/10)13:5¡389::AID-PER361¿3.0.CO;2-A.[23] C. Geiser, G. Lockhart, A comparison of four approaches to account for methodeffects in latent state–trait analyses., Psychological methods 17 (2) (2012) 255,doi:10.1037/a0026977. 20 ppendix A. Supplement Abbreviation ExplanationAUC Area under the curveAUROC Area under the ROC = AccuracyIS Internal standardLC–MS/MS Liquid chromatography tandem mass spectrometryROC Receiver operating characteristicSNR Signal to noise ratio
Table A.6: Abbreviations
Hair sample L C -MS/MS sampleLC-MS/MS sample Chromatogram pairsfor all 6 analytesChromatogram pairsfor all 6 analytes Expert 1AlgorithmExpert 1Expert 2Algorithm230 225221 13501326 Σ ∩ Σ ∩ Figure A.9: Visualization of sample number and dependence. 230 different hair samples were split andprepared to 446 LC-MS/MS samples. This Figure displays the resulting number of chromatogramsand how many of them were repeated measurements. Σ indicates the total sum including split samplesand ∩ the number of intersecting original samples. ● ● ●● ●● ●●●●● ●● ●●●●● ● ●● ●●●●●●●● ●● ●●●●● ● ●●●● ●●● ●●●●●●● ●●● ●●●● ●●●●● ●●●●●● ●● ●● ●●● ●●● ●●● ●●● ●● ●● ● ●● ● ●●●● ● ●●●●● ●● ●●●● ● ● ●●●●●●● ●●● ●●● ●●●● ●●●● ●●● ●● ●●●●● ●●● ● ●●●● ●●●●● ● ●● ●● ● ●● ●●● ●●●●● ● ● ●●●● ●●● ●●● ●●●● ●●●● ●●●●● ● ●●●● ●●●● ●●●● ●●●●● ●●●●●● ● ●● ●●●● ●● ● ●●● ●● ●●●●● ● ●● ●●● ● ●● ●● ●● ●●●● ● ● ●●● ●● ●●●●●●●●●●●●●● ●● ● ●●●●●● ●●● ●●●● ●●● ●● ●●●●● ●●● ●●● ●●●●● ●●● ● ● ●●●● ●●●● ●●●●●● ● ●●● ●●●● ●●●● ●● ●●●● ●●● ●●●● ●●●●● ●●●●●● ● ●●●●●● ●●● ●● ●● ●● ●● − − Algor i" (a) Cortisol. ● ● ● ●●● ●●●● ● ●●●● ●● ●●●●●● ● ● ●● ●●●●● ● ●●●● ●● ●●●● ●● ●●●●● ●● ● ●●●● ● ●●● ●●● ●●● ●● ● ●● ●●●● ●●● ●●● ●●● ●●● ●● ● ●● ●●● ●●● ●● ●● ●● ●●●●●● ● ●●●●● ● ● ●●●● ●● ● ●● ●●●●● ● ●●●●● ●●●● ● ●●●●●● ● ● ●● ●● ●● ●●● ●●●●● ●● ●●● ●●● ●●● ●●●● ●●● ● ●●● ●●●●● ●●●● ●● ● ●●● ●●● ●●●● ●●● ●●●●● ●●●● ●● ● ●● ●●● ● ●● ●●●●● ● ●●●●●● ●● ●● ●● ●●●● ●●● ●● ●●●● ● ● ●●● ● ● ●●●●●●●●●●●●● ●●● ●●●●●● ●●● ● ●●● ●●● ●● ●●● ●●● ●●● ●●● ●●● ●● ●●● ●●●● ●●● ●● ●● ●●●● ●●● ●●● ●●●● ●●● ●●● ●●● ●● ●●●●● ● ● ●●● ● ●●●● ●●●● ●●● ●● ● ●● ●● − − Algor (b) Cortisol Q. ● ● ●● ●●●●●●● ●●●●●●● ●●●●● ●●● ●● ●●●●●●●● ●●● ● ●● ●● ●● ●●● ● ● ● ● ●●●● ●●● ● ● ●●●●● ●● ●● ● ● ● ●●●●●●● ●●●●● ●● ●● ●●●● ● ●●●●● ●●● ●●● ●●● ●●● ●●●●● ●● ●●●●●●●●●● ●●● ●●● ●●● ●●●●●●● ● ●●●●● ●●● ● ● ●●● ●●● ●● ●●●● ●●●● ● ●●● ● ●●●●●● ●● ●● ●● ●● ● ●●●●●●●● ●● ● ● ●●● ●●●●●● ●●●● ●●●● ●●●●● ● ●●●●●● ●● ●●●●●●● ●●●●● ● ● ●●●●● ●●●●●●●●● ●●●● ●● ●●●● ●● ●●●●●● ●● ●●●●●●●●●● ●●●● ●●●●●● ●● ●●●● ●● ●●●●●● ●●●●●●●● ●● ● ●●●● ●●●● ●●● ● ●●●●● ● ●● ●●● ●● ●●● ●●● ●●● ●●●●●●●●●●● ●●●● ●●●●● ●●●●● ●●●● ●●●● ●● ●●●● ●●● ● ● ●●● ●●●● ●●●●● ●●● ●● ●●●● ●●● ● ●●●●● ●●●●●●● ● ●● − − Algor
TVWX YZ[\]^_‘abcdefghjklmn (c) Cortisone. ●●● ●● ● ●● ●● ●●●● ●● ●●●● ●● ●●● ●●● ●● ● ●●●●●● ● ●●● ●●●● ●●●● ●●● ●●●●● ●● ● ●● ●● ●●●● ● ●● ●●●●●● ● ●● ●●● ●● ● ● ●●● ●● ●● ●●● ●●●●● ●● ● ● ●●● ● ●● ●●●● ● ●●●● ● ●●●●● ●●●● ●●● ●● ●●● ●●● ●● ●● ● ● ●●●● ●●●● ●●● ●●●●●●● ● ●● ● ● ●● ●●● ● ● − − − Algor opqr suvwxyz{|}~(cid:127)(cid:128)(cid:129)(cid:130)(cid:131)(cid:132)(cid:133)(cid:134)(cid:135)(cid:136) (d) Testosterone. ●● ●●●●●● ●● ●●● ●●● ●● ●●● ●●●●●● ●●● ●● ●● ● ● ●● ●● ●●●●●●●●●● ●●●● ●● ●●●●●●●● ●● ●● ● ●● ●●●● ● ● ●● ●● ●● ●●● ●●●● ●● ●●● ●●●●● ●● ●●● ●● ●●● ●●● ● ●● ● ●●● ● ●●●● ●●●●● ●● ●●● ●● ●●●● ● ●●● ● ● ●●●●●● ●● ●● ●●● ●●●●●●● ●● ●●● ●●● ●● ●●●●● ●●● ●●●● ●● ● ●● ●● ●●●●●●● ●● ●● ● ● ● ● ●●● ●●● ●● ● ● ●●● ● ●●● ● ●● ● ●● ●● ● ● ● ● ●●●● ● ●●●● ●● ●● ● ●●●●●● ●●● ●●● ●● ●●●● ●●●●● ● ● ●●● ● ●●● ●● ●●●●● ●● ●●● ●● ● ●●● ● ●● ●●●●●● ●●●● ● ● ●●● ●● ●●●● ●● ●● ● ● ●● ●●● ● ●●● ●● ●●●● ●●● ●● ●●● ●●●● ●● ●●● ● ●● ● ●● ●●● ● ●●●● ●●● ● ●●● ● ●●● ● ●●●●● ●●●● ● ●●● ●●● ● ●●● ●● ● − − Algor (cid:137)(cid:138)(cid:139)(cid:140) (cid:141)(cid:142)(cid:143)(cid:144)(cid:145)(cid:146)(cid:147)(cid:148)(cid:149)(cid:150)(cid:151)(cid:152)(cid:153)(cid:154)(cid:155)(cid:156)(cid:157)(cid:158)(cid:159)(cid:160)¡ (e) Progesterone. ●● ●● ●● ●● ● ●● ● ●●● ● ● ●●● ● ●● ● ●● ●●● ● ●●● ●● ● ● ●● ●● ●●● ●●●● ● ●●●● ●● ●● ●●●●●● ●● ●●● ●●● ●●● ●●●● ●●● ●●●● ●● ●● ●●● ●● ●●● ●● ● ● ●●●● ●● ●●●●●● ●● ●●●● ●● ● ●●● ●● ● ●●●● ●●● ●● ●●●● ●●● ●●●●● ●●●● ●●● ●●●●● ●● ●● ●●●●●●● ● ●● ●●●● ●● ● ●●● ●● ● ●●●●●●● ●● ●● ●●●● ●●●●●● ●●●● ●● ●●●●● ●● ●●●●● ●● ●●● ● ●● ●●●●●● ●●●● ●●● ● ●●●● ● ● ● ●● ●●● ● ●●● ●● ● ●● ● ●●●●●●● ●● ●●●● ●●●●●● ●● ●● ● ●●●● ●● ● ●●● ●●●●●● ●● ●●● ● ●●● ●● ● ●●● ●●●● ●●●●● ● ●● ●●● ●●● ● ●●● ●● ● ●●● ●●● ●● ● ●●●●●●●● ● ●●● ●●●●● ● ● ●●●● ●● ● ●●● ● ●●●● ●●● ●● ●●● ● ●●●● ●● ●●● ● ● − − − Algor ¢£⁄¥ ƒ§¤'“«‹›fifl(cid:176)–†‡·(cid:181)¶•‚„” (f) DHEA.