Anal. Chem. 2004, 76, 3666-3674
Sample Classification Based on Bayesian Spectral Decomposition of Metabonomic NMR Data Sets Radka Stoyanova,*,† Jeremy K. Nicholson,‡ John C. Lindon,‡ and Truman R. Brown§
Fox Chase Cancer Center, 333 Cottman Avenue, Philadelphia, Pennsylvania 19111, Biomedical Sciences Division, Biological Chemistry, Faculty of Medicine, Imperial College London, Sir Alexander Fleming Building, South Kensington, London, SW7 2AZ, U.K., and Hatch Center for MR Research, Columbia University, 710 West 168th Street, New York, New York 10032
1H
NMR spectra of biofluids provides a wealth of biochemical information on the metabolic status of an organism. Through the application of pattern recognition and classification algorithms, the data have been shown to provide information on disease diagnosis and the beneficial and adverse effects of potential therapeutics. Here, a novel approach is described for identifying subsets of spectral patterns in databases of NMR spectra, and it is shown that the intensities of these spectral patterns can be related to the onset and recovery from a toxic lesion in both a time-related and dose-related fashion. These patterns form a new type of combination biomarker for the biological effect under study. The approach is illustrated with a study of liver toxicity in rats using NMR spectra of urine following administration of a model hepatotoxin hydrazine. 1H NMR spectroscopy of biofluids when coupled with the application of computer-based pattern recognition techniques (PR) have been successfully used to classify biofluids in terms of known pathological effects caused by the administration of toxic substances.1-5 1H NMR of urine allows the simultaneous assessment of thousands of resonances resulting from the presence of both endogenous biochemicals and drug metabolites related to the effects of a wide range of xenobiotic compounds.6 Studies have shown that urine samples obtained from rats following the administration of hepatotoxins, nephrotoxins, and testicular toxins
* To whom correspondence should be addressed. Fax: 215-728-1622. Email:
[email protected]. † Fox Chase Cancer Center. ‡ Imperial College London. § Columbia University. (1) Holmes, E.; Nicholls, A. W.; Lindon, J. C.; Connor, S. C.; Connelly, J. C.; Haselden, J. N.; Damment, S. J.; Spraul, M.; Neidig, P.; Nicholson J. K. Chem. Res. Toxicol. 2000, 13, 471-478. (2) Holmes, E.; Nicholls, A. W.; Lindon, J. C.; Ramos, S.; Spraul, M.; Neidig, P.; Connor, S. C.; Connelly, J.; Damment, S. J. P.; Haselden, J.; Nicholson, J. K. NMR Biomed. 1998, 11, 235-244. (3) 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. (4) Anthony, M. L.; Sweatman, B. C.; Beddell, C. R.; Lindon, J. C.; Nicholson, J. K. Mol. Pharmacol. 1994, 46, 199-211. (5) Nicholson, J. K.; Higham, D. P.; Timbrell, J. A.; Sadler, P. J. Mol. Pharmacol. 1989, 36, 398-404. (6) Lindon, J. C.; Nicholson, J. K.; Everett, J. R. Annu. Rep. NMR Spectrosc. 1999, 38, 1-88.
3666 Analytical Chemistry, Vol. 76, No. 13, July 1, 2004
can be distinguished on the basis of 1H NMR-detected perturbations in the urine.2,7 The need for early identification of drug toxicity in drug discovery and development, and the concept that toxins that induce lesions in the same tissue or via similar mechanisms will generally produce a similar profile set of biochemical changes in the urine 8 have motivated a large scale investigation in identifying the spectral changes caused by a range of toxins.9 High-resolution 1H NMR urine spectra present a significant challenge for the traditional peak-quantifying techniques for spectral analysis, because they contain thousands of partially overlapping resonances arising from hundreds of endogenous molecules. The success of detection and quantification of the changes in NMR spectroscopic metabolite patterns is due to the implementation of powerful pattern-recognition techniques that view each data set as a whole, and allowing the discrimination between different combinations of metabolites characteristic of a biological system and its dynamic changes. Since the first application of PR to the analysis of biological NMR data,10 numerous approaches have been developed for analysis of large spectral data sets: principal component analysis (PCA)-based metabolic trajectories 3, partial least squares (PLS)-based batch statistical process control,11 hierarchical cluster analysis and k-nearest neighbor approaches,12 and a density superposition approach, “CLOUDS”.13 The major goal of these techniques is the classification of biofluid samples on the basis of differences between the unaltered and the “aberrant” spectral patterns related to the administration of a particular toxin, thus providing new means for predicting the toxic effects of novel foreign compounds. These approaches have also been employed in other areas of (7) Gartland, K. P. R.; Beddell, C. R.; Lindon, J. C.; Nicholson, J. K. J. Pharm. Biomed. Anal. 1990, 8, 963-968. (8) Gartland, K. P. R.; Bonner, F. W.; Nicholson, J. K. Mol. Pharmacol. 1989, 35, 242-250. (9) 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. Toxicol. Appl. Pharmacol. 2003, 187, 137-146. (10) Gartland, K. P.; Sanins, S. M.; Nicholson, J. K.; Sweatman, B. C.; Beddell, C. R.; Lindon, J. C. NMR Biomed. 1990, 3, 166-172. (11) Antti, H.; Bollard, M. E.; Ebbels, T. M. D.; Keun, H. C.; Lindon, J. C.; Nicholson, J. K.; Holmes, E. J. Chemom. 2002, 16, 1-9. (12) Beckonert, O.; Bollard, M. E.; Ebbels, T. M. D.; Keun, H. C.; Antti, H.; Holmes, E.; Lindon, J. C.; Nicholson, J. K. Anal. Chim. Acta 2002, CAC2002 procedings. (13) Ebbels, T.; Keun, H.; Beckonert, O.; Antti, H.; Bollard, M.; Holmes, E.; Lindon, J.; Nicholson, J. Anal. Chim. Acta, in press. 10.1021/ac049849e CCC: $27.50
© 2004 American Chemical Society Published on Web 05/21/2004
application, such as clinical diagnosis,14 and have used other types of analytical technology, such as mass spectrometry.15 Here, we present a different approach for class separation and assignment: Bayesian spectral decomposition (BSD). This allows the identification of a series of subspectral features inherent in the total spectra and determination of their strengths over the course of an experiment.16 A major advantage of BSD is that the identified spectral patterns have a direct metabolic interpretation, since this decomposition of the spectral data set produces spectral shapes with positive values. If a pattern consists of a single metabolite, its strength in the decomposition will be proportional to its concentration in the samples. The patterns in metabonomic data are more complex, containing several metabolites, and thus, we will refer to the strength of a particular pattern as amplitude. The changes in the amplitude of the spectral shapes can be related to other information known about the samples, such as timecourse, dose, etc., and subsequently can be subjected to any clustering algorithm for class separation. A common issue when applying PR is how to decide on the number of the sought patterns, features, and classes. In this paper, we also explore the effect of the different number of BSD solutions on data interpretation. Spectral preprocessing, prior to application of PR, often predetermines the success of their performance. Typically, the biologically significant spectral changes are related to altered concentrations of the different analytes, and hence, any fluctuations in peak position, for example, can have a deleterious effect on PR performance. To avoid this, the PR approaches are often applied to “binned” data where spectral regions are aggregated. However, here BSD is applied to the original spectral frequency data, following an alignment procedure as described previously.17 The application of BSD is illustrated on 1H NMR urine spectra from rats treated with different doses of hydrazine over a 150-h period. Hydrazine has been commonly used as a model hepatotoxin in animal studies.18 Spectral shapes that are related to the unaltered and “aberrant” urine patterns were identified with BSD. The strength of these patterns is readily correlated with the dose and time of hydrazine administration, and these can also be successfully subjected to classification. EXPERIMENTAL SECTION Data Description. High-resolution 1H NMR spectroscopy of urine was used to investigate the systematic biochemical effects of oral hydrazine administration in male Han Wistar rats.18 Twenty male rats (200-250 g) were housed individually in metabolism cages. The rats were assigned to four dose groups: control, 75, 90, and 120 mg/kg hydrazine, and each animal was administered either hydrazine in sterile water or the vehicle alone via oral gavage. Urine samples were collected at predose (-8 to 0 h), and (14) Brindle, J. T.; Antti, H.; Holmes, E.; Tranter, G.; Nicholson, J. K.; Bethell, H. W.; Clarke, S.; Schofield, P. M.; McKilligin, E.; Mosedale, D. E.; Grainger, D. J. Nat. Med. 2002, 8, 1439-1444. (15) Plumb, R.; Granger, J.; Stumpf, C.; Wilson, I. D.; Evans, J. A.; Lenz, E. M. Analyst 2003, 128, 819-823. (16) Ochs, M. F.; Stoyanova, R. S.; Arias-Mendoza, F.; Brown, T. R. J. Magn. Reson. 1999, 137, 161-176. (17) Stoyanova, R.; Brown, T. R.; Lindon, J.; Nicholls, A.; Nicholson, J. Proc. Int. Soc. Magn. Reson. Med. 2001. (18) Nicholls, A. W.; Holmes, E.; Lindon, J. C.; Shockcor, J. P.; Farrant, R. D.; Haselden, J. N.; Damment, S. J. P.; Waterfield, C. J.; Nicholson, J. K. Chem. Res. Toxicol. 2001, 14, 975-987.
Table 1. Data Set Consisting of 57 Urine 1H NMR Spectra from Treated and Control Rats, Collected at the Time Points Marked with “+” hours group/dose
0
8
16
24
32
56
80
104
128
152
control control control 75 mg/kg 90 mg/kg 120 mg/kg 120 mg/kg
+ + + + + + +
+ + + + + +
+ + + + + +
+ + + + + +
+ + + + + +
+ + + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
at 0-8, 8-16, 16-24, 24-32, 48-56, 72-80, 96-104, 120-128, and 144-152 h postdose. Here, each treatment period is indicated by its upper limit; i.e., 0-8 h is 8 h. The samples were prepared as described previously,18 and 50 µL of a solution of sodium 3-trimethylsilyl-(2,2,3,3-2H4)-1-propionate (TSP) in D2O was added to the urine samples. The addition of TSP and D2O provides both a chemical shift reference (0.0 ppm) and a deuterium lock signal for the NMR spectrometer. One-dimensional 1H NMR spectra of urine were measured at 600.13 MHz on a Bruker DRX-600 spectrometer. Sixty-four FIDs were collected into 64 K points using a spectral width of 7003 Hz, an acquisition time of 4.68 s, and a total pulse recycle time of 7.68 s. The composition of the analyzed data set is presented in Table 1. Spectra were line-broadened by 0.3 Hz, transformed to the frequency domain, phased, and baseline-corrected. A typical 1H NMR spectrum from one of the control rats is shown in Figure 1. For normalization of the urine signal, the integral of the creatinine [4.00-4.1 ppm] (Figure 1) was used. The area of the peak was estimated after correction of local background by fitting a straight line into the creatinine region and subtracting it from the signal. Individual Peak Alignment. The spectral peaks in the data set were aligned following a procedure described earlier.17 When principal component analysis (PCA) is applied to a spectral data set, the peak derivative shapes appearing in the second- and higher-order principal components (PCs) are indicative of the presence of shifts in the corresponding peak regions.19 To identify the peak regions where these shifts occur, derivatives of Lorentzian triplets, doublets, and single lines were simulated with varying line widths and were correlated with the second- and higher-order PCs at each frequency so that any underlying derivative shapes are detected. The peak region with the highest correlation was isolated, and the frequencies of the highest points were aligned. The procedure was applied to all peaks with a significant (>0.80) correlation coefficient. Bayesian Spectral Decomposition. Bayesian spectral decomposition is described in detail elsewhere.16 Briefly, it is assumed that the spectra in a data set are not independent of one another, but that they are rather a mixture of a small number of basic spectra coming from different subcomponents of the studied sample with variable contributions to the individual spectra. It is assumed that the number of these basic shapes is k and F (k × (19) Stoyanova, R.; Brown, T. R. NMR Biomed. 2001, 14, 271-277.
Analytical Chemistry, Vol. 76, No. 13, July 1, 2004
3667
Figure 1. A typical 600-MHz 1H NMR spectrum of urine from a control rat. The water peak has been suppressed, and the resonance at 0 ppm is the reference signal from TSP. The total integral of the CH2 peak of creatinine [4.0-4.1 ppm] was used to normalize the spectral intensities.
m) is defined as a matrix containing these spectra, where m is the number of frequency points in the data. Then, BSD determines the following representation of the data matrix D (n × m), containing n spectra
D ≈ AF
(1)
where A (n × k) (k , m, n) contains the mixing coefficients of the basic spectra in each spectrum. The conditions for the decomposition in eq 1 are as follows: F are physical spectral shapes, and as such, they consist of positive peaks. Thus, the spectra in F are not necessarily orthogonal (as in the case of PCA). They are orthogonal only when each peak is present in only one of the shapes. The mixing coefficients, A, are also nonnegative, since a shape cannot have a negative contribution to the measured spectrum. Since neither A nor F is known, the number of possible solutions is large and cannot be identified analytically. Often in such cases, Markov chain Monte Carlo (MCMC) methods are used to sample the space of possible solutions to determine its properties. MCMC methods require probability measurements in each sampled point in the solution space, which in BSD is provided through a Bayesian approach. In the terms of the decomposition in eq 1, Bayes’ theorem may be rewritten as
p(A, F|D) ∝ p(D|A, F)p(A, F)
(2)
where p(A, F|D) is the conditional probability for the model to be true, given the data (the posterior); p(D|A, F) is the conditional probability that the data are true, given the model (the likelihood); and p(A, F) is the probability of the model independent from the data (the prior). Equation 2 provides the MCMC with the needed probability at a point in the solution space for the problem. The likelihood can be easily determined by comparing the data reconstructed by the model (the right-hand side in eq 2) to the data in D. The prior encodes any knowledge of the model. For example, a reasonable prior for a solution is that a higher probability is given to a series of narrow features, as compared to a flat offset. BSD utilizes the massive inference Gibbs sampler (MaxEnt Solutions Ltd., Cambridge, U.K.) to sample from the 3668
Analytical Chemistry, Vol. 76, No. 13, July 1, 2004
posterior distribution and for the encoding of the prior.20 Massive inference uses an “atomic” prior, which represents the model as a few point fluxes (atoms), with the highest probability assigned to the distributions with the fewest atoms. The algorithm also enforces positivity on the solutions.21 Determining the number, k, of the sought components is not straightforward. In Ochs et al. as well as in other PR approaches, the choice of k is based on the number of significant PCs in the data. Once A and F are determined, A can be related to any additional information known about the data: time-course, dose, sample composition, etc. The spectral patterns in F have been referenced on the usual NMR chemical shift parts-per-million scale for ease of identification of the peaks. Class Segregation. For the purposes of data classification, any standard classification procedure [hierarchical cluster analysis; k-nearest neighbor (kNN) classification, etc.] applied to the amplitudes of the patterns in matrix A will result in a reliable classification. Hierarchical cluster analysis (HCA), for example, is an unsupervised pattern recognition method that permits the grouping of data points that are similar by virtue of being “near” to one another in a multidimensional space and whose coordinates are defined by the correlations of the amplitudes in matrix A.22 The closest pair of points are reported with their separation distance, and then the two points are deleted and replaced with a single combined point. The process is then repeated iteratively until only one point remains. The reported connectivities are then plotted as a dendrogram (a tree-like chart, which allows visualization of clustering), showing sample-to-sample connectivities versus an increasing separation distance. A property of the dendrogram is that the branch lengths are proportional to the distances between the various clusters, and hence, the length of the branches linking one sample to the next is a measure of their similarity. This way, similar data points may be identified algorithmically. Alternatively, any other suitable classification algorithm could be used. (20) Geman, S.; Geman, D. IEEE Trans. Pattern Anal. Mach. Intell. 1984, 6, 721-741. (21) Sibisi, S.; Skilling, J. J. R. Stat. Soc. 1997, B 59, 217-235. (22) Everitt, B. Cluster Analysis; Halsted Press: New York, 1980.
Figure 2. The first four PCs from the spectral region between 2.20 and 3.63 ppm (a) before and (b) after application of the procedure for individual peak position alignment, together with their normalized eigenvalues. The percentages indicate the variance explained by each PC. The peaks, aligned by the procedure are marked with asterisks in the third and fourth PC. (c) A scree plot of the first 10 normalized eigenvalues of the spectral region before (open diamonds) and after (filled circles) the alignment procedure. (d) Plot of the spectra from the corrected-forfrequency-shifts hydrazine data set in the plane, defined by the first two principal components. The data from control and the three dose groups are indicated by different symbols.
RESULTS The region of the urine 1H NMR spectra that has been analyzed lies between 2.20 and 3.63 ppm (8051 points), because this is the richest with respect to metabolite peaks. PCA was applied to the real part of the data-matrix D (57 × 8051), and the PCs were calculated as eigenvectors of the data-covariance matrix, calculated around the origin.23 In Figure 2a, the first four PCs of this spectral region and their corresponding normalized eigenvalues are presented. It is clear that although the first two PCs are entirely related to changes in the peak amplitudes, the third and fourth are composed of derivative shapes, indicating the presence of peak shifts in the data. The procedure for individual peak-alignment was applied by detecting derivative shapes in the third PC. A total of six peak regions (marked with asterisks) were isolated and corrected. Subsequent PCA indicated that there were still derivative shapes (in the fourth PC; data not shown), and two more peaks were aligned (also marked) in Figure 2a. The results of (23) Stoyanova, R.; Kuesel, A. C.; Brown, T. R. J. Magn. Reson., Ser. A 1995, 115, 265-269.
PCA after these corrections are shown in Figure 2b. The total variance in the data set decreased by 12%; the first PC now explains 64%, versus 52% in the initial data set. A marked improvement was achieved in the shapes of the peaks in the first PC, particularly in the region of the citrate resonance. There are, however, very few changes in the second PC. The derivative shapes disappeared entirely in the higher order PCs. The plot of the fraction of the total variance explained by each eigenvalue versus the PC number, shown in Figure 2c, also indicates the effect of the alignment procedure. As mentioned above, the number of significant PC can be used to determine k, the number of sought patterns in BSD. In many cases, this number is uncertain, and the scree plot in Figure 2c indicates that this is particularly true for the uncorrected-for-frequency-shifts initial data set, for which the fraction of total variance, explained by each eigenvalue (open diamonds) declines in monotonic fashion. On the other hand, the sharp decline in the slope of the eigenvalues of the corrected data set (filled circles) indicates that at most, only a few significant PCs are needed to describe the data set. Analytical Chemistry, Vol. 76, No. 13, July 1, 2004
3669
Figure 3. (a) Amplitude of two spectral patterns, identified by BSD for k ) 2, across the treated and control rat urine NMR spectra as a function of time after administration of a single dose of hydrazine: filled diamonds are associated with the “normal” urine pattern, and open circles are associated with the “aberrant” pattern (last two points of treatment with 120 mg/kg are off the scale of the graph). (b) Amplitude of the “normal” pattern (left panel) and the “aberrant” pattern (right panel) in the treated rats as function of time. (c) Spectral patterns, identified by BSD in the spectral data set after the alignment correction, “normal” (upper) and “aberrant” (lower).
Below, results for k ) 2, 3, and 4 are described. In Figure 2d, the plot of the spectra from the corrected-for-frequency-shifts hydra3670 Analytical Chemistry, Vol. 76, No. 13, July 1, 2004
zine data set is presented in-plane, defined by the first two principal components. The data from the control and three dose groups
Figure 4. Analysis by BSD as for Figure 3, but instead for k ) 3. Patterns in addition to those seen in Figure 3 are identified by BSD for k ) 3. (a) Amplitude of the pattern in the treated rats as a function of time. (b) Corresponding spectral pattern.
are indicated by different symbols, and clear separation between the treated and control rats can be observed. BSD was applied to the real portion of the spectra in the corrected data set, initially looking for two solutions (k ) 2). The amplitudes of the two identified patterns, A (57 × 2), as a function of time for the control rats and for each dose group are presented in Figure 3a [the 0-h point of the second rat treated with 120 mg/ kg (last row of Table 1) is not shown]. First, the separation between the treated and control rats is well-demonstrated. Although the first shape (filled diamonds) maintains almost a constant high level in the nontreated rats, the second (open circles) is close to 0 throughout the experiment for this group. Thus, the first pattern is associated with the urine unaffected by the hydrazine, and the second is reflecting the changes due to the hydrazine administration. Here and in the rest of the paper, we will refer to the first pattern as “normal” and the second as “aberrant”. For the dosed groups, it can be seen that the “aberrant” pattern shows a peak with respect to time, whereas the “normal” pattern shows a corresponding drop at that time point. The amplitudes of the two patterns for the treated animals are shown in Figure 3b. The hydrazine effect “peaks” at 32 h for both the 75 and 90 mg/kg groups. Subsequently, the two groups recover to the control level at 100 h. The data from the 120 mg/kg dose indicated no sign of recovery at 56 h, at which point the animal was sacrificed.18 The spectral patterns F (2 × 8051) (Figure 3c), can be directly related to the biochemical changes due to the hydrazine effects: the “normal” pattern (upper row) contains
principally the peaks from citrate, succinate, 2-oxoglutarate, and trimethylamine-N-oxide (TMAO). Alternatively, the amplitude of the “aberrant” pattern consists mainly of 2-aminoadipic acid, β-alanine, creatine, and taurine. The BSD results for k ) 3 revealed very similar “normal” and “aberrant” patterns as in the previous solution plus an additional solution, the amplitude and shape of which are displayed in Figure 4. The behavior of the amplitude for the treated animals in Figure 4a suggests that this pattern increases with the uptake of hydrazine. The corresponding spectral patterns (Figure 4b) represent a mixture of some of the metabolites in the “normal” (citrate, succinate, and TMAO) and “aberrant” (creatine, taurine) spectral patterns from the k ) 2 solution. Finally, BSD with four solutions revealed the same “aberrant” patterns as in the case with three solutions; however, the “normal” pattern appears to be split. The associated spectral shapes and amplitudes for the treated groups are shown in Figure 5. In contrast with the “aberrant” pattern, they have very similar behavior among the three treated groups. Their spectral characteristics (Figure 5b) are mainly citrate and succinate in the upper one and 2-oxoglutarate in the lower. TMAO appears to be split among the two patterns. HCA was applied to the amplitudes in the matrix A for k ) 2, 3, and 4, and the corresponding dendrograms are presented in Figure 6. The control rats are denoted with NT (not treated) and colored in green, and the treated (D1, D2 ,and D3) correspond to the three administered hydrazine doses and colored in shades Analytical Chemistry, Vol. 76, No. 13, July 1, 2004
3671
Figure 5. Analysis by BSD as for Figure 3, but instead for k ) 4. (a) Amplitude of the two “normal” patterns obtained by BSD for k ) 4 in the treated rats as a function of time. (b) Corresponding spectral patterns.
of red. It is clear that for the cases with 3 and 4 solutions, the two major branches clearly separate the treated from the nontreated groups. HCA correctly assigns the initial measurement (0 h) taken prior to hydrazine administration, as well as the later data points of treatment (after 104 h) for the dosed rats, to the “nontreated” cluster. For k ) 2, there were five data points related to early uptake (8 h) and clearing (80 h) of the drug, which were grouped in the “no treatment” cluster. These samples actually correspond to time points when the manifested toxic effect of hydrazine is minimal. With the increase of the information in matrix A, the subgroupings in the “treatment” cluster indicate different stages of response to treatment. For example, the last two points (D3-32 and D3-56) from the data of treatment with 120 mg/kg, associated with the highest level of toxicity, fall within a separate branch. DISCUSSION In this paper, we present a novel application of Bayesian spectral decomposition for class separation and class discovery. In contrast to other PR techniques, this application readily 3672 Analytical Chemistry, Vol. 76, No. 13, July 1, 2004
identifies the spectral patterns associated with the different classes in the treatment and control data sets. Thus, the technique can be used not only for classification and for predicting the toxic effects of novel foreign compounds, but also as a tool for investigating the biochemical basis of the observed differences in the response to diverse toxins. BSD allows the decomposition of complex spectral data sets into underlying spectral shapes with positive spectral amplitudes corresponding to changes of combinations of individual compounds associated with the response to toxin administration. The BDS-identified spectral shapes can, thus, be directly correlated to the biochemical/metabolic mechanisms underlying the response, and hence, it may contribute to their further delineation. When two patterns are sought (k ) 2), BSD identifies a “normal” and an “aberrant” pattern, which exhibit a distinctly different temporal behavior in the treated animals, as compared to the nontreated controls. This difference is not limited only to the direction and amplitude of the changes, but also in their
Figure 6. Hierarchical cluster analysis applied to the amplitudes of BSD patterns for k equal to (a) 2, (b) 3, and (c) 4. The cells colored in green represent the control samples (NT, not treated) and the ones in shades of red, those treated (denoted with D1, D2, and D3 for doses of 75, 90, and 120 mg/kg and the corresponding times).
temporal occurrence and persistence: the “aberrant” pattern is characterized by relatively short-lived changes and rapid recovery, which may be explained by the character of the direct tissue damage. Notably, the difference in the amplitude of the “aberrant” pattern for the three treatment groups reflects the administered dose of hydrazine and, hence, its potentially direct effect in terms of toxicity and extent of tissue damage. In contrast, for the changes in the “normal” pattern, although also quite rapid immediately
following the treatment, their “recovery” phase is significantly delayed. The “normal” spectral pattern contains intermediates of the Krebs cycle, including the peaks of citrate, succinate, and 2-oxoglutarate, which indicates the deregulation of this pathway in response to the toxin and its slower recovery, which is associated with complex anabolic and catabolic processes. For k ) 3, a third pattern is identified with a characteristic mixed behavior, containing both “aberrant”- and “normal”-associated spectral peaks. In contrast to the first two, “aberrant” and “normal” patterns, the “mixed-aberrant” pattern appears to be doseindependent because it maintains a similar behavior for all treatment dose levels. This suggests that the “mixed-aberrant” pattern may be instead associated with biological processes that are indirectly and only moderately affected by the administered toxin, such as processes involved in the regeneration of cellular and tissue physiological integrity. Finally, for k ) 4, the “normal” spectral pattern is split into two: one consists entirely of 2-oxoglutarate, and the second consists of the remaining Krebs cycle intermediates. The observed differences in the temporal behavior of these two “normal” patterns, and specifically the observed immediate depletion of 2-oxoglutarate, occurring prior to the decrease of succinate and citrate, suggest that the administered toxin triggers, as expected, distinct toxicity- and regenerationassociated mechanisms in a temporally sequential manner. Although the depletion of Krebs cycle intermediates in the one spectral pattern is rapid but relatively transient, the depletion of 2-oxoglutarate is quite persistent. It was previously suggested that the early loss of 2-oxoglutarate in response to hydrazine, which precedes that of citrate and succinate, is induced by Krebs cycleindependent mechanisms.18 As shown earlier, hydrazine conjugates directly with 2-oxoglutarate to form the metabolite THOPC.24 Additionally, hydrazine is known to lead to sequestration of pyridoxal 5′-phosphate (PLP), a cofactor of aminotransferase and decarboxylase enzymes.25,26 PLP sequestration may lead to an upregulated conversion of pyridoxamine 5′-phosphate to PLP, which is dependent on the conversion of 2-oxoglutarate to glutamate,27 thus also leading to depletion or loss of 2-oxoglutarate. Most of the metabolic changes occurring following the treatment with hydrazine have been previously described qualitatively by employing a number of biochemical or PR methods. BSD, however, not only quantifies directly the temporal changes taking place, but it also readily provides the identification of metabolites associated with them. This is, in part, due to a thorough preprocessing of the data for alignment of the individual peaks, which thereby allows the analysis of all frequency points of the spectral region, rather than merely aggregating the intensity data into bins. The procedure for individual peak alignment, as demonstrated in the paper, is quite successful in removing the frequency variations in the data. Although other alignment techniques have been proposed,28,29 their main drawback is the (24) Sanins, S. M.; Timbrell, J. A.; Elcombe, C.; Nicholson, J. K. Arch. Toxicol. 1992, 66, 489-495. (25) Newsholme, E. A.; Leach, A. R. Biochemistry for the Medicinal Sciences; John Wiley and Sons: Chichester, 1995. (26) Seigel, G. J. Basic Neurochemistry: Molecular, Cellular and Medical Aspects; Raven Press: New York, 1993. (27) Alberati-Giani, D.; Malherbe, P.; Kohler, C.; Lang, G.; Kiefer, V.; Lahm, H. W.; Cesura, A. M. J. Neurochem. 1995, 64, 1448-1455. (28) Forshed, J.; Schuppe-Koisyinen, I.; Jacobson, S. P. Anal. Chim. Acta 2003, 487, 189-199.
Analytical Chemistry, Vol. 76, No. 13, July 1, 2004
3673
assumption of the presence of the same peaks in the spectra to be aligned. In many metabonomic studies, the spectra from perturbed systems can show complete loss of some peaks or appearance of new peaks. It should be noted that BSD application to the uncorrected spectral data (results not shown) did not yield such crisp separation between the “normal” and “aberrant” spectral patterns; the identified spectral shapes were related mainly to frequency shifts present in the data. PCA is an inseparable part of the investigation of large, complex data sets, both in terms of aiding the determination of the number of significant components in the data and indicating the sources of variation. “Turning” the abstract, mathematical principal components into factors with physically meaningful content is not trivial, especially in the absence of some intuitive criteria, such as simple structure, parsimony, etc.30 An important part of the methods to resolve spectral mixtures is the presence of a pure variable, defined as a variable, that has intensity contribution from only one of the components of the data set.31 Techniques such as the method for self-modeling mixture analysis (SIMPLISMA) can resolve overlapping structures, but they also require some user interaction.32 BSD “resolves” the spectral patterns only on the basis of positivity of the solutions. In the case of the hydrazine data set, it is clear that, in the two-source solution, the peaks of citrate and aminoadipic acid can serve to identify the signatures of the “normal” and “aberrant” patterns. In the cases with a higher number of sources, this is no longer true. Further, in the rather complex in vivo and in vitro NMR spectra, such prior information is not generally available, underscoring the importance of the exploratory findings of BSD. Once the spectral patterns are identified, a reasonable hypothesis about the association of particular biochemical state with the presence or absence of a given metabolite can be formulated and tested. In the absence of any other constraints but positivity of A and F, the uniqueness of the solutions for the decomposition in eq 1 cannot be proved analytically. In fact, in most cases of selfmodeling curve resolution methods, the following word of caution is obligatory: “Without a theoretical model to base the results on, an infinite number of reasonable solutions or ‘best fits’ exist; (29) Eilers, P. H. C. Anal. Chem. 2004, 76, 404-411. (30) Malinowski, E. R. Factor Analysis in Chemistry; John Wiley & Sons: New York, 1991. (31) Hamilton, J. C.; Gemperline, P. J. J. Chemom. 1990, 4, 1-13. (32) Windig, W.; Guilment, J. Anal. Chem. 1991, 63, 1425-1432.
3674
Analytical Chemistry, Vol. 76, No. 13, July 1, 2004
while constraints can narrow the set of feasible solutions somewhat, physically correct solutions are never guaranteed”.31 The question for the BSD is whether there are other solutions for a given k that maximize the posterior probability in eq 2. To address this issue for any k described in the paper, multiple runs of BSD initiated with different random seeds were executed. In each case, the solutions were practically identical. Although not a formal proof of the uniqueness of the solution, this approach indicates their stability. Exploration of the effect of the number of sought patterns on the retrieved information demonstrated that a larger number of solutions provides, in addition to a more accurate classification, a more detailed and interpretable description of the involved metabolic processes. For example, the differential time response to hydrazine shown for the two “normal” patterns, found for k ) 4, distinguishes the observed rapid decrease of 2-oxoglutarate from that of other Krebs cycle intermediates that alter over a different time scale. This almost certainly indicates the involvement of at least two independent mechanisms. It is feasible to expect that to fully explore the behavior in these series of spectra requires even more solutions. Yet, the solutions for k ) 4 appear to capture the main variation in response to the hydrazine administration. The presented results demonstrate the potential of BSD, when applied in combination with correct and rigorous preprocessing of the data, to identify underlying basic spectral shapes within a complex set of spectra and to connect them quantitatively with the biological end-point measurements, as well as to accurately classify the samples. This approach shows promise for deconvoluting complex metabolic responses from disease, drugs, and toxins and should find application to identification of biomarker combinations for disease diagnosis and for monitoring the beneficial and adverse effects of pharmaceutical compounds. ACKNOWLEDGMENT The authors thank Dr. A. Nicholls for his help with data preparation and Drs. C. Patriotis and F. Kappler for critical review of the manuscript. The work described in this report was supported in part by funds provided through NIH grant PO1CA41078 to T. R. Brown. Received for review January 26, 2004. Accepted March 26, 2004. AC049849E
