Anal. Chem. 2006, 78, 4430-4442
Targeted Profiling: Quantitative Analysis of 1H NMR Metabolomics Data Aalim M. Weljie,†,‡ Jack Newton,† Pascal Mercier,† Erin Carlson,† and Carolyn M. Slupsky*†,§ Chenomx Inc., Edmonton, Alberta, Canada, and Metabolomics Research Centre, University of Calgary, Calgary, Canada
Extracting meaningful information from complex spectroscopic data of metabolite mixtures is an area of active research in the emerging field of “metabolomics”, which combines metabolism, spectroscopy, and multivariate statistical analysis (pattern recognition) methods. Chemometric analysis and comparison of 1H NMR1 spectra is commonly hampered by intersample peak position and line width variation due to matrix effects (pH, ionic strength, etc.). Here a novel method for mixture analysis is presented, defined as “targeted profiling”. Individual NMR resonances of interest are mathematically modeled from pure compound spectra. This database is then interrogated to identify and quantify metabolites in complex spectra of mixtures, such as biofluids. The technique is validated against a traditional “spectral binning” analysis on the basis of sensitivity to water suppression (presaturation, NOESY-presaturation, WET, and CPMG), relaxation effects, and NMR spectral acquisition times (3, 4, 5, and 6 s/scan) using PCA pattern recognition analysis. In addition, a quantitative validation is performed against various metabolites at physiological concentrations (9 µM-8 mM). “Targeted profiling” is highly stable in PCA-based pattern recognition, insensitive to water suppression, relaxation times (within the ranges examined), and scaling factors; hence, direct comparison of data acquired under varying conditions is made possible. In particular, analysis of metabolites at low concentration and overlapping regions are well suited to this analysis. We discuss how targeted profiling can be applied for mixture analysis and examine the effect of various acquisition parameters on the accuracy of quantification. Measurement of small-molecule metabolites, either endogenous or exogenous, provides a chemical “snapshot” of an organism’s metabolic state. Simultaneous characterization of numerous metabolites, or metabolic profiling, is an emerging area of analytical chemistry exemplified by the fusion of analytical metabolite measurement with pattern recognition chemometric statistical analysis1,2 and can be considered as a biological subset of general mixture analysis problems. Commonly referred to as “metabolomics” or “metabonomics”, this technique has wide applicability to a number of fields including medicine, plant sciences, toxicology, and food sciences.3 Within the realm of medicine, such * Corresponding author. Tel: 780-492-8919. Fax: 780-492-5329. E-mail:
[email protected]. † Chenomx Inc. ‡ University of Calgary. § Current address: Metabolomics Centre, Department of Medicine, University of Alberta, Edmonton, Alberta, Canada.
4430 Analytical Chemistry, Vol. 78, No. 13, July 1, 2006
profiling has provided both diagnostic and prognostic information using a variety of sources including serum, urine, cereberospinal fluid. and tissue extracts (See Table 1 from Holmes and Antti2 and references therein for representative examples.). The essence of the technology lies in distilling chemical information from complex spectroscopic information, such as nuclear magnetic resonance (NMR) spectroscopy or mass spectrometry data. NMR has been pursued as an analytical platform for metabolomics due to its inherent quantitative nature and the wealth of chemical information it can provide about nuclei that are “NMR visible”. Typical 1H 1D NMR experiments involve multiple samples (N from ∼10 to 100s or 1000s) with between ∼25-75 (tissue) and >200 (urine) NMR-visible metabolites in each sample. The sensitivity of NMR analysis to the chemical environment also brings certain challenges; subtle differences in pH, ionic strength, temperature, protein content, etc., between samples will cause differences in the NMR-detected peak position and line widths of a given metabolite (Figure 1). Furthermore, each metabolite is differentially sensitive to these effects, and in many compounds, the NMR resonances are affected independently, thus making a global correction infeasible for a mixture of multiple metabolites. Finally, in complex spectra there is a high degree of overlap in certain regions of the NMR spectrum, which hampers analysis. While this complexity can be reduced using selective data acquisition (e.g., ref 4), much more common is chemometric analysis of complete 1D or 2D spectra. The most widely employed method to overcome issues of peak and line width changes is a form of data reduction referred to as “spectral binning”.5,6 The spectrum is subdivided into a number of regions, and the total area within each bin is considered in further analysis. The assumption is that by considering regions of the spectra, as opposed to individual data points, minor peak shifts and line width differences for the same compound can be accounted for across samples. The size of the bins can be of a fixed width or variably sized using manual inspection (best optimization) or automated algorithms.7,8 For example, a typical 64k point NMR spectrum is reduced to ∼250 variables (K) upon fixed width binning at 0.04 ppm. In theory, each of these variables (Ki where i is the ith bin) will contain the same latent chemical (1) Eriksson, L.; Antti, H.; Gottfries, J.; Holmes, E.; Johansson, E.; Lindgren, F.; Long, I.; Lundstedt, T.; Trygg, J.; Wold, S. Anal. Bioanal. Chem. 2004, 380, 419-429. (2) Holmes, E.; Antti, H. Analyst 2002, 127, 1549-1557. (3) Griffin, J. L. Philos. Trans. R. Soc. London, Ser. B 2004, 359, 857-871. (4) Sandusky, P.; Raftery, D. Anal. Chem. 2005, 77, 7717-7723. (5) Gartland, K. P. R.; Beddell, C. R.; Lindon, J. C.; Nicholson, J. K. Mol. Pharmacol. 1991, 39, 629-642. (6) Anthony, M. L.; Sweatman, B. C.; Beddell, C. R.; Lindon, J. C.; Nicholson, J. K. Mol. Pharmacol. 1994, 46, 199-211. 10.1021/ac060209g CCC: $33.50
© 2006 American Chemical Society Published on Web 05/17/2006
Table 1. List of Urine Metabolites Studied and Their Concentrations (µM) in Each Sample sample
1-methylhistidine
citrate
creatine
creatinine
fucose
glycine
glycolate
hippurate
isoleucine
lactate
taurine
1 2 3 4 5 6
991 352 707 736 520 383
997 425 357 1289 780 1129
251 127 229 206 149 166
7677 5333 7967 6377 5874 7753
84 70 71 91 65 69
924 1273 455 1327 407 470
221 172 211 266 226 251
819 619 1659 1487 734 754
12 10 10 12 9 13
79 101 129 111 61 98
149 78 169 248 104 164
Figure 1. Structures for the compounds analyzed in this study. Note that compounds added for chemical shift indication or pH referencing (DSS, imidizole) were excluded from the analysis, unless otherwise noted. Numbering corresponds to assignments in Table 2.
information for each sample (Nj, where j is the jth sample), i.e., the contribution of all resonances from the original spectral region. Success of subsequent multivariate pattern recognition analysis to elucidate metabolite patterns, biomarkers, or both is critically dependent on the assumption that the input variables (K) contain the same latent information across all samples (N). In practice however, the sensitivity to sample conditions, coupled with baseline and other spectral distortions that can arise during an experiment, means that the bin integrations will not necessarily reflect true changes in the spectral areas.9,10 Since pattern recognition techniques, such as principal component analysis (PCA), depend on linear combinations of input spectral bins, artifacts in the integrations of input bins will compromise the analysis. If the metabolites of interest are modulated in a subtle manner, peak position and line width differences between samples as well as instrumental and spectral artifacts may mask significant (7) ACD/Labs. Ref Type: Internet Communication. 2006. (8) Eads, C. D.; Furnish, C. M.; Noda, I.; Juhlin, K. D.; Cooper, D. A.; Morrall, S. W. Anal. Chem. 2004, 76, 1982-1990.
changes to their concentrations, especially for low-concentration metabolites, which may be important in disease processes. Statistical tools such as orthogonal signal correction can be applied to regression-type experiments, and these have been shown to eliminate components of the data that are not relevant to the analysis, such as the effects of physiological variation or instrumental effects.11 However, these analysis corrections are still dependent on the assumption that the latent information in each variable K is the same across N samples. Several studies have examined methodologies to align spectra prior to analysis using all points in the spectra or binning.12,13 Successful implementation of these approaches would ensure the same chemical information content for each point or bin in the spectrum. One significant limitation to these approaches is that overlapping peaks are highly problematic, especially when the peak shifts of the overlapping entities are influenced by heterogeneous factors. Alignment is usually successful for the most intense peaks and less so for lower concentration metabolites. Despite these advances to spectral preprocessing algorithms and pattern recognition methods for spectral binning data, little information is available with respect to individual metabolites and their concentrations in the biofluid. Underlying any statistical treatment of NMR spectra in metabolomics is the basic notion that metabolites are the actual variables of interest. The ideal statistical treatment would be based directly on the concentrations of all metabolites in the samples, since these data represent the underlying physical model that generated the observed data (the NMR spectrum). Such treatment would allow for analysis of compounds as an entity or selectively (such as pathway modeling). In this case, there would still be N samples; however, each variable K would be the concentration of a particular metabolite. An approach to data reduction, which addresses the quantitative concern, entails profiling a spectrum by comparison to NMR spectral signatures of individual metabolites found in a reference database. This technique works by reducing spectral data to quantified metabolites, which can then be used as input variables into pattern recognition tools such as PCA or partial least squares (PLS). An advantage of this approach is a reduction in the dimensionality of problem space compared to spectral binning, as assignment of all protons in a chemical compound will show (9) Crockford, D. J.; Keun, H. C.; Smith, L. M.; Holmes, E.; Nicholson, J. K. Anal. Chem. 2005, 77, 4556-4562. (10) Potts, B. C.; Deese, A. J.; Stevens, G. J.; Reily, M. D.; Robertson, D. G.; Theiss, J. J. Pharm. Biomed. Anal. 2001, 26, 463-476. (11) Beckwith-Hall, B. M.; Brindle, J. T.; Barton, R. H.; Coen, M.; Holmes, E.; Nicholson, J. K.; Antti, H. Analyst 2002, 127, 1283-1288. (12) Stoyanova, R.; Nicholls, A. W.; Nicholson, J. K.; Lindon, J. C.; Brown, T. R. J. Magn Reson. 2004, 170, 329-335. (13) Forshed, J.; Torgrip, R. J.; Aberg, K. M.; Karlberg, B.; Lindberg, J.; Jacobsson, S. P. J. Pharm. Biomed. Anal. 2005, 38, 824-832.
Analytical Chemistry, Vol. 78, No. 13, July 1, 2006
4431
all spectral regions correlated to a specific compound. As a result, a variety of approaches to targeted profiling have recently been developed for both in vivo and ex vivo NMR.9,14,15 This methodology is most useful when the compounds of interest are previously defined and can be targeted. NMR relaxation properties of metabolites are a key factor influencing the accuracy of metabolite quantification by NMR in targeted profiling approaches.16 In high-throughput metabolomics experiments, it is not feasible to allow complete relaxation of all nuclei as experimental times would be too long due to the relatively large longitudinal relaxation time (T1) of small molecules. In principle, if sample NMR spectra are acquired under the exact same conditions (solvent suppression, relaxation time, etc.) as a particular reference database, quantitation should be ideal as each one will undergo the same relaxation phenomenon. However, in practice, this is often not feasible. For example, sample data that were acquired prior to knowledge of reference database acquisition parameters, or sample data requiring Carr-Purcell-Meiboom-Gill (CPMG)-type experiments for the reduction of interfering protein or lipid baseline artifacts,17 will produce quantification errors. Here we present a novel method of quantitatively characterizing 1H NMR spectra, hereafter referred to as “targeted profiling”. In this approach, metabolites of interest are first chemically modeled using their peak center and J-coupling information. This information is stored in a database, which is accessed during the analysis of an unknown metabolite mixture spectrum, to create a mathematical model of each metabolite in a cumulative manner. Quantification is achieved through the use of an internal standard. We have previously successfully used this approach to identify and quantify overlapped glutamine/glutamate peaks from rat brain extracts.18 Here we demonstrate the utility of the method using multivariate statistical analysis in a mixture of metabolites at physiological concentrations. We compare this approach to the ideal case of manually created bins in a spectral binning approach. Specifically, we address the effect of water suppression techniques and acquisition parameters on pattern recognition using either spectral binning or targeted profiling to determine the stability of PCA models. The effect of presaturation,19 1D nuclear Overhauser enhancement spectroscopy (NOESY)-presaturation,20 and water suppression enhanced through T1 effects (WET)21 methods of solvent suppression has previously been studied with respect to spectral binning.10 Here we additionally examine the CPMG 22 pulse sequence often used in metabolic profiling of samples containing larger MW components. We expect that the different pulse sequences will cause differences in baseline due to the difference in water suppression techniques. However, these differences will not be quantifiable for a particular pulse sequence since spectral distortions are also a result of sample conditions and magnet shimming. The choice of scaling factor, autoscaling (14) Provencher, S. W. NMR Biomed. 2001, 14, 260-264. (15) Bamforth, F. J.; Dorian, V.; Vallance, H.; Wishart, D. S. J. Inherited Metab. Dis. 1999, 22, 297-301. (16) Evilia, R. F. Anal. Lett. 2001, 34, 2227-2236. (17) Lucas, L. H.; Larive, C. K.; Wilkinson, P. S.; Huhn, S. J. Pharm. Biomed. Anal. 2005, 39, 156-163. (18) McGrath, B. M.; Greenshaw, A. J.; McKay, R.; Weljie, A. M. Slupsky C. M.; Silverstone, P. H. Int. J. Neurosci. 2006, (In Press). (19) Hoult, D. I. J. Magn. Reson. 1976, 21, 337-347. (20) Kumar, A.; Ernst, R. R.; Wuthrich, K. Biochem. Biophys. Res. Commun. 1980, 95, 1-6.
4432
Analytical Chemistry, Vol. 78, No. 13, July 1, 2006
(unit variance), or Pareto scaling is also evaluated as an influencing variable on the PCA analysis; scaling has been shown to have a greater impact on analysis than bin size for spectral binning.23 The effect of variable acquisition times is considered for the evaluation of the influence of resolution and relaxation on pattern recognition. Furthermore, the influence of water suppression techniques and relaxation is investigated as a factor in the accuracy of quantification of metabolites using the targeted profiling approach. Finally, the results of the analyses are described with respect to the impact of each pulse sequence on different chemical groups and regions of the spectrum, and T1 properties of the metabolites characterized. METHODOLOGY Sample Preparation. To unambiguously analyze the spectral information, a simplified synthetic urine for six samples was prepared using the concentrations of metabolites as described in Table 1. Metabolite concentrations were chosen to represent actual concentrations that may be observed in urine. To determine representative metabolite concentrations, we averaged the measured metabolite concentrations from the urine of 6 healthy males over 25 days of 1-methylhistidine, citrate, creatinine, fucose, glycine, glycolate, hippurate, isoleucine, lactate, and taurine, using the Chenomx NMR Suite 3.1 software (Marrie et al., unpublished observations). All assignments were confirmed by comparison with pure reference compounds. All samples contained a mixture of salts,24 including 78 mM NaCl, 16 mM Na2SO4, 18.6 mM NH4Cl, 21.2 mM KCl, 0.02% NaN3, 3.26 mM CaCl2, 3.17 mM MgCl2, 20 mM phosphate, 10 mM creatinine, 400 mM urea, 10 mM imidizole (pH indicator), 10% D2O, and 0.50 mM DSS (chemical shift indicator) and were at pH 7.00 ( 0.05 uncorrected for HOD using an appropriately calibrated pH meter. NMR Data Acquisition and Processing. The NMR spectra for each sample were acquired using four water suppression methods: NOESY-presaturation (noesypr), presaturation (presat), CPMG, and WET on a four-channel Varian Inova 600 spectrometer with a Triax-gradient 5-mm HCN probe. Acquisition parameters were as follows: presat and noesypr water saturation of 0.99 s during prescan delay, preceded by a short delay of 0.01 s; noesypr mixing time of 100 ms with prescan delay of 1.00 s; CPMG saturation time of 80 ms with prescan delay of 1.00 s. All spectra were acquired with a 12 ppm sweep width, 4-s acquisition time, four dummy scans, and 32 transients. An additional set of four spectra for each sample were acquired with the noesypr pulse sequence using a sweep width of 20 ppm and acquisition times of 2, 3, 4, and 5 s, for a total of eight spectra per sample. All spectra were zero filled to 128k data points, Fourier transformed with a 0.5-Hz line broadening applied, and manually phased and baseline corrected using VNMR software. T1 Measurements. A separate sample containing the metabolites listed in Table 1 at 1.0 mM was prepared to measure T1 values for nonoverlapping resonances. Nine spectra were acquired (21) Smallcombe, S. H.; Patt, S. L.; Keifer, P. A. J. Magn. Reson. A 1995, 117, 295-303. (22) Meiboom, S.; Gill, D. Rev. Sci. Instrum. 1958, 29, 688-691. (23) Webb-Robertson, B. J.; Lowry, D. F.; Jarman, K. H.; Harbo, S. J.; Meng, Q. R.; Fuciarelli, A. F.; Pounds, J. G.; Lee, K. M. J. Pharm. Biomed. Anal. 2005, 39, 830-836. (24) Jones, D. S.; Djokic, J.; Gorman, S. P. J. Biomed. Mater. Res., B: Appl. Biomater. 2006, 76, 1-7.
with T1 relaxation delays of 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, and 128 s and a prescan delay of 180 s. The area for individual resonances was used and fit to the equation
Table 2. T1 Values as a Function of Compound and Individual Resonance compound
I ) Io(1 - e-TR/T1) where TR is the repetition time for one scan. Targeted Profiling. Quantification was achieved using the 600MHz library from Chenomx NMR Suite 3.1 (Chenomx Inc., Edmonton, Canada), which uses the concentration of a known reference signal (in this case DSS) to determine the concentration of individual compounds. The library is predicated on a database of individual metabolite spectra acquired using the noesypr sequence with a 4-s acquisition time, and 1-s recycle delay between transients, and contained 215 metabolites. Each reference compound was fit to record peak centers and homonuclear J-coupling constants at pH 7.00, and this information was stored in a fielddependent database. Profiling of the analyte mixture spectra was accomplished using the Profiler module. Essentially, a Lorentzian peak shape model of each reference compound is generated from the database information and superimposed upon the actual spectrum. The linear combination of all modeled metabolites gives rise to the total spectral fit, which can be evaluated with a summation line. While automated numerical optimization algorithms would be possible for the data in this study, and is the long-term goal of such work, quantification was achieved by ensuring that all nonexchangeable peaks from the reference model and the analyte spectrum were at the same level in the y-dimension (concentration) by visual inspection, as this is currently the best method for analyzing biofluids. Total fitting time was ∼2 min/ spectrum for 14 compounds (metabolites in Table 1, DSS, imidazole, and urea). Spectral Binning. Integral bins were created in such a manner as to ensure that each resonance was in the same bin throughout all spectra. Custom bin sizes were created for each resonance over all spectra, with overlapping resonances being considered together in one bin. The spectra were identical to those used for quantification using the targeted profiling approach. All bins were normalized to the area of the DSS methyl peak to provide a measure of the absolute contributions of particular resonances to the spectrum, as well as the total spectral area for a comparison corrected for dilution effects. Principal Component Analysis. PCA was performed using standard procedures as implemented in Simca P+ 10.5 (Umetrics, Umeå, Sweden). Input variables consisted of either compound concentrations (targeted profiling) or integral bin areas (spectral binning). Data were preprocessed by mean-centering and either autoscaling or Pareto scaling prior to analysis. Autoscaling, also known as unit variance scaling, divides each of the “x” variables by its standard deviation, while Pareto scaling entails a division by the square root of the standard deviation. The quality of the models was judged by the goodness-of-fit parameter (R2), and the goodness of prediction parameter based on the fraction correctly predicted in a 1/7th cross-validation (Q2),25 with each sample being used in the test set. In the case of targeted profiling data, PCA models were built for quantitative analysis with the input variables (K) being the (25) Multi- and Megavariate Data Analysis: Principles and Applications; Umetrics AB: Umeå, Sweden, 2001.
taurine taurne/fucose creatinine/creatine creatinine hippurate/methylhistidine hippurate
fucose fucose/isoleucine
fucose/unknown fucose/water citrate isoleucine
lactate 1-methylhistidine
glycine glycolate creatine DSS DSS/creatinine satellite imidazole a
resonancea 1a,2b 1b,2a H6 H5 H9 H3,H5 H4 H6,H2 H8 aH6 bH6 bH3 bH4,aH2,bH5,aH4 aH3 aH5 bH1 aH1 2a,4a 2b,4b H5 H6 4b H3 H2 H3 H2 3b 3a H6 H5 H4 H2 H2 H2 H1 H3 H4 H5
T1 (s)
bin index
bin width (ppm)
2.47
1 2 3 4 6
0.035 0.095 0.053 0.018 0.050
7 8 9 10 11 12 14 15 16 17 18 19 20 21 22 23 25 26 27 28 29 30 31 32 34 35 36 37 39 40 41 42 43
0.058 0.053 0.055 0.248 0.029 0.084 0.037 0.087 0.038 0.050 0.036 0.032 0.102 0.078 0.048 0.058 0.102 0.091 0.028 0.052 0.060 0.057 0.055 0.038 0.123 0.156 0.037
2.73
2.87 3.84 2.73 0.89 1.46 1.38
1.98 0.71 0.61 1.55 1.07 1.23 1.80 1.82 1.79 3.99 0.88 0.79 1.92 3.42 3.76 2.74 3.21 1.28 2.87 1.62 1.65 7.46 7.71
0.016 0.033 0.058 0.078 0.060
44 45
Note that numbering corresponds to structures in Figure 1.
percent error in concentration determination of each metabolite over all samples (N). All variables used for quantitative analysis were Pareto scaled. Separate models were built to examine the effect of pulse sequence and relaxation times. RESULTS Spectral Binning and T1 Analysis. A total number of 45 custom bins were created from the data representing the 11 compounds listed in Table 1 and depicted in Figure 1. Each bin was custom-sized considering all spectra as a group, such that peak and peak clusters would remain in the same bin over all spectra. The goal of the binning procedure was to ensure that peak areas from each resonance were assigned to a specific bin (and not divided across two or more bins). This was done to account for the minor change in pH and ionic strength between samples. This method of binning is simply a manually optimized implementation of variable bucketing algorithms 7,13 where local minima in spectra are recognized and buckets are adjusted to prevent peaks from being split. Typically, the buckets range from 0.02 to 0.06 ppm using this method. In our study, the bin sizes ranged from 0.016 (creatine peak at ∼3.91 ppm) to 0.248 ppm (exchanging hippurate peak at ∼8.40 ppm) (Table 2). As can be seen from Table 2, most buckets were less than 0.1 ppm in width. Analytical Chemistry, Vol. 78, No. 13, July 1, 2006
4433
Figure 2. Representative spectral bins showing the superposition of all spectra (6 samples, 48 spectra) collected. (A) Downfield 1-methylhistidine aromatic proton (bin 35) illustrates the large variation in chemical shift with slight changes in pH and ionic strength between samples. In this case, the bin width was set to 0.156 ppm to account for all variations. (B) Bin 6 contains resonances from 1-methylhistidine, hippurate, and glycolate. Hippurate and glycolate show little chemical shift movement whereas the resonance from 1-methylhistidine changes. In this case, the bin width was set to be much smaller (0.05 ppm) since the 1-methylhistidine resonance in this bin does not move as much with slight variations in pH and ionic strength as the resonance in bin 35. (C) Targeted profiling illustrates how overlapping resonances may be deconvoluted. The spectral line is shown in black, 1-methylhistidine in purple, hippurate in blue, glycolate in green, and summation of all three fit compounds in red. Note how the summation (red) line and black line superimpose well.
Interestingly, with very similar ionic strengths and pH’s within 0.1 pH unit, many peak chemical shifts change dramatically. This could unfortunately lead to a great deal of error when doing spectral PCA type of analysis unless these changes are taken into account. Figure 2A shows an example of a bin (bin 35) corresponding to the downfield aromatic proton from 1-methylhistidine. Shown is the superposition of all spectra (six different samples) containing this resonance. As may be seen, the bin width is relatively large due to large chemical shift changes of this resonance with the minor pH and ionic differences between the samples. Resonance overlap between and within compounds resulted in some bins consisting of more than one resonance and, in some cases, resulted in multiple resonances per bin. Figure 2B depicts an overlapping region (bin 6) containing 1-methylhistidine, glycolate, and hippurate resonances with assignment of peaks using spectral profiling software shown in Figure 2C. The peaks from glycolate (∼3.945 ppm) and hippurate (doublet at ∼3.959) do not shift significantly between samples; however, the resonance from 1-methylhistidine (doublet of doublets centered at ∼3.953 ppm) shows significant chemical shift changes. For T1 analysis, only bins with unique resonances originating from a single compound were analyzed. Table 2 lists the T1 values obtained for resonances that could be separated into unique bins. Citrate clearly demonstrates the fastest relaxation in the group of compounds considered with T1 values of 0.71 and 0.61 s for the two bins considered. The two bins containing imidazole had 4434
Analytical Chemistry, Vol. 78, No. 13, July 1, 2006
the longest T1 values of 7.46 and 7.71 s, respectively. The methyl group of DSS had a T1 value of 2.87 s. Bins containing peaks from DSS or imidazole were excluded from the principal component analysis procedure. PCA Analysis. PCA is a classical multivariate statistical method that has been widely used in all types of data analysis to reduce the dimensions of a single set of measured variables to determine which variables or set(s) of variables form coherent subsets that are relatively independent of one another. In other words, PCA transforms a number of correlated variables into a smaller number of uncorrelated variables called principal components. The first component accounts for the greatest variability in the data, and each succeeding component accounts for that much of the remaining variability. Often the first two principal components are represented in a PCA plot, and these represent the two largest components separating the data. However, PCA data may be best separated using other components. PCA has been widely used in the field of metabolomics (for a review, see refs 2 and 3) for comparing raw NMR spectral data obtained from various biofluids and tissue extracts to determine whether NMR spectral differences exist between groups such as “healthy” and “diseased”. Using NMR spectral data for this type of analysis can be tricky at best as 1H NMR chemical shifts are sensitive to pH and ionic strength as shown above. Often bins are set (rather than using discrete points) since peaks move around, but as shown, these bins need to be dynamic in nature. Even the variable-sized bins approach, however, has its limitations
Figure 3. PCA, using Pareto scaling, of the NMR spectral data arising from six different samples, with metabolites at the concentrations indicated in Table 1. 1-D NMR spectra were acquired using different water suppression schemes present in four pulse sequences: WET, prnoesy, presat, and CPMG. Additionally, four different acquisition times (2, 3, 4, and 5 s) with a constant 1-s relaxation delay were collected using the prnoesy pulse sequence. Each sample is represented by a unique color (samples: 1, black; 2, red; 3, blue; 4, green; 5, orange; 6, magenta), and selected pulse sequences have been indicated (2, CPMG; *, WET; 4, PRESAT; 9, PRNOESY; ×, PRNOESY-2s; 0, PRNOESY3s; b, PRNOESY-4s; O, PRNOESY-5s). The ellipses in all PCA plots represent the 95% confidence limit (Hotelling T2). (A) Spectral binning data normalized to DSS. (B) Spectral binning data normalized to total area (Table 3). (C) Targeted profiling data.
as the complexity of the spectrum increases. Furthermore, the use of bins complicates analysis as more than one bin can represent the same compound. In our case, we have broken up the spectra into 45 bins. Analyzing compounds with their corresponding concentrations directly (targeted profiling) greatly simplifies the analysis and reduces our 45 bins to 13 compounds making analysis easier to interpret. We are interested in a direct comparison of these two techniques of data reduction by comparing the use of different data acquisition schemes (pulse sequences and acquisition times). Figure 3 shows PCA plots using our variable spectral binning and targeted profiling approaches separating samples by concentration, pulse sequence, and spectral acquisition time used. The
number of components (A) and explained variance relating to the goodness of fit (R2) for each model are provided in Table 3. Each color in Figure 3 represents a different sample with varying concentrations of metabolites (as described in the figure legend), and each point represents a different acquisition condition (either pulse sequence or acquisition time). In Figure 3A, the spectral data have been normalized to DSS, whereas in Figure 3B, the spectral data have been normalized to the total area under the spectrum. Figure 3C illustrates the PCA plot using direct compound profiling. As may be observed, irrespective of the acquisition time and method of analysis, each set of data collected with the same pulse sequence tends to cluster together when using the targeted profiling method. Interestingly, the binned data Analytical Chemistry, Vol. 78, No. 13, July 1, 2006
4435
Table 3. Table of Models Pareto scaling model
A
binned, normalized to DSS binned, normalized to spectral area targeted profiling targeted profiling errors (pulse sequence) targeted profiling errors (acquisition time)
12 12 5 2 2
R2
Q2
autoscaling A
R2
Q2
0.996 0.881 5 0.776 0.517 0.995 0.839 6 0.841 0.492 0.993 0.853 6 0.962 0.768 0.606 -0.108 0.825 0.607
normalized to DSS are more susceptible to changes in acquisition time than the binned data normalized to total area. In addition, a comparison of the binned data normalized to the total area with the targeted profiling data was made. In addition
to the goodness-of-fit metric (R2), the fraction of data correctly predicted (Q2) using 1/7 cross-validation (Umetrics AB 2001) was determined. A value of 1 for both of these variables implies complete data correlation and 100% prediction for a given data set. Typically, Q2 is lower than R2; however, a significant difference between the two parameters is undesirable, and a negative Q2 implies that the model is not predictive. Figure 4A shows a plot of the contribution of each bin to the overall model, and Figure 4B shows each compound for the targeted profiling model. The R2 and Q2 values for all variables in the spectral binning model are relatively good except for bins 22-26, which correspond to isoleucine. In other words, all bins fit the model relatively well and are reasonably well predicted using cross-validation except
Figure 4. Comparison of contribution plots between (A) spectral binning (over 12 components separating the data in Figure 2B, Table 3) and (B) targeted profiling (over 5 components separating the data in Figure 2C). Bin number (for spectral binning) or metabolite ID (for targeted profiling, Table 3) are plotted as a function of R2 (gray) and Q2 (black). R2 represents the goodness of fit, and Q2 represents the fraction of data correctly predicted in a 1/7 cross-validation and provides an estimate of the predictive ability of the model. Negative Q2 indicates a nonpredictive model for bins 22-26. 4436 Analytical Chemistry, Vol. 78, No. 13, July 1, 2006
for the bins representing this compound. In addition, bins 10, 12, 14, 29, and 36 show a higher than desirable gap between R2 and Q2, prompting further investigation. Bin 10 derives from a slow exchanging peak from hippurate near 8.4 ppm, bin 12 from an overlap of fucose and isoleucine near 1.2 ppm, bin 14 a fucose resonance at 3.61, bin 29 from the lactate peak near 4.05 ppm, and bin 36 the glycine resonance near 3.5 ppm. These latter three bins suggest some residual artifacts from solvent suppression which the Pareto scaling was not able to completely eliminate. This phenomenon has been observed by other researchers with respect to differences in water suppression methods.10 It is worth noting that compound-based profiling does not produce the same result for isoleucine (Figure 4B), where it contributes positively to Q2. In addition, none of the other compounds appears affected by solvent suppression. This result implies that it is not the integral values from the compound itself causing the deviation in the binning model, but some other artifact. Further analysis of the bins comprising peaks related to isoleucine demonstrates that the combination of low concentration and high J-multiplicity of the aliphatic resonances produces signal very close to the baseline. It is the baseline that may be the source of error in this case. In the targeted profiling model, since all of the regions in isoleucine combine to form a single variable, these variations are eliminated. Effect of Scaling Parameters on PCA Analysis. A number of different scaling parameters are available for preprocessing data prior to PCA analysis. The choice of scaling parameter is important as it defines the relationship between variables. Autoscaling forces all x values to have equal weight, irrespective of the starting intensity, and hence is prone to distortion from poor baseline and other spectral artifacts.10 On the other hand, Pareto scaling gives a greater weight to variables that start with a larger valuesso, for example, intense peaks in an NMR spectrum will have greater weight than weak ones. Figure 5A shows the pattern that emerges when binned data normalized to the DSS methyl peak are subject to PCA using autoscaling. There is a separation of the six samples through PC1. However, in PC2, data separation is caused by a factor other than sample difference. Inspection of the outliers indicates that acquisition with CPMG, and to a lesser degree WET pulse sequences, was responsible for the PC2 separation. The same analysis with Pareto scaling instead of autoscaling (Figure 3A) demonstrates better clustering, although the effect of the CPMG pulse sequence is still evident. While the separation is improved in PC2 with Pareto scaling, the dispersion through PC1 and PC2 still suggests that a combination of intrinsic sample differences and intrasample acquisition condition variations are responsible for diffuse clustering. The bins resulting from the integration procedure were further analyzed, after normalization to the entire total spectral area, using both autoscaling (Figure 5B) and Pareto scaling (Figure 3B). In this instance, autoscaling still provides components with mixed inter- and intrasample variation. The most obvious intersample separation of all binned data, reducing the differences to just acquisition parameters in PC2, is with the total area normalization using Pareto scaling. Figure 6A,B illustrates loading plots for the Pareto scaled and autoscaled data normalized to total spectral area (Figures 3B and 5B, respectively). The compounds causing variation in the PCA plots are highlighted in Figure 6A, although
in some cases (e.g., bins 2, 3, 6, 30, and 31), there are overlapping resonances possibly confounding the analysis. There is a substantial change in the loadings plot due to scaling. In the case of the autoscaling, bins from fucose, lactate, creatine, and isoleucine between 3.65 and 4.21 ppm are distinct from the remaining bins (Figure 6B), suggesting that the acquisition parameters (specifically the CPMG pulse sequences and, to a lesser extent, WET) cause a significant difference in this region (Figure 7). Pareto scaling (Figure 6A) reduces this effect and allows resonances from individual compounds to be discriminated with relative ease as they cluster into moderately defined regions. Compound-based targeted profiling results in transforming integral regions from the various resonances of a single compound into a single parameter, its concentration (Figure 2C). One would therefore expect PCA results from targeted profiling to be similar to those from spectral binning. This result is shown in Figures 5C (Autoscaling) and 3C (Pareto scaling). As with the binned data, Pareto scaling of targeted profiling data results in tight clustering. Interestingly, the effect of scaling is marginal, only causing a slight change in both PC1 and PC2. In both cases, the effect of acquisition parameters is clearly minor relative to the intersample variation. This result is reflected in the loadings plots (Figure 6C and D, for the Pareto scaled and autoscaled models, respectively). Since the initial input to PCA analysis was compound concentrations, the loadings plot is relatively easy to interpret. A change in scaling causes the variation in both components to be interpreted slightly differently. However, key features such as the difference between citrate and hippurate in PC2 are preserved. This is in contrast to spectral binning (Figure 6B) where autoscaling indicated larger differences in PC2 due to pulse sequence selection rather than variation between citrate and hippurate, which appears to be obscured by artifacts. Effect of Acquisition Parameters. One distinct advantage of using a targeted profiling approach is that the data of interest are quantitative. As a result, we were interested in assessing the accuracy of quantitation and the underlying factors from this study that might influence absolute concentration determination. The effect of the various acquisition parameters on the analysis was evaluated using the difference between the measured concentrations from the targeted profiling approach with the actual known concentrations of metabolites. This method was chosen since the difference between the measured and actual concentrations represents a vector in multivariate space that provides valuable information about how the particular condition (acquisition time, pulse sequence) will change the measured concentration. In essence, this is a pattern recognition analysis on the error associated with concentration determination; larger vectors will be associated with increased error. Because this is a measure of error in quantitation, all metabolites, including imidazole, were included in the analysis. Typical errors for the prnoesy and presat sequences were
