Application of a Bayesian Deconvolution Approach for High

May 6, 2010 - ... the Metabolic Effects of Acute Phenobarbital Exposure in Liver Tissue .... A review of blind source separation in NMR spectroscopy. ...
0 downloads 0 Views 273KB Size
Download PDF
Anal. Chem. 2010, 82, 4479–4485

Application of a Bayesian Deconvolution Approach for High-Resolution 1H NMR Spectra to Assessing the Metabolic Effects of Acute Phenobarbital Exposure in Liver Tissue Denis V. Rubtsov,† Claire Waterman,† Richard A. Currie,‡ Catherine Waterfield,‡ Jose´ Domingo Salazar,‡ Jayne Wright,‡ and Julian L. Griffin*,†,§ The Department of Biochemistry and the Cambridge Systems Biology Centre, Tennis Court Road, University of Cambridge, Cambridge, CB2 1QW, U.K., and Human Safety, Syngenta, Jealott’s Hill International Research Centre, Bracknell, Berkshire RG42 6EY, U.K. High-resolution 1H NMR spectroscopy is frequently used in the field of metabolomics to assess the metabolites found in biofluids or tissue extracts to define a metabolic profile that describes a given biological process. In this study, we aimed to increase the utility of NMR-based metabolomics by using advanced Bayesian modeling of the time-domain high-resolution 1D NMR free induction decay (FID). The improvement over traditional nonparametric binning is twofold and associated with enhanced resolution of the analysis and automation of the signal processing stage. The automation is achieved by using a Bayesian formalism for all parameters of the model including the number of components. The approach is illustrated with a study of early markers of acute exposure to different doses of a well-characterized nongenotoxic hepatocarcinogen, phenobarbital, in rats. The results demonstrate that Bayesian deconvolution produces a better model for the NMR spectra that allows the identification of subtle changes in metabolic concentrations and a decrease in the expected false discovery rate compared with approaches based on “binning”. These properties suggest that Bayesian deconvolution could facilitate the biomarker discovery process and improve information extraction from high-resolution NMR spectra. Metabolomics aims to understand biological processes by examining the global metabolic changes associated with them. It has found applications in drug toxicology,1,2 understanding gene * To whom correspondence should be addressed. Phone: +44 1223 764922. Fax: +44 1223 766022. E-mail: [email protected]. † Department of Biochemistry, University of Cambridge. ‡ Syngenta. § Cambridge Systems Biology Centre, University of Cambridge. (1) Lindon, J. C.; Nicholson, J. K.; Holmes, E.; Antti, H.; Bollard, M. E.; Keun, H.; Beckonert, O.; Ebbels, T. M.; Reilly, 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. (2) Nicholson, J. K.; Connelly, J.; Lindon, J. C.; Holmes, E. Nat. Rev. Drug Discovery 2002, 1, 153–161. 10.1021/ac100344m  2010 American Chemical Society Published on Web 05/06/2010

function,3 and monitoring disease processes.4 High-resolution 1H NMR spectroscopy is one technique used widely to assess the metabolites found in biofluids or tissue extracts of such metabolomic studies.5 In most NMR-based metabolomic experiments, the signals from metabolites are normally present as a mixture of overlapping resonances, making quantification difficult. However, deconvolution of an NMR spectrum into its individual components has proven to be a nontrivial task. Commercially available software resorts to peak fitting, which often perform badly for complex resonances and broad signals (as typical in spectra of blood plasma, for example). While 2D homonuclear and heteronuclear NMR investigations will always remain essential to de novo structure elucidation of unassigned NMR resonances, current requirements in terms of time for performing most 2D experiments preclude their routine use in most metabolomic studies. Furthermore, at the moment, extracting relevant information from NMR data sets can be a laborious manual process that affects reproducibility and introduces additional uncertainty and errors into downstream analysis. Another challenge of modern NMR-based metabolomics is associated with the need for increasingly high-throughput and automated sample scanning that leads to a dramatic increase in the scale of metabolomic experiments, now routinely dealing with hundreds of biological samples as a part of epidemiological-scale studies, for example.6,7 As a single high-resolution NMR spectrum reports tens of thousands of variables, the automation of information extraction becomes imperative in order to retain the benefits of the high-throughput NMR analysis. At the same time it is important to ensure the highest possible quality of data processing (3) Griffin, J. L. Philos. Trans. R. Soc., B 2004, 359, 857–871. (4) Kirschenlohr, H. L.; Griffin, J. L.; Clarke, S. C.; Rhydwen, R.; Grace, A. A.; Schofield, P. M.; Brindle, K. M.; Metcalfe, J. C. Nat. Med. 2006, 12, 862– 862. (5) Robertson, D. G.; Reily, M. D.; Sigler, R. E.; Wells, D. F.; Paterson, D. A.; Braden, T. K. Toxicol. Sci. 2000, 57, 326–337. (6) Dumas, M. E.; Maibaum, E. C.; Teague, C.; Ueshima, H.; Zhou, B.; Lindon, J. C.; Nicholson, J. K.; Stamler, J.; Elliott, P.; Chan, Q.; Holmes, E. Anal. Chem. 2006, 78, 2199–2208. (7) Smith, L. M.; Maher, A. D.; Want, E. J.; Elliott, P.; Stamler, J.; Hawkes, G. E.; Holmes, E.; Lindon, J. C.; Nicholson, J. K. Anal. Chem. 2009, 81, 4847–4856.

Analytical Chemistry, Vol. 82, No. 11, June 1, 2010

4479

so that the gains from using the new and powerful equipment is not canceled out by crude and ad-hoc processing methods. Currently one of the most popular and relatively time-efficient approaches for NMR data processing relies on interval-averaging or binning. Interval-averaging is usually performed prior to statistical analysis to limit the noise added by small pH and ionic strength-dependent variations in peak shift and decrease the number of variables.8,9 Binning is a nonparametric approach and is easy to implement but it has several disadvantages. Lowintensity peaks can be obscured by noise variation included in the same bin that limits the sensitivity of the approach, while the resolution is restricted when several overlapping peaks are summarized by a single bin making it impossible to distinguish between their individual variations. As the bins do not generally correspond to individual resonances, the process of identification invariably includes confirmation of the apparent changes by referring back to the original spectra. Further complications can arise as the inclusion of bins containing merely noise can prevent the building of successful multivariate models as described in ref 10. A number of algorithms have been developed to address some of these concerns by locating bin edges in the local minima of the spectra, for example.11,12 Other focus on peak alignment13 and use curve-fitting14 and targeted profiling.15 They generally perform better than equidistant binning as the variables generated by these methods represent the NMR spectra more accurately. Principal components analysis (PCA) has also been used to quantify peak areas.16,17 The methods described above are nonparametric and therefore ignore a significant amount of prior knowledge about the nature of the NMR signal and, thus, are not optimal for obtaining certain peak attributes, such as chemical shift or line width. It is recognized that simple nonparametric methods have limited accuracy when dealing with spectral regions densely populated by almost coresonant or overlapping peaks that commonly occur in certain types of high-resolution 1H NMR spectra.17-19 A number of advanced approaches have been developed that use model-based quantification (see ref 19 for a recent overview). One class of algorithms includes popular methods such as linear (8) Beckwith-Hall, B. M.; Nicholson, J. K.; Nicholls, A. W.; Foxall, P. J.; Lindon, J. C.; Connor, S. C.; Abdi, M.; Connelly, J.; Holmes, E. Chem. Res. Toxicol. 1998, 11, 260–272. (9) Holmes, E.; Nicholls, A. W.; Lindon, J. C.; Ramos, S.; Spraul, M.; Neidig, P.; Connor, S. C.; Connelly, J.; Damment, S. J.; Haselden, J.; Nicholson, J. K. NMR Biomed. 1998, 11, 235–244. (10) Halouska, S.; Powers, R. J. Magn. Reson. 2006, 178, 88–95. (11) Davis, R. A.; Charlton, A. J.; Godward, J.; Jones, S. A.; Harrison, M.; Wilson, J. C. Chemom. Intell. Lab. Syst. 2007, 85, 144–154. (12) De Meyer, T.; Sinnaeve, D.; Van Gasse, B.; Tsiporkova, E.; Rietzschel, E. R.; De Buyzere, M. L.; Gillebert, T. C.; Bekaert, S.; Martins, J. C.; Van Criekinge, W. Anal. Chem. 2008, 80, 3783–3790. (13) Witjes, H.; Melssen, W. J.; Zandt, H. J. A. I.; van der Graaf, M.; Heerschap, A.; Buydens, L. M. C. J. Magn. Reson. 2000, 144, 35–44. (14) Crockford, D. J.; Keun, H. C.; Smith, L. M.; Holmes, E.; Nicholson, J. K. Anal. Chem. 2005, 77, 4556–4562. (15) Weljie, A. M.; Newton, J.; Mercier, P.; Carlson, E.; Slupsky, C. M. Anal. Chem. 2006, 78, 4430–4442. (16) Stoyanova, R.; Brown, T. R. NMR Biomed. 2001, 14, 271–277. (17) Elliott, M. A.; Walter, G. A.; Swift, A.; Vandenborne, K.; Schotland, J. C.; Leigh, J. S. Magn. Reson. Med. 1999, 41, 450–455. (18) Vanhamme, L.; van den Boogaart, A.; Van Huffel, S. J. Magn. Reson. 1997, 129, 35–43. (19) Poullet, J. B.; Sima, D. M.; Van Huffel, S. J. Magn. Reson. 2008, 195, 134– 144.

4480

Analytical Chemistry, Vol. 82, No. 11, June 1, 2010

prediction singular value decomposition (LP SVD), Hankel SVD (HSVD), the Prony method, and the matrix pencil.20-23 However, their performance can be strongly affected by the choice of the model order, which gives rise to a detection problem. Bayesian inference has also been previously applied to the problem.24-28 These advanced techniques function effectively over various experimental setups and can provide considerable improvement in the analysis of comparatively sparse 31P and 13C or in vivo 1H NMR spectra. Nevertheless, despite the considerable attention, the problem of detection and estimation remains generally unsolved for difficult cases when SNR is low and frequencies of components are closely spaced (coresonant).23 High-resolution in vitro 1H NMR spectra are particularly affected by overlapping peaks and low sensitivity. In addition, the important aspects of practical routine application of the advanced methods such as signal preprocessing and imposing constraints on the parameters are rarely addressed in an inclusive manner that limits their wider adoption. In this article, we demonstrate the practical application of advanced data processing techniques using Bayesian statistics, as described in ref 29, in a real-life metabolomic study and give a comprehensive description of a framework for the routine NMR data processing that seeks to minimize user intervention and at the same time provide a high level of information recovery. The performance of the approach is demonstrated on a toxicological study of the effects of phenobarbital (PB) exposure on the metabolic processes in the rat liver examined by high-resolution 1 H NMR spectroscopy. At the core of the framework is the Bayesian parametric modeling approach that utilizes knowledge about properties of the free induction decay signal. We demonstrate the improvement in biomarker detection that is based on improved information recovery offered by the Bayesian model as compared to traditional nonparametric approaches. In particular, the improved resolution allowed for the identification of changes in the concentration of glycine that would otherwise be undetected with a traditional binning approach. Generally the deconvoluted peaks provided a better basis for both univariate statistical testing and multivariate modeling, with a 2-fold reduction in the expected false discovery rate (FDR) and better predictive properties than binned spectra. MATERIALS AND METHODS Animal Handling and Sample Collection. This study has previously been described in ref 30. The 9-week old male Fischer (F344) rats were randomly divided into 4 groups of 20 animals. Each group was exposed to a different concentration (0, 50, 500, (20) Hoffman, R. E.; Levy, G. C. Prog. Nucl. Magn. Reson. Spectrosc. 1991, 23, 211–258. (21) Barone, P.; Guidoni, L.; Ragona, R.; Viti, V.; Furman, E.; Degani, H. J. Magn. Reson., Ser. B 1994, 105, 137–146. (22) Hua, Y. B.; Sarkar, T. K. IEEE Trans. Signal Process. 1991, 39, 892–900. (23) Vanhamme, L.; Sundin, T.; Hecke, P. V.; Huffel, S. V. NMR Biomed. 2001, 14, 233–246. (24) Bretthorst, G. L. J. Magn. Reson. 1990, 88, 571–595. (25) Bretthorst, G. L. J. Magn. Reson. 1990, 88, 552–570. (26) Bretthorst, G. L. J. Magn. Reson. 1990, 88, 533–551. (27) Dou, L. X.; Hodgson, R. J. W. Inverse Probl. 1996, 12, 121–137. (28) Dou, L.; Hodgson, R. J. W. Inverse Probl. 1995, 11, 1069–1085. (29) Rubtsov, D. V.; Griffin, J. L. J. Magn. Reson. 2007, 188, 367–379. (30) Waterman, C. L.; Currie, R. A.; Cottrell, L. A.; Dow, J.; Wright, J.; Waterfield, C. J.; Griffin, J. L. BMC Genomics 2010, 11, 9.

Table 1. Classification Errors of Orthogonal Partial Least Squares -Discriminate Analysis (OPLS-DA) Modelsa day 1 binning deconvolution

day 3

day 7

day 14

correct class

error score

correct class

error score

correct class

error score

correct class

error score

52% 68%

3000 2050

63% 63%

2600 1650

55% 60%

5350 3000

60% 85%

3450 1450

a The percentage of correctly classified test samples and the summary error score. All OPLS-DA models consisted of two components and were fitted with the leave-one-out cross-validation procedure. The percent of correctly classified samples was measured using the left-out samples. The error score was computed as a sum of differences between the true and predicted PB dose in order to reflect the higher risk associated with a bigger difference between predicted and true class labels, particularly for the highest (carcinogenic) dose. The higher correct class prediction and lower error score are associated with a better model and highlighted in bold.

or 1000 ppm) of the sodium salt of phenobarbital which was added to the standard laboratory diet and milled to homogeneity. Five animals from each group were killed by exsanguination under halothane anesthesia after 1, 3, 7, and 14 days of PB exposure. Livers were removed immediately, and a sample was frozen for subsequent analysis (∼75 mg). Tissue extraction was performed as reported in ref 30, and the aqueous fraction examined by NMR spectroscopy. High-Resolution 1H NMR Spectroscopy. Dried aqueous extracts from 500 µL of the aqueous layer were dissolved in phosphate buffered D2O (pH 7.0, 600 µL, 33 mM Na2HPO4, 6.7 mM NaH2PO4, 0.1% sodium azide, 0.17 mM TSP). 1H NMR spectra were acquired on a Bruker DPX400 spectrometer operating at 400.13 Hz using a 5 mm probe. The standard “noesypr1d” pulse sequence (Bruker Analytik GmbH, Germany) was used for solvent suppression (relaxation delay ) 2 s, t1 ) 3 µs, mixing time ) 150 ms, solvent presaturation applied during the relaxation time and the mixing time). A total of 128 transients were collected into 32K data points over a spectral width of 20 ppm at 37 °C. Bayesian Modeling of a Time-Domain NMR Signal. The Bayesian deconvolution of the NMR signal is a parametric approach that seeks to estimate both the number of harmonics and their parameters from a complex-valued time-domain NMR signal. Its detailed description can be found in ref 29. The free induction decay (FID) signal is modeled as a mixture of exponentially damped sinusoids: k

y(t) )

∑ a z (t) + ε(t) i i

(1)

i)1

where zi(t) ) ej2πωit-Rit are signal poles, ai ) Ai ej2πφi are complex amplitudes, t are time points, K is the number of model components, φ is the phase shift, R is the decay rate, ω is the frequency, and ε is the random noise. In the vector-matrix form: y ) Da + ε

(2)

where y ) [yt]T is a vector of N FID data points, D ) [Zti] is the N × K matrix of the signal poles, a )[ai]T, the vector of K complex amplitudes. The number of the sinusoids is usually not known beforehand and is estimated from data along with their parameters. The joint posterior function is highly nonlinear and cannot be resolved in closed form. The details of the Bayesian modeling process can

be found in ref 29, while the details of automated phase correction, phase alignment, normalization, and frequency domain transformations are given in the Supporting Information. Statistical Analysis and Interpretation. Peak amplitudes were normalized to unit variance. Per-peak two-tailed Wilcoxon rank sum tests31 were used for tests of difference between two classes. Kruskal-Wallis tests31 were used for one-way within time point comparisons between all dose groups. The multiple comparison problem was addressed with the FDR approach.32-34 FDRs were calculated using the Bootstrap q-value estimation algorithm of ref 35. The software package SIMCA-P+ (version 11, Umetrics, Umeå, Sweden) was used for multivariate statistical analysis as previously described in ref 30. RESULTS Predictive Models of PB Exposure. In order to examine the ability of both approaches to predict the PB dose from NMR data, a orthogonal partial least squares (OPLS) regression model was built with the PB dose as a target variable and binned data or amplitudes of the deconvoluted Lorenzian peaks as independent variables. The measure of predictive error is reported in Supplementary Table 1 in the Supporting Information. Both approaches demonstrated good predictive potential with the deconvolution having slightly better results for all the time points but day 1. A classification model based on orthogonal partial least squares-discriminate analysis (OPLS-DA) was built that classified the NMR spectra into one of four classes according to their PB exposure. Only the fourth class corresponding to the 1000 ppm dose is assumed to carry a high risk of developing liver cancer over the long-term, and thus, an approach that can distinguish this group from the others is particularly important. The deconvoluted spectra show consistently better classification results, especially when the difference in PB dose levels between the true and predicted classes is considered. The predictions also perform better if one considers that the misclassification to a significantly different PB dose class carries more risk than a misclassification to just one dose level up or down (Table 1). The examination of the coefficients of the OPLS-DA models allowed for the identification of the major changes in the data responsible for separation. The results obtained from both (31) Conover, W. J. Practical Nonparametric Statistics, 3rd ed.; Wiley: New York, 1999. (32) Benjamini, Y.; Hochberg, Y. J. R. Stat. Soc., Ser. B 1995, 57, 289–300. (33) Storey, J. D.; Tibshirani, R. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 9440– 9445. (34) Tusher, V. G.; Tibshirani, R.; Chu, G. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 5116–5121. (35) Storey, J. D. J. R. Stat. Soc., Ser. B 2002, 64, 479–498.

Analytical Chemistry, Vol. 82, No. 11, June 1, 2010

4481

Table 2. Metabolite Changes Associated with Exposure to Phenobarbitala relative increase with increasing dose of PB time point

binning

day 1

succinate glutamine glutathione

day 3

succinate glutamine glutathione taurine

day 7

succinate glutamine glutathione taurine

day 14

succinate phosphocholine

deconvolution succinate glutamine glutathione ATP/UTP glycine glutamate NAD succinate glutamine glutathione taurine ATP/UTP glycine 3-hydroxybutyrate glutamate NAD succinate glutamine glutathione glycine ATP/UTP taurine lipids succinate phosphocholine ATP/UTP glutamine glutathione glycine lipids

relative decrease with increasing dose of PB binning glucose glycogen adenosine leucine

leucine valine choline glucose glycogen adenosine

deconvolution glucose glycogen adenosine leucine/isoleucine valine uridine serine(?) leucine choline glucose glycogen adenosine myoinositol serine(?) GPC

glucose glycogen adenosine

glucose glycogen adenosine serine GPC nicotinamide

glucose glycogen adenosine nicotinamide adenosine

glucose glycogen adenosine nicotinamide adenosine malate

a Metabolite changes in the liver in response to an increasing dose of PB as detected by the binning and deconvolution processing of 1H-NMR spectra across the time course of the study. The bold font highlights the changes identified in both data sets. The question marks identify the need for further confirmation of the assignment. GPC is glycerophosphocholine.

approaches agreed very well (Table 2). Deconvoluted resonances offer a good starting point for NMR spectra interpretation using NMR standard databases or cross-correlation for peak assignment. The peaks were grouped according to similarity of their profiles across different samples as the highly correlated peaks usually belong to the same metabolic fingerprint.36 The assignment of each group of peaks was performed based on the reference spectra from BMRB (www.bmrb.wisc.edu) and HMDB (www.hmdb.ca) databases and the Chenomx software package (www.chenomx.com). The deconvolution method provided similar indications for the biggest changes associated with PB exposure compared with the integral integration approach including a dose dependent decrease in glucose and glycogen at all time points in rats exposed to PB. There was also an increase in succinate, glutamine, and glutathione and a decrease in adenosine with an increasing dose of PB at all time points. In addition, the Bayesian modeling provided identification of more subtle changes such as an increase in glycine, which is discussed in more detail below. The indicative findings of multivariate models were further tested by univariate statistics. The Kruskal-Wallis test was applied in order to identify the amount of variation associated with the phenobarbital dose, time course, and their interaction. The FDR level was set to 0.05 to control for multiple testing (see below). Significant changes associated with the increase of PB intake were (36) Holmes, E.; Cloarec, O.; Nicholson, J. K. J. Proteome Res. 2006, 5, 1313– 1320.

4482

Analytical Chemistry, Vol. 82, No. 11, June 1, 2010

found at all time points. Interestingly many metabolites demonstrated similar response to a PB dose irrespective of the time course. One of the persistent changes is a decrease in adenosine and increase in adenosine triphosphate (ATP) and uridine 5′monophosphate (UMP), which is observed at all time points. Increases in succinate, glutamine, glutathione, glutamate, glycine, 3-hydroxybutyrate, and lipids are also evident. There was also a decrease in glycogen, glucose, and isoleucine. The time-dependent changes include increases in phosphocholine at day 14. Comparison of Bayesian Deconvolution and Binning. The general examination of the quality of the information recovery suggests that the Bayesian modeling provides better representation of the variability in the NMR spectra with a median R2 ) 0.95 across all models as compared to equidistant binning, median R2 ) 0.65, even though the number of the model components (351) is significantly less than the number of bins (941). The reason for such a remarkable improvement is better resolution and sensitivity offered by the parametric model. Even for a relatively fine bin size of 0.01 ppm, there are more than 50 bins that contain two or more separate peaks, thus reducing the resolution of the signal. At the same time, the sensitivity is hampered by averaging over the noise data points as well as genuine signal contribution that causes a less accurate estimation of low-amplitude resonances. Considering a change in glycine concentration as an example of a potentially interesting biomarker in this study, it was detected with the Bayesian modeling but overlooked by the binning

Figure 1. Identification of the change in the 3.55 ppm glycine resonance.Part a shows deconvoluted peaks in the region of 3.51-3.59 ppm with the assigned glycine resonance. The edges of the bins are shown with vertical broken lines. The solid line boxes illustrate average values of the free bins covering most of the intensity of the glycine peak. Parts b, c, and d show the boxplots summarizing the intensity of the bins containing the glycine resonance across different PB dose levels. The ANOVA test shows that there is a significant (p-values of 0.0007, 0.022, and 1.3 × 10-12, respectively) decrease between the control and at least one of the treatment doses driven by the changes in glucose and glycogen. Part e shows that a modeled resonance for glycine in fact demonstrated the opposite direction of change to the trends detected in the binned data, with a p-value of 0.0004, associated with an increase between the control and the 1000 ppm PB dose. The increase in glycine was confirmed by GC/MS analysis of the aqueous phase of the liver tissue extract.

approach. The glycine peak is situated at ∼3.55 ppm in a region where glucose and glycogen resonances are dominant (Figure 1a). In the case of binning, the profiles of the three bins of the 3.54-3.57 ppm region are highly correlated across the samples treated with different doses of PB, with Pearson correlation coefficients between 0.59 and 0.93. If we examine the amplitude of the peak in the Bayesian model provisionally assigned to glycine and neighboring resonances, there is a noticeable difference in behavior with glycine increasing in the treated samples compared to the controls and both glucose and glycogen decreasing (maximum absolute correlation of -0.25). An ANOVA test demonstrated that in the case of binning there is a significant decrease in amplitude with an increase of the PB dose for all three bins. The ANOVA test for deconvoluted peaks demonstrated a completely different picture, with the glycine peak significantly increasing in intensity while neighboring peaks decrease (Figure 1b-e). The increase in glycine with the increase of PB exposure was confirmed by parallel gas chromatography/mass spectrometry (GC/MS) analysis of aqueous fraction of liver tissue as detailed in ref 30. The ability of the Bayesian model to improve resolution can be further demonstrated with examination of resonances associated with L-gluthathione. The multiplets at 2.16 and 2.53 ppm are highly visible and can serve as “internal controls” for testing the ability of an approach to identify concentration changes while the multiplet at 3.77 ppm is obscured by glucose resonances and therefore often “lost” from the analysis. The 2.16 and 2.53 ppm multiplets were confidently detected with both approaches, with a statistically significant increase associated with an increase in the PB dose. However, the profiles of the bins containing the multiplet at 3.77 ppm do not allow the detection of a change in the concentration of glutathione (Figure 2). The Bayesian model successfully deconvoluted the gluthathione resonances from glucose resonances and allowed for the identification of a statistically significant increase in the treated samples. False Discovery Rates. There is a multiple comparison problem associated with large numbers of univariate tests as might

be used to determine which variables are different in a metabolomic experiment. It is not possible to obtain specificity and sensitivity characteristics for the two sets of tests as the true positives are not known. However, we can estimate and compare expected sensitivity properties with the false discovery rate (FDR) approach (exchangeable with “q-values” here after33) that is used to control the number of false positives in the multiple testing. In order to compare the ability of nonparametric spectral binning and Bayesian quantification to provide quality candidates for biomarker identification we applied the robust Wilcoxon rank sum tests to the integral values produced by binning and the Bayesian estimates of peak amplitudes corresponding to the basal group (n ) 5) and 1000 ppm dose group (n ) 5) on days 1, 3, 7, and 14. The q-value versus p-value plots demonstrate that the quantification provides a consistently lower FDR than binning (Figure 3) with respect to a chosen p-value cutoff. It is worth noting that the difference in q-values is not a consequence of the difference in the total number of variables but rather the difference in the corresponding p-values’ distributions.33 The day 14 data set was chosen for further study of the properties of the two approaches. A total of 93 peaks and 179 buckets were marked as significantly different. However, the FDR of the quantified peaks was no more than 8% while the significant bins had a FDR at 16%. This suggests that the quantification provides more power for statistical testing as the equivalent p-value, while allowing for fewer false positives. While it is not possible to directly compare the differently expressed variables between the two approaches as one spectral bin can correspond to several peaks and vice versa, we still can compare the relative “quality” of each variable by comparing corresponding q-values. Figure 4 shows that the quantification has less risk of false positives in terms of the proportion of the variables accepted as significantly different between the control and dosed groups. For example, if we accept the most differentially expressed 10% of variables as significant then we will incur an FDR of 4.2% in the case of binning and 2.7% FDR in the case of quantification. Further Analytical Chemistry, Vol. 82, No. 11, June 1, 2010

4483

Figure 2. Identification of a change in the 3.77 ppm gluthathione resonance. Part a shows deconvoluted peaks in the region of 3.73-3.80 ppm. The solid line boxes illustrate average values of the free bins covering most of the intensity of the assigned gluthatione peak at 3.77 ppm. Parts b, c, and d show the boxplots summarizing the intensity of the bins containing a glutathione resonance at 3.77 ppm across different PB dose levels. The ANOVA test shows that there is no change (p-value of 0.1154) or a significant (p-values of 7 × 10-5 and 1 × 10-7, respectively) decrease between the control and at least one of the treated samples. Part e shows that a modeled resonance for glutathione in fact demonstrates the opposite direction of change, with a p-value 0. 0.0016, for an increase between the control and the 1000 ppm PB dose groups. The peak was assigned to glutathione based on an NMR spectrum of a standard and correlation with the 2.16 ppm multiplet (0.824 with p-value of 1 × 10-20) and correlation with the 2.53 ppm multiplet (0.77 with a p-value of 5 × 10-17).

Figure 3. False discovery rates corresponding to different p-value cut offs. The figures show the expected number of false positives (q-value) for a particular p-value cut off for all four time points.

exploration reveals that 54 of the most significant peak amplitudes have FDR less or equal to 0.027 while the minimal q value for significant bins is 0.041 (Figure 4) DISCUSSION In this article we present a novel application of Bayesian deconvolution for 1D NMR spectra to an actual metabolomic study of liver toxicity response to acute PB exposure. The motivation for this work has been twofold. One reason has been a necessity to improve the quality of NMR processing and therefore increase the utility of NMR as a biomarker discovery tool. At the same time a higher level of processing automation is needed to match increasingly bigger scale metabolomic studies with hundreds of samples routinely analyzed. In our view using the Bayesian modeling allows us to both improve the quality of analysis while 4484

Analytical Chemistry, Vol. 82, No. 11, June 1, 2010

minimizing user input. The identified changes in many metabolites agreed well with published results. In particular, the change in glycine, which was confirmed by GC/MS had not been detected by NMR spectroscopy using classical binning methods.30 A major quality improvement in the analysis pipeline over the widely used binning approach is demonstrated in terms of increased resolution and better statistical properties of the resulting set of variables that facilitates a more efficient biomarker discovery process. The Bayesian approach is model-based and therefore more powerful than nonparametric ones, especially in the case of closely spaced peaks. Prior knowledge about peak shape informs the deconvolution process and allows the extraction of information in an automated way and with higher resolution. The use of a trans-dimensional sampling technique permits

Figure 4. Comparison of the false discovery rates for Wilcoxon rank sum tests for the day 14 samples. The y-axis shows a number of variables assumed to be significantly different between the control and 1000 ppm dose groups. The x-axis shows the corresponding FDR. The bold lines and filled circle markers highlight a p-value cut off of 0.05 as a significance threshold and the empty circles mark a stricter threshold of 0.01. The expected FDR levels for Bayesian deconvolution and binning are 0.027 and 0.04, respectively.

automation at an even higher level as it does not require user decision on the important subject of model order. A previous study has shown that the proposed method is superior to such wellknown criteria as minimum description length (MDL) and Akaike information criterion (AIC) for model order determination.29 The cost of improved analysis is a higher computational burden of fitting parametric models. However, this is alleviated by an increasingly powerful computational capability becoming widely available. The Bayesian processing of the series of the localized models and subsequent generation of a global model for one 32K spectrum took about 30 min on a Pentium 5 PC with 2 GB of RAM. Actual benefits of improved resolution and the potential reduction in valuable user time needed for analysis may further justify a longer quantification process. The software program is written in Matlab script and available from the authors by request. In order to reduce the complexity of modeling, a frequencyselective approach was applied similar to the ones described in refs 37 and 38. The frequency window can potentially affect quantification of broad peaks when their power is spread across several windows. However, in metabolomic experiments, the sharp lines representing small molecules are usually of primary interest and decaying sinusoids are reasonably well localized in the (37) Mandelshtam, V. A. Prog. Nucl. Magn. Reson. Spectrosc. 2001, 38, 159– 196. (38) Sandgren, N.; Selen, Y.; Stoica, P.; Li, J. J. Magn. Reson. 2004, 168, 259– 272. (39) Stoyanova, R.; Nicholls, A. W.; Nicholson, J. K.; Lindon, J. C.; Brown, T. R. J. Magn. Reson. 2004, 170, 329–335. (40) Tomasi, G.; van den Berg, F.; Andersson, C. J. Chemom. 2004, 18, 231– 241.

frequency domain for small damping factors. The boundaries of the windows can be either chosen manually or automatically based on local minima. Different ways of choosing a set of windows did not affect the deconvolution results as long as the windows were overlapping and had reasonable width to cover the bulk of intensities for the peaks in the central, non down-weighted region. We used a separate restricted maximum likelihood procedure to obtain the final point estimation of the peak amplitudes simultaneously for all frequency windows combined. This allows one to take into account possible interference between neighboring sections of a spectrum. Consistency in chemical shift positions of resonances across a data set, although not affecting the Bayseian deconvolution itself, is essential for further identification and statistical analysis. We found that simple automated spectral alignment worked well for the relatively stable spectra produced by the liver tissue samples (details are given in the Supporting Information). However, variation in peak positions in less consistent samples such as urine can require the application of more sophisticated methods such as PCA-aided alignment39 or correlation warping.40 Given the reasonably consistent peak positions, further acceleration of the processing is achieved by applying the deconvolution algorithm to a sum of all spectra in a data set so a resulting set of peaks provides a common basis for comparison of different samples and facilitates assignment and interpretation. The proposed set of methods worked well for the 1D NMR spectra of liver tissue samples, but its applications are not limited to tissue samples and toxicological applications. If the problem of peak shifts is addressed as described above, the same approach could be applied to biofluid or bacterial culture samples as well. The combination of frequency selective filters and accurate peak modeling makes it an attractive method for targeted filtering and quantification of particular peaks of interest as, for example, the lactate methyl doublet which is coresonant with part of the lipoprotein fraction in blood plasma spectra. It therefore may find application in toxicological, epidemiological, and dietary studies as well as disease diagnostics and monitoring. ACKNOWLEDGMENT This work was funded by the BBSRC, U.K. (Grant BB/ H013539/1) and Syngenta. J.L.G. was also supported by a Royal Society University Research fellowship. SUPPORTING INFORMATION AVAILABLE Additional information as noted in text. This material is available free of charge via the Internet at http://pubs.acs.org.

Received for review February 7, 2010. Accepted April 5, 2010. AC100344M

Analytical Chemistry, Vol. 82, No. 11, June 1, 2010

4485