Comparing Capillary Electrophoresis−Mass Spectrometry Fingerprints

Oct 27, 2008 - A set of interactive exploratory tools, utilizing principal component analysis (PCA) and analysis of variance (ANOVA) results based on ...
0 downloads 0 Views 1MB Size
Download PDF
Anal. Chem. 2008, 80, 8946–8955

Comparing Capillary Electrophoresis-Mass Spectrometry Fingerprints of Urine Samples Obtained after Intake of Coffee, Tea, or Water Erik Allard, Daniel Ba ¨ ckstro¨m, Rolf Danielsson, Per J. R. Sjo ¨ berg, and Jonas Bergquist* Department of Physical and Analytical Chemistry, Analytical Chemistry, Biomedical Centre, Uppsala University, P.O. Box 599, SE-75124 Uppsala, Sweden Metabolomic fingerprinting is a growing strategy for characterizing complex biological samples without detailed prior knowledge about the metabolic system. A twoway analysis system with liquid separation and mass spectrometric detection provides detail-rich data suitable for such fingerprints. As a model study, human urine samples, obtained after intake of coffee, tea, or water, were analyzed with capillary electrophoresis electrospray ionization time-of-flight mass spectrometry (CE-ESI-TOFMS). In-house-developed software (in Matlab) was utilized to manage and explore the large amount of data acquired (230 CE-MS runs, each with 50-100 million nonzero data points). After baseline and noise reduction, followed by suitable binning in time and m/z, the data sets comprised 9 and 14 million data points in negative and positive ESI mode, respectively. Finally, a signal threshold was applied, further reducing the number to about 100 000 data points per data set. A set of interactive exploratory tools, utilizing principal component analysis (PCA) and analysis of variance (ANOVA) results based on a general linear model, facilitated visual interpretation with score plots (for group assessment) and differential fingerprints (for “hot spot” detection). In the model study highly significant differences due to beverage intake were obtained among the 10 first principal components (p < 10-6 for two of the components in both ESI modes). Especially, the contrasts between “coffee” and “tea or water” indicated several “hot spots” with highly elevated intensities (e.g., for uncharged masses 93, 94, 109, 119, 123, 132, 148, 169, 178, 187, 190, and 193) suitable for further analysis, for example, with tandem MS. Metabolomics has been defined as “the quantitative measurement of the dynamic multiparametric metabolic response of living systems to pathophysiological stimuli or genetic modification”.1 Studies of the metabolome can be divided into three major different approaches as done by Shualev:2 (i) targeted analysis, assessment of a limited number of target analytes; * Corresponding author. Phone: +46 18 4713675. Fax: +46 18 4713692. E-mail: [email protected]. (1) Nicholson, J. K.; Lindon, J. C.; Holmes, E. Xenobiotica 1999, 29, 1181– 1189. (2) Shulaev, V. Brief. Bioinf. 2006, 7, 128–139.

8946

Analytical Chemistry, Vol. 80, No. 23, December 1, 2008

(ii) metabolite profiling, investigating a set of metabolites associated to a specific metabolic system or pathway; (iii) metabolic fingerprinting, using tools to acquire and explore metabolic fingerprints rather than specific metabolites. So far, most metabolic studies are focused on targeted analysis and, to some extent, metabolite profiling. However, the number of published methods suitable for collecting and interpreting metabolic fingerprint data is increasing, and we have reasons to believe that the number of applications in this category will increase over the next few years. Metabolic fingerprinting, which is the focus of this work, can be applied to any living system, such as cell cultures, bacteria, plants, or animals, without detailed a priori knowledge about the metabolic system. By comparing fingerprints from different groups or time points obtained from, for example, humans or animals, it has the potential to serve as a diagnostic tool for diseases, assist in drug screening and toxicology studies, or to address underlying pathophysiological and pharmacodynamical processes.3,4 As a further step it is often required that the chemical content of the “hot spots” in the fingerprint are adequately identified. Metabolic changes in the body are potentially reflected more or less clearly in various body fluids. In this study we made use of urine samples, where the analytes are present in a complex matrix consisting mainly of salts and polar compounds. One benefit of performing a study based on urine samples is that they can be noninvasively collected from humans. Urine contains a wide range of endogenous as well as exogenous metabolites and can therefore be assumed to provide an extensive insight in the metabolic state of an individual. High-resolution data in combination with sophisticated pattern recognition and data exploration tools are required to fully utilize the potential of the metabolic fingerprint strategy. Previous work in metabolomics has used a wide variety of analytical tools to acquire metabolic information, and a survey of these has been given in the recent reviews.5,6 Nuclear magnetic resonance (NMR), direct infusion mass spectrometry, gas chromatography/ mass spectrometry (GC/MS), and liquid chromatography-mass spectrometry (LC-MS) are all fairly common techniques used (3) (4) (5) (6) (7)

Kell, D. B. Drug Discovery Today 2006, 11, 1085–1092. Lindon, J. C.; Holmes, E.; Nicholson, J. K. FEBS J. 2007, 274, 1140–1151. Bedair, M.; Sumner, L. W. TrAC, Trends Anal. Chem. 2008, 27, 238–250. Issaq, H. J.; Abbott, E.; Veenstra, T. D. J. Sep. Sci. 2008, 31, 1936–1947. Pierce, K. M.; Hope, J. L.; Hoggard, J. C.; Synovec, R. E. Talanta 2006, 70, 797–804. 10.1021/ac801012y CCC: $40.75  2008 American Chemical Society Published on Web 10/28/2008

for plant metabolomics as well as studying bacteria and animals.2,4-8 Pattern recognition of complex samples benefits from utilizing two-way data generated by hyphenated techniques like GC/MS and LC-MS, as the separation step divides the total mass spectrum (obtained with direct infusion) in a number of time separated contributions. When studying urine samples, capillary electrophoresis (CE) has proven to have a number of advantages over other separation techniques, providing high resolving power of charged hydrophilic molecules with high sample throughput and low sample volume requirements. The ability of CE to handle crude samples with low memory effects (provided that the capillary is properly flushed between the runs) is also a benefit compared to, for example, liquid chromatography, where columns may clog, or gas chromatography, where the derivatization step is both time-consuming and likely to introduce variations to the analysis. For detection of the fast transient signals that are generated by an efficient CE separation system, time-of-flight mass spectrometry (TOF-MS) is well suited. In view of its merits, CE-TOFMS with electrospray ionization (ESI) in both positive and negative mode was the choice of technique for this study. Since successful exploring of metabolic fingerprints often demands replicates with respect to individuals as well as sampling occasions in order to increase the confidence, data evaluation tools that can handle large data sets are highly needed. There are several commercial chemometric tools available that can be used for data evaluation in metabolomic studies. Most of them rely on extracting features such as specific peaks and then comparing these features in order to define, and identify differences between, the groups. A problem with many of the commercial tools is, however, the limited amount of data that can be handled, and usually it is not possible to compare more than two samples (fingerprints) at a time. Furthermore many of them are expensive to purchase, even though there are examples of software that can be downloaded free of charge from the Internet.9,10 In the literature there are several reports on different ways of exploring sets of two-way metabolic fingerprints. A first obstacle when comparing data from different runs is the unwanted but inevitable shifts in retention or migration times. Different peak alignment strategies have been suggested, either by selection of landmarks in the fingerprints to be aligned11,12 or by optimizing the positions for a number of time segments in order to obtain best agreement.13-15 Although the data set forms a three-way structure (time × m/z × run) each two-way fingerprint is usually unfolded into rows with sequential mass spectra or, equivalently, time profiles for each mass channel. The whole data set is thus represented as a two-dimensional data table, which then may be (8) Cubbon, S.; Bradbury, T.; Wilson, J.; Thomas-Oates, J. Anal. Chem. 2007, 79, 8911–8918. (9) Katajamaa, M.; Miettinen, J.; Oresic, M. Bioinformatics 2006, 22, 634– 636. (10) Smith, C. A.; Want, E. J.; O’Maille, G.; Abagyan, R.; Siuzdak, G. Anal. Chem. 2006, 78, 779–787. (11) Malmquist, G.; Danielsson, R. J. Chromatogr., A 1994, 687, 71–88. (12) Krebs, M. D.; Tingley, R. D.; Zeskind, J. E.; Holmboe, M. E.; Kang, J.-M.; Davis, C. E. Chemom. Intell. Lab. Syst. 2006, 81, 74–81. (13) Nielsen, N.-P. V.; Carstensen, J. M.; Smedsgaard, J. J. Chromatogr., A 1998, 805, 17–35. (14) Eilers, P. H. C. Anal. Chem. 2004, 76, 404–411. (15) Pierce, K. M.; Wright, B. W.; Synovec, R. E. J. Chromatogr., A 2007, 1141, 106–116.

explored by multivariate data analysis tools like principal component analysis (PCA) and discriminant analysis based on partial least-squares regression (PLS-DA). A primary task for the data exploration is to evaluate similarities and dissimilarities between the samples or runs, i.e., identifying clusters and potential outliers. This procedure is usually guided by some beforehand-given information about the objects or samples, for example, class assignments according to clinical observations, object treatments, age and gender, etc. Principal components analysis has frequently been used as a tool to reduce the large number of variables (time and m/z combinations) into a much smaller set of score vectors, which then are subject to further evaluation with respect to the a priori information.7,8,16-18 Sometimes a visual inspection of the score plots, with the samples labeled according to class assignment, suffers to reveal which principal components (if any) are relevant for the phenomena at study. The loadings for the selected principal components can be used to pinpoint the variables involved, possibly facilitating identification of the underlying chemical components. With PLSDA the formation of scores and loadings is guided by the known class assignment. This supervised method, especially suited for binary classification, has also been applied in metabolic studies.7,8,18 An important issue in the multivariate analysis is the scaling of the variables, in order to enhance the influence of variations in low-abundant variables. Reported procedures have been standardization to unit variance, logarithmic conversion, and square rooting, and by different scaling the relevance of the variables will probably be judged differently. In some recent studies the individual variables were directly evaluated and selected with respect to their class-related information.19,20 The variable selection was performed according to a statistical measure of the class separation, the Fischer ratio, which is applicable also in case of more than two classes. With this procedure the variables are scaled according to the pooled within-group variance rather than the overall variance. The work presented in this paper is a part of on ongoing project aiming at a set of generic Matlab tools for exploration of two-way fingerprints of complex samples, obtained by CE/GC/LC separation with MS detection. The guiding star has been to use original raw data when possible and to perform rapid initial data analysis by visually guided and interactively controlled procedures. A challenge has been the limitations in computer performance about data storage and calculation speed, and much effort has been paid to efficient data reduction and representation. As an application for this study, urine samples were collected from a group of 13 individuals 2 h after intake of coffee, tea, or water at three separate occasions for each beverage. The variations between samples were due to two factors, beverage intake and individuals, with random “error” variations occurring between occasions. In total 115 samples were analyzed using CE-ESI-TOF(16) Ullsten, S.; Danielsson, R.; Backstrom, D.; Sjoberg, P.; Bergquist, J. J. Chromatogr., A 2006, 1117, 87–93. (17) Johannesson, N.; Olsson, L.; Ba¨ckstro ¨m, D.; Wetterhall, M.; Danielsson, R.; Bergquist, J. Electrophoresis 2007, 28, 1435–1443. (18) Gika, H. G.; Theodoridis, G. A.; Wilson, I. D. J. Sep. Sci. 2008, 31, 1598– 1608. (19) Pierce, K. M.; Hoggard, J. C.; Hope, J. L.; Rainey, P. M.; Hoofnagle, A. N.; Jack, R. M.; Wright, B. W.; Synovec, R. E. Anal. Chem. 2006, 78, 5068– 5075. (20) Johnson, K. J.; Synovec, R. E. Chemom. Intell. Lab. Syst. 2002, 60, 225– 237.

Analytical Chemistry, Vol. 80, No. 23, December 1, 2008

8947

MS in both positive and negative mode. Following the steps outlined in previous works,16,17,21-23 a “fuzzy” correlation matrix (115 × 115) was obtained for each mode. The PCA scores obtained directly from the correlation matrix were subject to two-way analysis of variance (ANOVA) based on a general linear model (GLM), and the most informative pair of components (according to the F value for beverage intake) was selected for score plots. When ANOVA indicates significant variation between groups, the next step is usually to evaluate pairwise differences either between two groups or between two group combinations. The calculation of such contrasts (following the GLM model) for each variable have been implemented as a way of identifying and presenting “hot spots” in the metabolic fingerprints. Such areas can be selected by mouse control for a “close-up” data analysis, with possible substance identification or at least as a guidance for refined instrumental analysis (e.g., MS/MS). The zoomed areas can also be selected to monitor substances known beforehand, and some attempts were made to verify other findings in the literature concerning caffeine metabolites24 and potential endogenous biomarkers for coffee consumption.25 In this work two-way GLM-ANOVA was implemented in order to enhance the effects of beverage intake in presence of large individual variations. This was possible due to the current design, where “individuals” was a crossed factor (all individuals tried all beverages). A more common situation is a nested design, where different groups of individuals are subject to different treatments. However, there are often other crossed factors like gender and/ or age, and also in such cases the data evaluation should benefit from the multiway GLM-ANOVA approach. EXPERIMENTAL SECTION Chemicals and Materials. Formic acid, acetic acid, ammonium acetate, 2-propanol (all pro analysi), sodium hydroxide (Titrisol-ampule) (Merck, Darmstadt, Germany), ammonium formate (>99.5% pure; BDH Laboratory supplies, Poole, U.K.) were used. Doubly distilled water was produced using a Milli-QPLUS water purification system (Millipore, Bedford, MA). Fused-silica capillary tubing was obtained from Polymicro Technologies Inc. (Phoenix, AZ). Urine Samples. Thirteen healthy individuals fasted for 12 h prior to intake of beverage consisting of coffee, tea, or water. Each individual collected urine 2 h after intake for 3 days for each beverage. Collected urine samples were directly frozen and stored at -18 °C. Before CE-MS analysis, urine samples were thawed at room temperature and gently shaken. Prior to injection each sample was diluted 5 times with MilliQ water. Instrumental Setup. Separations were performed using an Agilent Technologies HP3D CE instrument (Agilent Technologies, Waldbronn, Germany) coupled to an Agilent Technologies series 1100 LC/MSD TOF mass spectrometer (Agilent Technologies, Santa Clara, CA) equipped with a sheath flow ESI interface (21) Danielsson, R.; Backstrom, D.; Ullsten, S. Chemom. Intell. Lab. Syst. 2006, 84, 33–39. (22) Backstrom, D.; Moberg, M.; Sjoberg, P. J. R.; Bergquist, J.; Danielsson, R. J. Chromatogr., A 2007, 1171, 69–79. (23) Daszykowski, M.; Danielsson, R.; Walczak, B. J. Chromatogr., A 2008, 1192, 157–165. (24) Schneider, H.; Ma, L.; Glatt, H. J. Chromatogr., B 2003, 789, 227–237. (25) Ito, H.; Gonthier, M.-P.; Manach, C.; Morand, C.; Mennen, L.; Re´me´sy, C.; Scalbert, A. Br. J. Nutr. 2007, 94, 500–509.

8948

Analytical Chemistry, Vol. 80, No. 23, December 1, 2008

(G1603A Agilent Technologies). A 50 cm fused-silica capillary with 50 µm i.d., 365 µm o.d. was used for the separations. CE Conditions. For every 20 injections the capillary was removed from the electrospray interface and conditioned for 20 min with 1.0 M sodium hydroxide, 5 min with MilliQ water, and finally for 5 min with the separation buffer consisting of 50 mM ionic strength ammonia/ammonium acetate, pH 9. Injections were made hydrodynamically at 50 mbar for 5 s without turning off the nebulizing gas. Separation was carried out using a field strength of 420 V cm-1, with a total of 21 kV applied to the injection end. Between injections the capillary was conditioned for 1 min with separation buffer. The first 30 cm of the capillary was thermostatted to 25 °C, whereas the last 20 cm was left without temperature control at ambient temperature. MS Conditions. A sheath flow consisting of 2-propanol/buffer solution (pH 9), 80:20 (v/v) was added at a flow rate of 2 µL/min to the CE/MS interface by a Jasco PU-980 HPLC pump (Jasco Inc., Easto, MD) to aid the electrospray formation. Detection was carried out in both positive and negative ion mode using the following mass spectrometer parameters: capillary voltage -4 kV (in negative ion mode 4 kV), gas temperature 130 °C, drying gas 5 psi, nebulizer 2 psi, fragmentor 175 V, skimmer 80 V, OCT RFV 250 V. Data was collected between m/z 50 and 2500 u with 2000 transients per spectrum, resulting in 3.54 Hz. DATA PROCESSING According to the design, the number of CE-MS runs should have been 13 × 3 × 3 × 2 ) 234; for the 13 individuals, since there were three replicate occasions for each of the three beverages with MS analyses in both positive and negative mode. However, there was one missing sample for two of the beverages, resulting in four missing runs and the total number of 230. Due to the vast amount of data generated in this study, a great deal of effort was paid to procedures for data handling and data exploration. The main goal was to establish and describe possible systematic differences between samples, related to the intake of beverage. This section describes the different steps involved, in some places with reference to illustrative figures later discussed in detail in the Result and Discussion section. Collection and Preprocessing Raw Data. For each run the raw data (consecutive MS spectra during the CE separation) was stored in the proprietary WIFF format, later converted to NetCDF26 by an Analyst software utility. The analyses within this study comprise 230 such files, each holding about 350 MB of data (50-100 million nonzero data points per run). Each NetCDF file was read into Matlab,27 where the data was organized in a regular mesh with respect to m/z and migration time, respectively. The mesh was defined by one vector of m/z values in steps of 0.01 u and one vector of time values for the consecutive spectra (about 0.3 s per spectrum). The data points were stored in a Matlab data file as a list of intensities (s) together with corresponding lists of m/z and time indices (i and j, respectively). In addition, the mesh vectors (masses and spectra) and identification data were stored in the data file. In the next step, the data was preprocessed run by run in order to remove background, noise, and spikes as described previ(26) NetCDF. http://www.unidata.ucar.edu/software/netcdf/ (July 2008). (27) Matlab. http://www.mathworks.com/ (July 2008).

Figure 1. Metabolic fingerprint for one of the samples (after coffee intake). The lower left panel shows the fingerprint for the total time and m/z region. The small rectangle, visible as a slightly thicker dash, indicates the zoomed area in the upper left panel. The intensities in the zoomed area are shown as a summed mass spectrum (lower right panel) and a time profile (upper right panel).

ously.16 The amount of data was hereby substantially reduced, to about 1-5 million data points per run, and the runs were collected in two data sets (positive and negative ESI mode, respectively). Each data set (115 runs) was stored in a single file, this time in the HDF5 format that allows Matlab to write and read the results from individual runs.28 The different data blocks are addressed by paths similar to ordinary file addressing. Finally, a coarse linear time alignment was performed for all runs within each data set. This manual procedure was guided by the correlation plot between a selected master run and the run to be time adjusted, as recently described in ref 22. Starting with the NetCDF files, the procedures for data collection and preprocessing were performed with in-housedeveloped Matlab routines, run on a dual Opteron computer (2.0 GHz clock frequency, 6 GB RAM). For each run the data processing time was almost the same as that for the CE-MS analysis, i.e., about 10 min. Compiling Data Sets. So far, the procedures were performed for each run separately, which made it possible to handle all the data generated during the run. Even though the baseline subtraction and noise elimination results in considerable data reduction, it was impracticable to apply the exploratory tools directly on the whole data set. There are different ways to reduce the amount of data, by (i) selecting a subset of runs, (ii) restricting the m/z and time range, and (iii) decreasing the m/z and time resolution (binning or bucketing). In this study all runs (115 for each of the two data sets) were used, while the m/z range was confined to 50-500 u with resolution 0.1 u and the time range to 50-350 s with resolution 2 s. By this the total number of nonzero data points was reduced to 14.1 million and 9.2 million in positive and negative ESI mode, respectively. For each run the data were normalized to unit sum, thus neutralizing variations in, for example, urine concentration and ESI efficiency. When compiled, the two data sets have an internal three-way structure (samples × m/z × time). However, each data set was

represented as an ordinary two-way data matrix (denoted X1) in Matlab by unfolding each run (m/z × time) into a single vector, forming one row of the matrix. When loading the data matrices, the amount of data was further reduced by deleting data points with intensities below a certain threshold. In the following this threshold was set to 0.001 times the maximum intensity found in the data set, which reduced the number of nonzero data points to about 90 000 in positive and 110 000 in negative ESI mode. The individual fingerprints can be presented and visually examined as shown in Figure 1. No peak-picking algorithm was applied, but roughly 100-200 peaks corresponding to clusters of 3-10 data points could be observed in Figure 1. Correlations or Similarities between Fingerprints. The primary task for the data exploration was to elucidate the relation between the samples and to investigate the interrelation with the known experimental factors (individual and beverage). The relation, or rather correlation, between the runs can be found by multiplying the data set matrix X1 with its transpose X1T. The result is a correlation or similarity matrix, which can be normalized to hold values between 0 (no correlation) and 1 (full correlation). Prior to the calculations, however, the dynamic range of the intensities was compressed in order to reduce the influence of high-intensity peaks on the correlation measure. The intensities (s) were transformed according to s(1/n) where n is the compression parameter. With n ) 1 the intensities are unchanged, whereas n ) 2 results in square rooting. Increasing n leads to higher compression, in the limit resulting in binary correlation where the results only depend on the presence or absence of intensities at corresponding positions. A positive side effect of data compression is the reduced influence of “closure” caused by the normalization to unit sum of intensities (or squared intensities).22 In this work the compression parameter was set to n ) 3. Another obstacle in getting informative correlations is that the effective CE time scale varied from run to run, even though the coarse prealignment and the time binning greatly reduced this effect. Analytical Chemistry, Vol. 80, No. 23, December 1, 2008

8949

Figure 2. Total ion counts (all samples) for positive (left panel) and negative (right panel) ionization mode. The migration time is divided in three regions with respect to the migrating species: (i) positive ions, (ii) neutrals, and (iii) negative ions.

To deal with the alignment problem we applied the “fuzzy” correlation concept introduced in previous papers.17,21 In short, the X1 matrix is multiplied with the transpose of a “blurred” version X2, where the individual time profiles (for each m/z channel) are “blurred” by a moving average filter. The average is taken over a certain number, the fuzziness parameter, of time bins on each side of the current time bin. The fuzziness parameter should be selected to cover the time shifts occurring between the runs within the data set. Because even minor time shifts could result in varying peak distribution between adjacent time bins, a minimum value of 1 should be used. The same holds true for possible shifts in the m/z directions, and the individual mass spectra are “blurred” in the same way. The “blurring” procedure implies that when a certain run (i.e., row in X1) is correlated (multiplied) with another run (i.e., row in X2) the cross product will include the product of the two peak values originating from the same ion (i.e., compound). This occurs in spite of possible shifts in time and m/z within the range defined by the fuzziness parameters. The cost for the reduction of detrimental shift effects is a loss of resolution; other peaks occurring in the same fuzziness region, as defined by the bin sizes and fuzziness parameters, will give contributions to the correlation. Thus, the original peak resolution in time (for the individual electropherograms) and m/z (for the individual mass spectra) will put a limit to the amount of shift that can be handled. However, with the high mass resolution obtained with TOF mass spectrometers there are seldom multiple peaks in the individual electropherograms, and less so within the same “blurring” interval. The time and m/z fuzziness parameters were 5 and 1, respectively, resulting in an effective time resolution of 22 s compared to the investigated 300 s migration window. Visualization of the Similarities by Score Plots, Selected by Means of ANOVA. A table (matrix) S of all mutual similarities between the runs is obtained after normalization of the matrix product X1(X2)T. The eigenvectors of S, scaled by the square root of the corresponding eigenvalue, provide score values in the same sense as PCA. When based only on the similarities, and not directly on the original variables, one may think of PCA as a way to locate the runs in a low-dimensional space in accordance with the given similarities. With ordinary PCA, however, the similarities are calculated from “unblurred” data alone, namely, as X1(X1)T, with the adherent need of peak alignment. Another divergence from ordinary PCA is that, in our case, the underlying 8950

Analytical Chemistry, Vol. 80, No. 23, December 1, 2008

variables (the columns in X1 and X2) are noncentered. This implies that our first eigenvector has quite a different significance than what usually is connected to the first principal component. Since we use the singular value decomposition (SVD) procedure in Matlab, an eigenvector will be referred to as a singular vector (SV) rather than a principal component (PC). In the principal component context, SV2 corresponds to PC1, SV3 to PC2, etc., while SV1 has its own merits. A score plot where SV1 is involved has proved to serve as a good tool for identifying outliers, i.e., runs with low similarities to the other runs. However, no such outliers were found in this study. A score plot involves selection of two score vectors, and when looking for evidence of systematic differences due to a selected factor (in our case “beverage” or “individual”) it may be a timeconsuming task to find the most informative pair. Even though we restrict the number of score vectors to 10 there will be 45 different pairs to choose from, and as a guide for this selection we apply a two-way ANOVA to each score vector, with type of beverage and individuals as factors. In our case it is the influence of the “beverage” factor that is of primary importance, and for each score vector this influence is quantified as the F ratio or the corresponding p value. A high F ratio implies an appreciable variation due to the type of beverage, further corroborated by a sufficiently low p value. Hence the score plot obtained by the two score vectors with the highest F values will yield the most convincing evidence of systematic “beverage” effect on the CE-MS fingerprint. Visual inspection can also reveal clustering according to the type of beverage if the different runs (“fingerprints”) are labeled accordingly. The ANOVA modeling of the score vectors, preferably as a GLM due to possible imbalance in the design, also enables the computation of least-squares means (LSM) for the factor levels (different beverages as well as different individuals). The true beverage influence may be hard to discern in the full score plot including all the runs, as variations due to different individuals and different sampling occasions also give rise to some scattering along the selected score vectors (illustrated in Figure 3, upper panels). Taking the means for each individual and beverage over different sampling occasions, and subtracting the average “individual” contribution over all beverages, results in less scattered score values. Such score plots, with one point for each combination of “beverage” and “individual”, is a kind of

Figure 3. Score plots for the two most prominent singular vectors (SV) according to ANOVA, for positive (left panel) and negative (right panel) ionization mode. The upper row shows score values for all samples, while the middle row mean shows scores for each individual and beverage (adjusted for individual means). The lower row shows the least-squares means (LSM) for each beverage as Studentized t values, with confidence regions based on the least significant difference (LSD).

interaction plot showing the beverage effect for each individual. Systematic differences between the different beverages will appear more clearly in this type of presentation (cf., Figure 3, middle panels). Finally, the means can be taken for each beverage over all individuals as well, resulting in the average deviation (from the overall mean) for each beverage. These means can further be divided with the associated standard deviation, resulting in t values in conformance with the Student’s t test. When comparing the means for two beverages, the least significant difference (LSD) can be used as a criterion. In a two-dimensional plot with the three mean points for each beverage, ellipsoid regions around each point

based on the LSD measures indicates which beverages that can be viewed as significantly different with respect to the fingerprints (cf., Figure 3, lower panels). The conventional statistical significance, however, should not be stressed. When judging the p values obtained by ANOVA for each score vector, we should require very low p values (far below the usual 5% or 1%, or even 1‰) in order to rank a score vector as relevant. Initially the score vectors were transferred to Minitab29 where the ANOVA calculations were performed using GLM procedures. (28) HDF5. http://www.hdfgroup.org/HDF5/index.html (July 2008). (29) Minitab. http://www.minitab.com/ (July 2008).

Analytical Chemistry, Vol. 80, No. 23, December 1, 2008

8951

Figure 4. Differential fingerprint (positive contrasts) between coffee and tea or water in positive ESI mode, displayed and explored in the same manner as the individual fingerprints (cf., Figure 1) The right-most panel shows the contributions from the samples (each individual for each beverage) to the summed contrasts in the zoomed area (14.5% of all contrast values). The contrast centroid is found at m/z 94.03.

Later these procedures were implemented in Matlab and validated by comparison with results obtained with Minitab. Differential Fingerprints (Contrasts). Although score plots and scorewise ANOVA may indicate distinguished systematic “fingerprint” differences due to the intake of different beverages, we have no clues to what peaks (ions) that constitute these differences. In ordinary PCA the next step would be to investigate the loadings for the selected PCs. In previous papers16,25 we have indicated how such loadings could be obtained by correlation of X2 with the score vectors. However, if we are looking for variables with variation in conformance with a priori knowledge (here type of beverage taken) we could directly correlate X2 with an appropriate target vector. This approach is well suited for the case when two levels (beverages) are to be compared. With more than two levels we have to make pairwise comparisons between levels and/or groups of levels. For a balanced design the target vector can be constructed in a rather intuitive way, as demonstrated in recent papers.17,22 Correlating the target vector with the data vectors (variables) then corresponds to calculate the difference between the two means for the two levels (or group of levels). For a large study, however, there are often missing samples or samples that must be excluded as outliers. In such cases it is far from obvious how to calculate the various means, but again the GLM approach provides a simple solution. The target vector t is composed of the coefficients used to calculate the so-called “contrast” between two levels (or groups of level), which is defined as the difference between the corresponding mean values (as LSM). Hence the matrix multiplication (X2)tT results in a contrast vector that may be refolded into an m/z × time matrix. This matrix constitutes a differential “finger8952

Analytical Chemistry, Vol. 80, No. 23, December 1, 2008

print” with respect to the two levels (or groups of levels) under comparison (presented in Figure 4). Due to the large number of intensity variables we could expect that rather high contrast values (positive as well as negative) may be obtained by mere chance, even with only random intensity variations. In order to find a suitable limit for sorting out contrasts of potential interest, the target vector is randomly rearranged into 100 “random targets”. All the corresponding contrasts, 100 for each intensity variable, are compiled, and the 5% and 95% percentiles are selected as negative and positive thresholds, respectively. Only contrasts outside these limits are retained, resulting in a more distinct differential fingerprint. The contrasts are also divided into positive and negative contrasts, and from a sorted list the most prominent variables can be identified. The m/z position for these “hot spots” may be important clues in, for example, the evaluation of metabolic pathways or the search for biomarkers. Exploring Fingerprint Details by ANOVA. When “hot spots” or interesting fingerprint areas have been identified from the contrasts (or defined by a priori information), there is a need to turn back to the intensity values (i.e., the fingerprints) and relate them to the factor levels in a more direct way. Most efficiently this is performed for selected columns in X2 rather than X1, as the “blurring” ensures that there are columns with the relevant information in spite of possible shifts in time (and m/z). Moreover, for this procedure it might be better to use a version of X2 based on uncompressed data. The options we have chosen are that a fingerprint area, sized from the full fingerprint down to a single matrix cell, is selected interactively by a graphical tool. For each fingerprint (run) we obtain a total intensity, summed over time as well as m/z, as a new variable. This variable is analyzed by

Table 1. Extension of the Three Retention Time Regions for Positive Ions, Neutrals, and Negative Ions, Respectively, Are Given Together with the Relative Distribution of the Number of Data Pointsa

positive ESI mode

retention time region distribution of data points beverage separation negative ESI mode retention time region distribution of data points beverage separation a

positive ions

neutrals

negative ions

50-145 s 78.7% F ) 33.2 (SV7) F ) 18.2 (SV5) 50-165 s 53.3% F ) 5.0 (SV5) F ) 3.0 (SV3)

145-165 s 1.5% F ) 8.9 (SV10) F ) 5.8 (SV9) 165-220 s 19.0% F ) 4.5 (SV2) F ) 3.7 (SV8)

165-350 s 19.8% F ) 77.3 (SV6) F ) 14.3 (SV7) 220-350 s 27.7% F ) 318.3 (SV4) F ) 3.8 (SV6)

As a measure of the separation between the beverages, the F ratios for the two most prominent singular vectors (SV) are presented.

ANOVA in the same way as previously the scores; with F and p values, t values for the three beverage means, and means over the replicates (different occasions). This partly graphical information will assist in the critical evaluation of the relation between the beverage intake and the selected fingerprint detail, for example, an intense peak, the electroosmotic flow (EOF) migration band, or a mass range. RESULTS AND DISCUSSION Figure 1 demonstrates the graphical tool to present the metabolic fingerprint, here for one of the samples after intake of coffee, obtained in positive ionization mode. The total fingerprint is shown as a matrix image in the lower left panel. The cells in the data matrix are displayed as rectangular fields; for clarity a number of the original bins in m/z and time are combined to provide fields of visible size in the figure. The resulting intensity is indicated by the field color, here according to a grayscale. The marked rectangle, selected interactively, defines a subsection that is shown in the zoomed window (upper left). The intensities within that window are summarized over the m/z bins, displayed as a migration time profile (upper right), as well as over the time bins, displayed as a mass spectrum (lower right). The most intense peak in the selected area represents 4.3% of the total intensity for the run, and the center of gravity along the m/z axis (the centroid) is 114.26 u. Overview of the Collected Data. The intensities can be summed over all m/z bins in all runs, which results in a migration time profile that represents the whole data set. Figure 2 presents these profiles for negative (left) as well as positive (right) ESI mode. By visual inspection of the elution profiles, three distinct regions were defined in relation to the EOF. Neutral (and nearly neutral) molecules are eluted in the small time segment around the EOF position, being the cause of the sharp high peak in the middle of the time profiles. The left-hand region covers positive ions in the solution, migrating before the EOF, while the righthand region covers negative ions migrating after the EOF. The distribution of the intensities over these regions is summarized in Table 1. SVD Scores and ANOVA. With the use of the full m/z × time window, the correlation (similarity) matrix was computed and subject to SVD analysis. The occurrence of metabolomic patterns related to beverage intake was elucidated by ANOVA on the SV scores. Strong evidence for such patterns was found, for positive as well as negative ionization mode. The ANOVA results for the first 10 SVs are shown in Table 2. With positive ESI mode, the strongest group separation is seen in SV7 and SV4 (F ratios of 23.2 and 21.0, respectively), with some beverage influence also

Table 2. Top 10 Singular Vector (SV) Components for Both Positive and Negative ESI Mode for the Whole Retention Time Region (50-350 s)a positive ESI mode

negative ESI mode

SV no.

F ratio

p value

F ratio

1 2 3 4 5 6 7 8 9 10

0.2 1.4 1.0 21.0 0.6 3.4 23.2 9.7 6.6 0.3

0.85 0.25 0.38