Statistical Search Space Reduction and Two-Dimensional Data

May 13, 2006 - ... time-of-flight mass spectrometry (UPLC Tof MS) data on biofluids is presented and exemplified using a study on hydrazine-induced li...
0 downloads 3 Views 348KB Size
Download PDF
Anal. Chem. 2006, 78, 4398-4408

Statistical Search Space Reduction and Two-Dimensional Data Display Approaches for UPLC-MS in Biomarker Discovery and Pathway Analysis Derek J. Crockford,*,† John C. Lindon,† Olivier Cloarec,† Robert S. Plumb,‡ Stephen J. Bruce,† Severine Zirah,† Paul Rainville,‡ Chris L. Stumpf,‡ Kelly Johnson,‡ Elaine Holmes,† and Jeremy K. Nicholson†

Biological Chemistry, Division of Biomedical Sciences, Sir Alexander Fleming Building, Imperial College, London, SW7 2AZ, U.K., and Waters Corporation, 34 Maple Street, Milford, Massachusetts 01757

A new analytical strategy for biomarker recovery from directly coupled ultra-performance liquid chromatography time-of-flight mass spectrometry (UPLC Tof MS) data on biofluids is presented and exemplified using a study on hydrazine-induced liver toxicity. A key step in the strategy involves a novel procedure for reducing the spectroscopic search space by differential analysis of cohorts of normal and pathological samples using an orthogonal projection to latent structures discriminant analysis (O-PLS-DA). This efficiently sorts principal discriminators of toxicity from the background of thousands of metabolic features commonly observed in data sets generated by UPLC-MS analysis of biological fluids and is thus a powerful tool for biomarker discovery. There is great interest in the development of new analytical techniques and strategies for biomarker discovery in a range of -omics research programs.1,2 Gas chromatography coupled to mass spectrometry (GC-MS) has proven to be a valuable tool in metabonomic and metabolomic studies owing to its high sensitivity and analytical robustness and the ready availability of reference libraries for chemical identification.3 However, for the analysis of more polar metabolites GC-MS is less suitable as it requires its analytes to be volatile. While this is readily achieved using various derivatization methods, these require additional sample preparation steps that are undesirable in high-throughput applications, and discovery capability for new metabolite structures is limited. High-performance liquid chromatography coupled to mass spectrometry (HPLC-MS) can handle less volatile analytes directly but also has limitations, the most taxing of which is differential * Corresponding author. Tel: +44-207-594-3107. E-mail d.crockford@ imperial.ac.uk. † Imperial College. ‡ Waters Corporation. (1) Dumas, M. E.; Canlet C.; Debrauwer, L.; Martin, P.; Paris, A. J. Proteome Res. 2005, 4, 1485-1492. (2) Nicholson, J. K.; Connelly, J.; Lindon, J. C.; Holmes, E. Nat. Rev. Drug Discovery 2002, 1, 153-161. (3) Bino, R. J.; Hall, R. D.; Fiehn, O.; Kopka, J.; Saito, K.; Draper, J.; Nicolau, B.; Mendes, P.; Roessner-Tunali, U.; Beale, M. H.; Trethewey, R. N.; Lange, B. M.; Syrkin Wurtele, E.; Sumner, L. W. Trends Plant. Sci. 2004, 9 (9), 418-425.

4398 Analytical Chemistry, Vol. 78, No. 13, July 1, 2006

ion suppression.4 The ion suppression problem arises due to the relatively large sample injections required, so that the most ionizable components of a complex mixture take a disproportionate fraction of the available ion current, and results in loss of signal from the less ionizable components. This can make results nonquantitative and poorly reproducible, particularly for samples of high inherent complexity and variable composition, where the context is exploratory and the measurement problem is not wellcharacterized. Metabolic profiling methods rely heavily on multivariate analysis techniques such as Principal Components Analysis (PCA) and Projection to Latent Structures (PLS) to handle the massive quantities of data typically arising from applying information-rich analytical techniques to large sample sets, and these have been supplemented recently by methods based on statistical correlation.5,6 However, to achieve optimal characterization of the samples and enable efficient biomarker detection, GC-MS and HPLCMS data sets generally require several steps of data pre-processing prior to any multivariate analysis. Considerable effort has gone into the development of algorithms to align signals from the same compound in data acquired from different samples.7,8 This problem has been partially solved for NMR spectroscopic data (which can be subject to peak position variation problems) by dividing the spectrum into frequency windows and integrating signal intensity within these segments9,10 or by means of advanced PCA techniques,11 and analogous approaches have been applied to MS data. (4) Avery, M. J. Rapid. Commun. Mass Spectrom. 2003, 17 (3), 197-201. (5) Cloarec, O.; Dumas, M. E.; Craig, A.; Barton, R. H.; Trygg, J.; Hudson, J.; Blancher, C.; Gauguier, D.; Lindon, J. C.; Holmes, E.; Nicholson, J. K. Anal. Chem. 2005, 77, 1282-1289. (6) Crockford, D. J.; Holmes, H.; Lindon, J. C.; Plumb, R. S.; Zirah, S.; Bruce, S. J.; Rainville, P.; Stumpf, C. L.; Nicholson, J. K. Anal. Chem. 2006, 78 (2), 363-371. (7) Jonsson, P.; Johansson, A. I.; Gullberg, J.; Trygg, J.; A, J.; Grung, B.; Marklund, S.; Sjo¨stro¨m, M.; Antti, H.; Moritz, T. Anal. Chem. 2005, 77, 5635-5642. (8) Idborg-Bjo ¨rkman, H.; Edlund, P.-O.; Kvalheim, O. M.; Schuppe-Koistinen, I.; Jacobsson, S. P. Anal. Chem. 2003, 75, 4784-4792. (9) Farrant, R. D.; Lindon, J. C.; Rahr, E.; Sweatman, B. C. J. Pharm. Biomed. Anal. 1992, 10, 141-144. (10) Holmes, E.; Foxall, P. J.; Nicholson, J. K.; Neild, G. H.; Brown, S. M.; Beddell, C. R.; Sweatman, B. C.; Rahr, E.; Lindon, J. C.; Spraul, M. Anal. Biochem. 1994, 220, 284-296. 10.1021/ac060168o CCC: $33.50

© 2006 American Chemical Society Published on Web 05/13/2006

Figure 1. Plot of log (intensity) for a UPLC-MS data set on a single urine sample, showing approximately 37 000 peaks.

In one study, a series of prominent marker peaks12 has been used to divide chromatograms into time windows in which the retention time scale is then adjusted piecewise by an interpolation scheme. Another approach13 divided chromatograms into time windows separated by regions showing a baseline response in all samples, constructed the equivalent direct injection mass spectra, and then extracted the m/z values that varied significantly using alternating regression. Recently, the new chromatographic technique of ultraperformance liquid chromatography has been developed, and it has been coupled to mass spectrometry (UPLC-MS) and applied to metabonomic studies.14,15 UPLC can provide better chromatographic resolution than HPLC due to the use of a reduced particle size in the column packing and allow for more direct and reliable multivariate analysis of metabonomic datasets derived from polar analytes. Ongoing studies show that mammalian metabolic complexity and interactions are greater than has been previously supposed due to combinatorial interactions of mammalian and symbiotic gut microbiotal systems that lead to the production of extensive series of microbe-mammalian co-metabolites that are often not accounted for in mammalian biochemical pathways.16,17 As it would (11) Stoyanova, R.; Nicholls A. W.; Nicholson J. K.; Lindon J. C.; Brown T. R. J. Magn. Reson. 2004, 170 (2), 329-35. (12) Yang, J.; Xu, G.; Zheng, Y.; Kong, H.; Wang, C.; Zhao, X.; Pang, T. J. Chromatogr. A 2005, 1084 (1-2), 214-221. (13) Jonsson, P.; Bruce, S. J.; Moritz, T.; Trygg, J.; Sjo ¨stro ¨m, M.; Plumb, R.; Granger, J.; Maibaum, E.; Nicholson, J. K.; Holmes, E.; Antti, H. Analyst 2005, 130, 701-707. (14) Wilson, I. D.; Nicholson, J. K.; Castro-Perez, J.; Granger, J. H.; Johnson, K. A.; Smith, B. W.; Plumb, R. S. J. Proteome Res. 2005, 4, 591-598. (15) Plumb, R. S.; Granger, J. H.; Stumpf, C. L.; Johnson, K. A.; Smith, B. W.; Gaulitz, S.; Wilson, I. D.; Castro-Perez, J. Analyst 2005, 130 (6), 844- 849. (16) Nicholson, J. K.; Connelly, J.; Lindon, J. C.; Holmes, E. Nat. Rev. Drug Discovery 2002, 1, 153-161.

appear that pathological processes caused by drugs, diseases, and genetic defects influence either directly or indirectly this cometabolic axis,18 there is much potential for novel biomarker discovery in such co-metabolic datasets. Partly as a result of these co-metabolic interactions, biofluids such as urine contain tens of thousands of metabolites, and many of these may be revealed by modern analytical techniques such as UPLC-MS.14,15 Paradoxically, this creates a new statistical problem related to biomarker discovery that concerns the recovery of key discriminating biomarker information in the background of thousands of detectable metabolic variables (Figure 1), the majority of which do not actually contribute to class separation or diagnosis. In this paper, we approach the processing of UPLC-MS data acquired on urine samples from an exemplar toxicology study by applying O-PLS-DA to matrices formed by binning the raw data over windows in both retention time and m/z. The sample set analyzed contains urine specimens from both control animals and from those exhibiting a strong response to administration of the model toxin hydrazine, which elicits a steatotic effect in the liver.19 The aim of the study is to demonstrate the stratification of the large number of UPLC-MS metabolite peaks according to their relevance to toxicity (and in principle to any other pathological or disease state). In this way, subsequent metabolite identification and quantitation efforts can be focused on the most important peaks. The approach efficiently avoids the vexed problems of deconvolution or peak alignment between samples, and the calculation of relative concentrations before applying the multivariate analysis, and allows significant peaks that are close to the (17) Nicholson, J. K.; Holmes, E.; Lindon, J. C.; Wilson, I. D. Nat. Biotechnol. 2004, 22, 1268-74. (18) Nicholson, J. K.; Holmes, E.; Wilson, I. D. Nat. Rev. Microbiol. 2005, 3, 431-8. (19) Jenner, A. M.; Timbrell, J. A. Arch. Toxicol. 1994, 68, 240-245.

Analytical Chemistry, Vol. 78, No. 13, July 1, 2006

4399

Figure 2. Analytical strategy (σ ) O-PLS weight).

detection limit to be recovered statistically without time-consuming inspection of individual spectra. We show that, given a suitable choice of window size, useful peak selection can be obtained rapidly and that these peaks can be treated as biomarker candidates for further analysis and structure characterization. The use of a cross-validated method like O-PLS-DA provides greater confidence than simple correlation analysis of variable against class. In addition, the method detects and models separately any variation that is orthogonal to class, so this can be separately analyzed and does not affect the class dependency results. A summary of the analytical strategy employed is given in Figure 2. Following binning of data, the selection step consists of filtration of variables by their weight in an O-PLS-DA model, which gives a measure of correlation between peak intensity and dosage. These weights, color coded by value, are displayed in twodimensional diagrams labeled by retention time bin and m/z bin, which are expanded as required. Once the discriminating peaks are found, a library search3 is conducted to identify the possible compounds involved. Given that a single bin may contain more than one peak and that accurate masses and retention times are required, it is still necessary to interrogate the original UPLCMS spectra to find the exact values for the biomarker candidates. In this study, the library search against accurate m/z was done manually for the list of peaks of interest, but this could be automated. A final step compares MS/MS fragment patterns of pure compounds with urine samples from the study to either confirm or exclude candidate compounds highlighted by the search procedure. Overall, we demonstrate the use of an efficient O-PLS-DA statistical approach to optimize extraction of candidate biomarker data from highly complex biofluid UPLC-MS data sets. METHODS Animal Studies and Sample Collection. Rat urine samples were obtained from a study of a model toxin (hydrazine) 4400

Analytical Chemistry, Vol. 78, No. 13, July 1, 2006

conducted as part of COMET,20 a major collaborative toxicology project. All animal experiments were conducted according to appropriate national guidelines. The hydrazine study was chosen because it elicits a reliable response with well-defined toxic effects.21 In this study male Sprague-Dawley rats were randomly allocated to groupsscontrol, low, and high dose (n ) 10 per group). Here we only consider the high dose group, which received hydrazine hydrochloride in 0.9% saline administered orally at 90 mg/kg. The time of dosing was defined as zero hours. Rats were housed in individual metabolic cages under controlled temperature, humidity, and light cycles. Urine samples were collected for time intervals from -16 to 168 h.22 Data used for this paper were obtained from 30 samples, all from the same 10 high dose animals, of which 15 corresponded to predose time points (-16 and 0 h), which are effective controls, and 15 corresponded to time points (48 and 72 h) when a strong toxic effect was present. UPLC-MS Samples and Spectra. UPLC-MS analysis of urine was performed using a Waters23 ACQUITY UPLC Separations Module coupled to a Waters MicroMass Q-Tof Premier mass spectrometer equipped with an electrospray interface and leucine enkephalin lock-spray. Urine samples were diluted 4-fold with LCgrade water, and 2 µL injections were made from a glass autosampler vial maintained at 4 °C into a 10 cm × 2.1 mm Waters ACQUITY BEH C18 1.7 µm UPLC column maintained in an oven at 40 °C. The column was eluted under gradient conditions at a flow rate of 500 µL/min; mobile phase component A consisted of 0.1% aqueous formic acid, and phase B consisted of 0.1% formic acid in acetonitrile. The composition was linearly increased from 0% B to 20% B over 4 min, then to 95% B over the next 5 min, returned to 100% A over 1 min, and then reequilibrated over a final 2.5 min prior to injection of the next sample. The mass spectrometer was operated in positive ion mode with a capillary voltage of 3 kV, cone voltage of 25 V, nebulizer gas flow of 600 L/h, desolvation temperature of 250 °C, and source temperature of 120 °C. The instrument was set to acquire over the mass range m/z 50-850 with acquisition time of 300 ms and an interacquisition delay of 100 ms. MS/MS Experiments. HPLC-MS/MS analysis of solutions of pure compounds and mixtures of urine samples was performed on a Waters 2795 HPLC Separations Module linked to a MicroMass Q-Tof Micro mass spectrometer with an electrospray interface and leucine enkephalin lock-spray. Reference compound solutions were prepared by 100-fold dilution of 5 mg/mL stock solutions of pure compounds. The solvent was 0.1% aqueous formic acid, either pure or mixed with 0.1% formic acid in acetonitrile to a maximum proportion of 1:1, depending on the (20) Lindon, J. C.; Nicholson, J. K.; Holmes, E.; Antti, H.; Bollard, M. E.; Keun, H.; Beckonert, O.; Ebbels, T. M.; Reily, M. D.; Robertson, D.; Stevens G. J.; Luke, P.; Breau, A. P.; Cantor, G. H.; Bible, R. H.; Niederhauser, U.; Senn H.; Schlotterbeck, G.; Sidelmann, U. G.; Laursen, S. M.; Tymiak, A.; Car, B. D.; Lehman-McKeeman, L.; Colet J.-M.; Loukaci, A.; Thomas, C. Toxicol. Appl. Pharmacol. 2003, 187, 137-146. (21) Nicholls, A. W.; Holmes, E.; Lindon, J. C.; Shockcor, J. P.; Duncan Farrant, R.; Haselden, J. N.; Damment, J. P.; Waterfield, C. J.; Nicholson, J. K. Chem. Res. Toxicol. 2001, 14, 975-978. (22) Bollard, M. E.; Keun, H. C.; Beckonert, O.; Ebbels, T. M. D.; Antti, H.; Nicholls, A. W.; Shockcor, J. P.; Cantor, G. H.; Stevens, G.; Lindon, J. C.; Holmes, E.; Nicholson, J. K. Toxicol. Appl. Pharmacol. 2005, 204 (2), 13551. (23) WATERS. http://www.waters.com.

solubility of the compound. Lesser dilutions were prepared as necessary according to the strength of MS signal obtained. Urine samples consisted of mixtures of equal parts of all urine specimens analyzed by UPLC-MS, the 15 predose specimens forming a single “control” analyte and the 15 maximum response specimens forming a single “maximum response” analyte. HPLC was conducted so as to broadly match the above UPLC conditions, but only the UPLC retention times are quoted in this paper. The mass spectrometer parameters employed were substantially as described above. MS/MS was triggered by detection of an ion at the nominal m/z of the expected parent ion. Fragmentation was performed in a collision cell. Data Analysis. UPLC-MS data were converted from the standard commercial format into text files using the Waters DataBridge utility. These files were read into MATLAB24 where the data for each sample were integrated (“binned”) to produce two-dimensional matrices of MS intensity. The log of one of these matrices is given in Figure 1. The matrices are stacked over sample number to produce the three-dimensional data matrix X, where the dimensions are labeled by retention time bin, m/z bin, and sample number. The bins are converted into their center values (in minutes and Thomson) when results are represented graphically. No normalization or peak matching steps were applied. All data were mean-centered before pattern recognition. The pattern recognition method used was O-PLS-DA.25 This derives separate models for the variation in X (the UPLC-MS data) and Y (the dosage class), linked by the scores matrix T. The part of the variation in X that is orthogonal to the variation in Y is given by TYoscPYosc, where TYosc is the orthogonal scores matrix and PYosc is the corresponding loading. An O-PLS-DA model can be written

Model of X: X ) TWt + TYoscPYosc + E Model of Y: Y ) TCt + F Prediction of Y: Y′ ) TCt

where W and C are joint orthonormal loading matrices, E and F are the respective residual matrices, Y′ is the predicted value of Y, and the subscript t indicates that the matrix is transposed. This method provides similar prediction to traditional PLS-DA but is more directly interpretable because the variation in X orthogonal to Y is modeled separately. In particular, the O-PLS-DA loading for a given variable can be readily converted to an O-PLS-DA weight σ, and in the case of two classes σ2 gives the proportion of the variable’s variance that it has in common with the class descriptor.26 As σ2 shows considerable variation, its value allows us to filter for variables (i.e., retention time and m/z intervals) (24) MATLAB. http://www.mathworks.com/products/matlab. (25) Trygg, J.; Wold, S. J. Chemom. 2003, 18, 53-64. (26) Cloarec, O.; Dumas, M. E.; Trygg, J.; Craig, A.; Barton, R. H.; Lindon, J. C.; Nicholson, J. K.; Holmes, E. Anal. Chem. 2005, 77, 517-526.

that are important to class discrimination and, hence, likely to be of biological significance. The σ2 value can be converted to a t-statistic, and all the peaks were validated by t-test at p ) 0.01% (unless stated otherwise) against the hypothesis that the peak is selected by chance where no real relation to class exists. To test the validity of models against over-fitting, the prediction accuracy (or cross-validation parameter, Q2Y′) was computed as:

Prediction accuracy in Y:

∑(TC - Y) )1∑Y t

2

QY'

2

2

Other model metrics presented are the variation in Y explained by the model, R2Y, the variation in X explained by the model, R2X, and the explained variation in X that correlates with the variation in Y, R2Xcorr. The distance to model metric was also calculated together with a critical threshold value of 99% and was used to check for outlying samples that might lead to spurious results. As a comparison, the same mean-centered input data was univariance scaled and subjected to PCA analysis, and the loadings plots for the main pairs of components that distinguished between classes were inspected in an attempt to find the most significant variables. Computer and Software. All processing and pattern recognition was conducted on a Dell personal computer with a 2.8 GHz Xeon processor and 2 GB of RAM. Processing of HPLC-MS and UPLC-MS data was performed using MassLynx23 and in-house code written in MATLAB 6.5.1.24 RESULTS Identification of UPLC-MS Peaks Relevant to Toxicity. A typical intensity plot of the UPLC-MS signal for a urine sample is given in Figure 1, where the intensity is denoted by color and has been converted to a logarithmic scale in order to accommodate the wide differences in signal levels. Given that urine contains a large number of small, polar compounds and because the C18 column employed is not optimal for polar analytes, the bulk of the signals are at low retention time and low m/z. The broad signal at 9.5 min retention time is due to highly retained material that is eluted off the column only when the mobile phase gradient approaches 100% acetonitrile. Owing to the sensitivity and resolution of UPLC-MS, there are a very large number of peaks, representing a formidable data interpretation challenge, but this can be considerably simplified by further analysis. The O-PLS-DA results are presented as plots of UPLC-MS peak position in the two-dimensional space of retention time and m/z with a color code corresponding to their σ2 value. To filter the peaks for relevance to discrimination, only those with σ2 exceeding a given threshold are included. For example, those peaks with σ2 exceeding 0.5 are plotted in Figure 3A (note that the highly retained material mentioned above no longer appears). The diagram offers a huge simplification over the raw UPLCMS data (Figure 1) and shows that the significant peaks are clustered at low retention times and masses, as would be expected for a urine analyte. The number of peaks present is still large but may be reduced further by increasing the σ2 threshold to 0.7 to obtain more strongly discriminating biomarker candidates (Figure Analytical Chemistry, Vol. 78, No. 13, July 1, 2006

4401

Figure 3. (A) Significant LC-MS peaks at a σ2 threshold of 0.5 (807 peaks). (B) Significant LC-MS peaks at a σ2 threshold of 0.7 (160 peaks).

3B). The red color indicates the increased σ2 values, and there are still a significant number of ionic compounds eluting close to the solvent front at 0.7 min. In addition, there are several ions that elute at multiple retention times, which might be explained as fragments of compounds that are structurally and, hence, biologically related. In the comparison PCA analysis, reasonable separation of classes was obtained in the score plots for PCA components 1 and 3, 2 and 3, and 3 and 4, but in each case the corresponding loadings plot yielded a dense pseudo-continuum of points that was 4402 Analytical Chemistry, Vol. 78, No. 13, July 1, 2006

very difficult to interpret in terms of variable importance (diagrams not shown here). In Figure 4A we expand the filtered region shown in Figure 3B and plot σ values so that we can distinguish whether the peaks increase or decrease with respect to dosage: those that increase are in blue while those that decrease in red. Horizontal “streaks” can arise if a peak lies close to the edge of a retention time bin, or if there are simply several significant peaks with the same mass, which is not unlikely given the unit m/z resolution employed. Setting a σ2 threshold of 0.9 and expanding further, as in Figure

Figure 4. (A) Expansion of Figure 3B color-coded by σ value to indicate effect of dosage. (B) Expansion of Figure 4A for σ2 threshold 0.9 showing nominal m/z (compare Table 1).

4B, allows us to select a reasonable number of peaks to investigate further using MS/MS. The peaks so obtained, along with their fragmentation patterns, are given Table 1. The resolution in Table 1 is finer than that of the O-PLS-DA model, and to obtain these values it was necessary to inspect carefully some of the raw spectra for various retention times. The Tof m/z values are derived from the original measurements made on a more accurate Tof instrument, while the QTof values come from subsequent MS/MS measurements. The QTof samples consisted of two mixtures, one consisting of urine taken from all

samples used in the high response class, and the other from all samples in the control class (see MS/MS Experiments section). Clearly, to identify peaks that increase with dosage we analyze the high response mixture, while to investigate those that decrease we analyze the control mixture. There is a possibility that a given bin may contain multiple peaks, so there is some potential for confusion over which peaks are significant. However, it is straightforward to re-run the O-PLSDA model at a finer resolution over a more limited range of retention time and m/z and hence distinguish exactly which peaks Analytical Chemistry, Vol. 78, No. 13, July 1, 2006

4403

Table 1. MS Results for Ions with Largest Absolute σ2 Value According to O-PLS-DA RT/min

Tof m/z

QTof MS/MS fragment m/z

putative identity

σ2

1.27 1.12 0.79 0.71 1.27 0.67 0.88 0.95 0.69 1.52 1.25

402.10 193.04 164.07 351.19 210.06 245.16 260.09 358.11 311.11 290.12 147.07

phosphopantothenyl cysteine (402.09) citrate, isocitratea (both 193.03) 4-hydroxyglutamic acid (164.06) prostaglandin-E3a (351.21) unknown N1-acetylsperminea (245.23), uridinea (245.08), biotin (245.10) glucose-6-phosphate, glucosamine-6/1-phosphate (both 260.05) L-stercobilinogen fragment 5-ribosylpyrophosphate, ribulose-1,5-bisphosphate (both 310.99) N-succinyl-2-amino-6-oxopimelate (290.09) lysinea (147.10), glutaminea (147.07), 2-dehydropantoate (147.06)

0.94 0.93 0.92 0.91 0.91 0.91 0.90 0.90 0.90 0.90 0.90

0.69 0.67 0.71 0.96 0.96 2.64 1.57

162.06 116.07 203.13 98.05 265.05 205.07 159.06

weak signal, ambiguous weak signal, ambiguous 164.13, 147.04, 123.05, 119.07 351.16, 132.08 not acquired 245.16, 185.14, 132.08 not acquired not acquired 311.11, 245.14 weak signal, ambiguous 147.09, 146.09, 129.07, 105.07, 87.04, 84.07 162.09, 144.05, 116.06, 98.06 116.08, 98.07, 70.09 129.13, 112.12 98.08, 70.09, 55.05 265.09, 219.07, 136.05 205.10, 158.09, 144.07 159.02, 116.04, 98.99, 61.06

0.84 0.83 0.80 0.80 0.70 0.69 0.68

1.57

201.09

2-aminoadipateb (162.08) prolinea (116.08) spermineb (203.22) unknown phenylacetylglutamine tryptophana (205.10) propionyl carnitine fragment, allantoina (159.05), dihydroorotic acida (159.03), propionyl carnitine fragment

0.74 2.99 0.96

132.06 190.03 217.09

creatineb (132.08) kynurenic acidb (190.06) glucuronate-Na+ adduct (217.03)

0.64 0.64 0.63

0.74 0.96 1.91 1.57

90.05 219.08 144.06 218.11

L-alanine,b

0.63 0.56 0.52 0.50

a

201.09, 183.09, 159.08, 114.07, 70.09, 60.07 132.05, 90.06 190.06, 172.05, 162.07, 144.05 217.13, 201.09, 183.08, 159.08, 115.10 90.06 sample too dilute 144.07, 98.08, 84.98 218.13, 201.09, 159.08, 70.09

β-alanine, sarcosineb (all 90.06) N-acetylserotonin (219.10) unknown propionyl carnitine (218.13)

0.67

Candidates inconsistent with observed MS/MS (see Table 2). b Candidates consistent with observed MS/MS of pure compounds.

Figure 5. Expansion of peak at nominal m/z 144 for a model of retention time range 0.6-0.9 and m/z range 100-200 at an m/z bin of 0.01 and σ2 threshold 0.3 (p ) 0.05). The range contains peaks that correlate oppositely with dose and tend to cancel at nominal m/z resolution. The less intense isotopic peaks are clearly visible.

in a bin are statistically relevant. This has been done for the polar retention times between 0.6 and 0.9 min with an m/z range between 100 and 200, and an m/z bin size of 0.01. Concentrating on a peak at nominal m/z 144, for example, which was not selected 4404 Analytical Chemistry, Vol. 78, No. 13, July 1, 2006

in the full model (Table 1), we see in Figure 5 that the range contains two strongly correlated peaks that vary in opposite senses with dose and therefore tend to cancel at nominal m/z resolution. It is interesting to note that the smaller isotope peaks are visible

Table 2. Nominal m/z Values for Selected Pure Compounds Using QTof MS/MS at Collision Energies Giving Best Match to Table 1 compound

collision expected energy/eV parent

2-aminoadipate

10

162

allantoin β-alanine citric acid isocitric acid creatine dihydroorotic acid

7.5 10 10 10 10 10

159 90 193 193 132 159

glutamine kynurenic acid

10 15

147 190

L-alanine

lysine N1-acetyl spermine

10 10 10

90 147 245

proline prostaglandin-E3

15 10

116 351

sarcosine spermine

10 15

90 203

tryptophan

10

205

uridine

10

245

ions detected (with count) 162(1060), 144(1090), 116(660), 98(500) 116(28), 61(12) 90(110), 72(161) weak signal weak signal 132(660), 90(400) 159(11), 131(5), 113(4), 87(9), 71(3) 147(18), 130(116), 84(25) 190(197), 172(141), 162(282), 144(381) 90(38), 89(3), 62(2) 147(6), 130(24), 84(4) 245(263), 171(39), 129(47) 116(23), 70(75) decomposed during analysis 90(105), 72(3) 203(28), 129(291), 112(205), 84(12) 205(32), 188(247), 146(10) 245(7), 133(24), 113(243)

for both monoisotopic peaks, showing that the model discriminates on relevance to dose and can detect less intense but relevant peaks alongside the more intense. The σ2 value for each biomarker candidate is given in the last column of Table 1. All candidates are significant at p ) 0.01% against the hypothesis that they arise by pure chance. A sensible strategy for identification of unknowns is to begin with the highest values and progress downward (i.e., to work down Table 1). The possible identities cited in these tables were obtained by manual searching of various resources, but this process could be automated if a comprehensive database were available, giving the method potential application in a high-throughput environment. The nominal fragment m/z values for some of the proposed matches, sufficient to distinguish compounds by their fragmention patterns, have been given in Table 2 and served to confirm or rule out some of these identities, as indicated. Table 1 contains several compounds that are well-known metabolites from NMR studies,21 such as 2-aminoadipate, creatine, and the alanine isomers. However, due to the complementary nature of NMR spectroscopy with MS there are many that for sensitivity reasons are not frequently detected in NMR, such as spermine and kynurenic acid. There are many compounds that warrant further investigation, and in a comprehensive study peaks that elude database searches or are not confirmed by MS/MS could be isolated and subject to more structurally informative analytical techniques such as 2D NMR spectroscopy. In addition, only the most correlated peaks are given in Table 1, but further studies that progress down the significance scale might be expected to yield further interesting and unsuspected compounds. Sensitivity of O-PLS-DA Model to Bin Size and Orthogonal Components. The prediction accuracy obtained by 7-fold crossvalidation of the model for various combinations of bin size is illustrated in Table 3. The search for relevant metabolites was

Table 3. Variation of O-PLS-DA Prediction Accuracy (Q2Y′) with Binning Parameters RT bin size/min m/z bin size

0.1

1.0

12.5

1 2 5 10 20 50 100 200 500 1000

0.89 0.89 0.89 0.88 0.87 0.82 0.76 0.76 0.65 0.65

0.88 0.88 0.86 0.82 0.78 0.66 0.58 0.54 0.51 0.48

0.83 0.79 0.69 0.40 0.08 -0.27 -0.26 -0.24 -0.17 -0.17

done with a retention time bin of 0.1 min and a m/z bin of 1, which was the highest practicable resolution using a high specification personal computer, and yielded a very high accuracy of 0.89. The diagram indicates that at this resolution the prediction accuracy is very stable to loss of resolution in retention time or m/z individually, showing that they both contribute significant information to the model. In fact, setting the retention time window to 12.5 min, the equivalent of a direct injection mass spectrum, still gives a prediction accuracy of 0.83. With m/z window 1000, the equivalent of a total ion count chromatogram, we obtain a prediction accuracy of 0.65. The curve for the 12.5 min retention time window demonstrates that the pseudo-direct injection spectra cease to have any predictive power when the m/z window reaches 20. The general behavior of the model is presented in Figure 6, where panel A characterizes the model as a function of the number of components of X included which vary orthogonally to Y, and panel B measures the distance to model27 for the data from each sample, with a 99% confidence threshold for outliers marked as a dashed line. It can be seen that the model contains no outlying data points at this threshold; therefore, we do not expect any biomarker candidate to arise due to the spurious effects of outlying samples. Panel A is especially interesting because it shows that the prediction accuracy Q2Y′ is scarcely affected by number of orthogonal components. The explained variations in Y and X, R2Y and R2X, respectively, do increase, but without an accompanying increase in prediction accuracy this is indicative of over-fitting. Like the prediction accuracy, the explained variation in X that correlates to Y, R2Xcorr, is also unenhanced by additional model components, and in fact decreases slightly. Overall, the diagram shows that, for the choice of samples in this study, there are no significant sources of orthogonal variation, and for this reason the results presented are for a model with no orthogonal components included. DISCUSSION General Approach. LC-MS-based metabolic profiling studies as applied to the characterization of disease and other pathological processes are becoming increasingly common. Yet many of these studies identify only one or two key molecules as candidate biomarkers or suggest putative markers based on a single m/z (27) Eriksson, L.; Johansson, E.; Kettaneh-Wold, N.; Wold, S. Multi- and Megavariate Data Analysis; Umetrics: Switzerland, 2001; ISBN 91-9737301-X.

Analytical Chemistry, Vol. 78, No. 13, July 1, 2006

4405

Figure 6. O-PLS-DA model metrics for bins of retention time 0.1 min and m/z 1.0 (Q2Y′ ) prediction accuracy; R2Y ) explained variation in Y; R2X ) explained variation in X; R2Xcorr ) variation in X that correlates to Y).

ratio and library match without confirmation of structure. A onestep LC-MS analysis followed by PCA is not sufficient to identify robust sets of biomarker candidates and may be unduly influenced by the ionizable metabolites present in high concentration, which can overshadow the contribution of low concentration metabolites that may be more diagnostic of a particular pathology. The strategy suggested here offers an efficient solution for reducing the search space of biomarker candidates while selecting metabolites based more on their correlation with pathology rather than gross concentration or propensity for ionization. In addition, the inbuilt cross-validation in the O-PLS-DA algorithm protects against overfitting the data and selection of spurious peaks. Having selected a relatively small number of metabolites which are significant in discriminating between the two classes, in this case samples from hydrazine-affected and control animals, it is feasible to implement a second analytical phase where MS/MS of analytes allows fragment pattern matching to confirm or eliminate librarybased metabolite hits. We have used a combination of LC-MS followed by MS/MS to provide information on the identity of some of the metabolites characterizing hydrazine toxicity (Table 4). In the following sections we demonstrate that this strategy is capable of capturing significant biomolecular and pathway information. In this discussion we have discounted those compounds ruled out by the MS/ MS analysis step but included those that were confirmed by MS/ MS or for which MS/MS results were unavailable. Inhibition of Amine Oxidases and Transaminases. Hydrazine and its derivatives are known to inhibit multiple enzymes, including amine oxidases28 and tranaminases.21 This is consistent (28) Gilvarg, C. Annu. Rev. Biochem. 1961, 239-268.

4406 Analytical Chemistry, Vol. 78, No. 13, July 1, 2006

with several of the perturbed urinary metabolites observed (Table 4). For example, serotonin sulfate has been detected in the urine of humans treated with the hydrazine derivative iproniazid, consistent with monoamine oxidase inhibition.28 Moreover, an increased urinary excretion of spermine may result from the inhibition of spermine-N1-acetyltransferase, a key enzyme in the catabolism of spermine.29 Similarly, there is an increase in the level of 4-hydroxyglutamate in the urine of the hydrazine-treated animals. The conversion of this compound to 2-oxoglutarate is mediated by 4-hydroxyglutamate transaminase, an enzyme that is inhibited by phenylhydrazine.30 Hence this increase is also consistent with hydrazine-induced transaminase inhibition. The elevation of 2-aminoadipate (2-AA), a well-known biomarker of hydrazine toxicity, can be explained by transaminase inhibition leading to disruption of the lysine catabolic pathway. This pathway relies on the enzyme 2-AA aminotransferase,21 which is inhibited by hydrazine, and this is associated with the neurotoxic effects of hydrazine due to pathway disruption in the hippocampus.31 In the brain, reduced kynurenic acid levels have been linked to the elevation of 2-aminoadipate caused by hydrazine.21 LCMS studies of kidney toxins32,33 have identified a reduction in urinary kynurenic acid in response to dosing. However, it has been established that levels of both kynurenic acid and serotonin in the brain increase in response to stress in the rat,34 which is more (29) Casero, R. A.; Pegg, A. E. FASEB J. 1993, 7, 653-661. (30) Goldstone, A.; Adams, E. J. Biol. Chem. 1962, 237, 3476-3485. (31) Wu, H. Q.; Ungerstedt, U.; Schwarcz, R. Eur. J. Pharmacol. 1995, 281, 55-61. (32) Lenz, E. M.; Bright, J.; Knight, R.; Wilson, I. D.; Major, H. Analyst 2004, 129, 535-541. (33) Lenz, E. M.; Bright, J.; Knight, R.; Wilson, I. D.; Major, H. J. Pharm. Biomed. Anal. 2004, 35, 599-608.

Table 4. Effect of Hydrazine Dosing and Possible Pathway Involvement for Some Ions Identified in Table 1 m/z 402.10 358.11 351.19 311.11 290.12 265.05 260.09 245.16 219.08 218.11 217.09 203.13 190.03 164.07 162.06 147.07 132.06 90.05 a

putative identity

effect of dosinga

pathway involvement

phosphopantothenyl cysteine fragment prostaglandin-E3 5-ribosylpyrophosphate/ribulose-1,5-bisphosphate N-succinyl-2-amino-6-oxopimelate phenylacetylglutamine glucose-6-phosphate/glucosamine-6-phosphate biotin N-acetylserotonin propionyl carnitine glucuronate-Na+ adduct spermine kynurenic acid 4-hydroxyglutamic acid 2-aminoadipate 2-dehydropantoate creatine L-alanine/sarcosine

V v V V V v v v v v v v v v v V v v

coenzyme-A biosynthesis bile pigment excretion antiinflammatory hormone purine biosynthesis in liver uncertain uncertain nitrogen excretion energy metabolism uncertain carboxylase cofactor neurotransmitter fatty acid metabolism detoxification cellular metabolism regulation neural regulation arginine & proline metabolism lysine catabolism coenzyme-A biosynthesis biomarker of liver disfunction amino acid urea cycle

L-stercobilinogen

From sign of σ value: v ) elevated by dosing; V ) depleted by dosing.

consistent with the observed increases in urine with dosage of the liver toxin hydrazine. Metabolic Perturbations Associated with Liver Toxicity. In addition to inhibition of transaminase enzymes, several metabolic changes in response to hydrazine administration were indicative of liver injury. One of the pathways reported to be disrupted in steatosis is the tryptophan f niacin pathway. The administration of hydrazine caused an increase in 5-ribosylpyrophosphate, which is a key intermediate of this pathway.35 This metabolite has not been detected in NMR-based metabolic profiling studies on hydrazine toxicity; however, other metabolic intermediates in the same pathway, such as N-methylnicotinate,6 have been identified. An increase in urinary creatine has been reported as a consistent feature of multiple types of liver injury and has been shown to be elevated after hydrazine dosing in NMR-based studies. This is also apparent here in the MS signature and is likely to result from increased cysteine synthesis associated with hydrazine-induced perturbation of the sulfur amino acid pathway.36,37 The first step in the detoxification of hydrazine in the liver is its acetylation,38 which requires coenzyme-A (CoA), and at the hydrazine doses used in this study could be expected to cause its depletion. Given that 2-dehydropantoate and phosphopantothenylcysteine (PPC) are both involved in the biosynthesis of CoA39 the increased demand for CoA production could lead to lower free levels of these compounds and hence reduced concentrations in urine. It is interesting to note that PPC is also (34) Pawlak, D.; Takada, Y.; Urano, T.; Takada, A. Brain Res. Bull. 2000, 52 (3), 197-205. (35) Fukuwatari, T.; Morikawa, Y.; Sugimoto, E.; Shibata, K. Biosci. Biotechnol. Biochem. 2002, 66 (6), 1196-1204. (36) Kenyon, S. H.; Waterfield, C. J.; Asker, D. S.; Kudo, M.; Moss, D. W.; Bates, T. E.; Nicolaou, A.; Gibbons, W. A.; Timbrell, J. A. Biochem. Pharmacol. 1999, 57, 1311-1319. (37) Clayton, T. A.; Lindon, J. C.; Everett, J. R.; Charuel, C.; Hanton, G.; Le Net, J.-L.; Provost, J.-P.; Nicholson, J. K. Arch. Toxicol. 2004, 78, 86-96. (38) Preece, N. E.; Nicholson, J. K.; Timbrell, J. A. Biochem. Pharmacol. 1991, 41 (9), 1319-1324. (39) KEGG. http://www.genome.ad.jp/kegg/pathway/map/map00770.html.

involved in the synthesis of taurine,40 the excretion of which is enhanced in hydrazine toxicity,21 and has been characteristic of several types of liver toxicity particular when modulated in conjunction with creatine.40,41 It has been suggested that creatine is synthesized as a consequence of increased cysteine synthesis in response to the requirement for taurine or glutathione as protectants against induced liver toxicity.37 Increased excretion of bilirubin is a known consequence of liver injury, such as hepatitis.42 Stercobilinogen is a breakdown product of bilirubin, which is absorbed in small quantities from the ileum via the hepatic circulation and then excreted in the urine. An increased urinary concentration of stercobilinogen is therefore completely consistent with liver toxicity. It is interesting to note that, although several other signals putatively assigned from the initial LC-MS data were not structurally confirmed by MS/MS fragmentation, they are nevertheless biologically consistent with the set of retained biomarkers. These compounds include N1-acetylspermine, uridine, lysine, glutamine, proline, and tryptophan. In particular, the “proline” fragment spectrum has one extra peak compared to authentic compound, so it could be an isomer with a related structure. However, since the chemical identities of these compounds could not be established by MS/MS analysis, it is not appropriate to comment further without additional experiments. CONCLUSIONS These results indicate the potential of UPLC-MS as a probe in metabonomics studies and show that the data produced can provide meaningful results from multivariate pattern recognition methods such as O-PLS-DA with a minimum level of preprocessing. It has been possible to build a model that predicts dosage class from data with a very high accuracy, and a number (40) Waterfield, C. J.; Turton, J. A.; Scales, M. D. C.; Timbrell, J. A. Arch. Toxicol. 1993, 67, 244-254. (41) Clayton, T. A.; Lindon, J. C.; Everett, J. R.; Charuel, C.; Hanton, G.; Le Net, J.-L.; Provost, J.-P.; Nicholson, J. K. Arch. Toxicol. 2003, 77, 208-217. (42) Paulev, P.-E. Textbook in Medical Physiology and Pathophysiology; Copenhagen Medical Publishers: Copenhagen, 1999.

Analytical Chemistry, Vol. 78, No. 13, July 1, 2006

4407

of significant peaks, both familiar and less familiar from NMR metabonomics studies, have been highlighted. The overall strategy allows the considerable efforts required to identify the unknown compounds in MS studies to be concentrated on the most significant unknowns first, and is amenable to automation. We believe that this approach will be of considerable value in handling the very large data sets that arise from modern chromatography and mass spectrometry, which can be efficiently analyzed for novel disease biomarkers and disrupted pathway activity. Furthermore, given that the structure of 2D gel electrophoresis data sets is

4408

Analytical Chemistry, Vol. 78, No. 13, July 1, 2006

similar to 2D LC-MS data sets, this method may also be applicable to biomarker detection in proteomic databases. ACKNOWLEDGMENT D.J.C. acknowledges support as part of a DTI Beacon Project. The authors are also grateful for access to urine samples obtained by the Consortium for Metabonomic Toxicology.20 Received for review January 25, 2006. Accepted April 17, 2006. AC060168O