Accepted Manuscript Title: A Method for Comparative Metabolomics in Urine Using High Resolution Mass Spectrometry Author: Padma Ramakrishnan Sreenath Nair Kannan Rangiah PII: DOI: Reference:
S0021-9673(16)30225-4 http://dx.doi.org/doi:10.1016/j.chroma.2016.02.080 CHROMA 357354
To appear in:
Journal of Chromatography A
Received date: Revised date: Accepted date:
7-1-2016 11-2-2016 29-2-2016
Please cite this article as: Padma Ramakrishnan, Sreenath Nair, Kannan Rangiah, A Method for Comparative Metabolomics in Urine Using High Resolution Mass Spectrometry, Journal of Chromatography A http://dx.doi.org/10.1016/j.chroma.2016.02.080 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A Method for Comparative Metabolomics in Urine Using High Resolution Mass Spectrometry Padma Ramakrishnan, Sreenath Nair, Kannan Rangiah*
[email protected] Metabolomics Facility, Centre for Cellular and Molecular Platforms, National Centre for Biological Sciences, GKVK, Bellary Road, Bangalore-560065, India. *
Corresponding author at: Metabolomics Facility, Centre for Cellular and Molecular Platforms,
GKVK, Bellary Road, Bangalore-560065, India. Fax: + 91-80-2363 6662.
1
Highlights
A workflow to perform comparative metabolomics in urine sample was established
Two different chromatographic methods in two ionization modes were used
Spiked reserpine or taurocholate were used for data extraction and normalization
Statistical analysis showed clustering of the two groups and the QCs
Comparative metabolomics analysis can provide meaningful insights for clinical diagnosis
2
Abstract Developing a workflow for metabolite profiling from biological fluids using mass spectrometry is imperative to extract accurate information. In this study, urine samples from smokers (n=10) and nonsmokers (n=10) were analyzed using an ultrahigh performance liquid chromatographyhigh resolution mass spectrometry (UHPLC-HRMS) system. For the analysis, two different chromatographic methods [Reversed phase liquid chromatography (RPLC) and Hydrophilic interaction liquid chromatography (HILIC)], in two ionization modes (positive and negative) were used. Spiked reserpine (positive ion mode) or taurocholate (negative ion mode) were used for data extraction and normalization. Quality controls (QCs), prepared by pooling urine samples from both smokers and non-smokers (each n=10), were used to assess the reproducibility of the method. The final data output from SIEVE 2.2 after applying a cut-off for QC coefficient of variation (CV) <20%) and p-value (<0.05) showed 165, 83, 177 and 100 unique components in RP positive/negative, HILIC positive/negative modes, respectively. Statistical analysis showed clustering of the two groups and the QCs, while the variable importance in projection (VIP) scores for the top fifteen metabolites in each of the four modes indicated the metabolites most responsible for the differences. Application of the developed workflow for comparative metabolomic analysis of urine in different diseased models will be of great use in the field of clinical metabolomics.
Abbreviations UHPLC-HRMS:
ultrahigh
performance
liquid
spectrometry RPLC: Reversed phase liquid chromatography HILIC: Hydrophilic interaction liquid chromatography QCs: Quality controls VIP: variable importance in projection HMBD: Human Metabolome Data Base KEGG: Kyoto Encyclopedia of Genes and Genomes 3
chromatography-high
resolution
mass
PCA: Principal component analysis PLS-DA: partial least squares discriminant analysis SPE: Solid Phase Extraction ROC: receiver operating characteristic AUC: area under the curve CV: Coefficient of Variation.
Keywords: Metabolite profiling; urine; mass spectrometry; HRMS; data analysis; data normalization.
4
1. Introduction Metabolomics is a study of the detection and quantification of small molecules (metabolome), which are the functional readout of cellular phenotype and the end products of upstream gene expression [1]. It is a cutting-edge biochemical approach at the interface of metabolic phenotype and genotype [2]. Metabolite changes have long been observed in diseased individuals either as a primary cause or a secondary indicator. The concept that individuals might have a metabolic profile that could be reflected in the makeup of their biological fluids is well known [3]. Until recently, the quantification of all possible metabolites from the biological fluids has been a challenging task. New developments in high-resolution mass spectrometers (HRMS) in the past few years have enabled simultaneous detection of metabolites from biological fluids [4]. There are a growing numbers of studies using liquid chromatography-mass spectrometry (LC-MS) as a tool for biomarker discovery. The high correlation between metabolites and phenotype has created a surge of interest that is reflected in the increasing number of publications in the field of metabolomics. In the last 4-5 years LC-MS based non-targeted metabolomics has been used for clinical applications mainly as a tool for disease diagnosis [5-8].
Metabolomics has its own challenges in terms of analyzing different chemical classes of molecules and quantifying them all together in a single method. As of now analysis of metabolites in a single method is impossible, but the two major chromatographic techniques namely
reversed
phase
chromatography
(RPC)
and
hydrophilic
interaction
liquid
chromatography (HILIC) are known to cover most of the metabolites in a given biological fluid [9]. There are reports using HRMS with two chromatographic columns for analysis of sera, urine and saliva samples to demonstrate the changes in metabolome of normal versus disease samples [9-11]. Producing data from this kind of analysis is comparatively easy, whereas extraction of information from the produced data is still a challenging problem. Metabolomics data analysis relies heavily on advances in bioinformatics tools that are required for analysis and on electronic databases for metabolite identification [12]. Currently, a number of databases are available like Human Metabolome Data Base (HMDB) [13-14], Kyoto Encyclopedia of Genes and Genomes (KEGG) [15], METLIN [16], PubChem [17], MassBank [18] etc. Yet, annotation of all
metabolites is quite cumbersome and an on-going process. There are many freely available (XCMS, MZmine) and vendor specific (SIEVE-Thermo Fisher, MarkerView-AB Sciex, 5
MassLynx-Waters and Mass Profiler-Agilent) software for data extraction and analysis, but problems still exist in extracting complete metabolite information. Also, normalization of data is an issue much debated currently. This makes metabolomic analysis of biological fluids a highly complicated process. Integration of “all omics” (genomics, epigenomics, transcriptomics, proteomics, glycomics, metabolomics, microbiomics and phenomics) will be highly useful in constructing molecular networks [19], helping us understand the complex biochemical processes in a better way, but analyzing the datasets all together is a great challenge. It is therefore, necessary to establish standardized protocols of metabolomic analysis from different biological fluids for all researchers in the field of metabolomics. There is a great need to improve the data quality and data mining strategies to obtain more and accurate information.
Here we have shown a workflow for urine metabolite profiling by using UHPLC-HRMS system. To demonstrate, we have taken urine samples from smokers/nonsmokers and profiled metabolites. Smokers/nonsmokers were chosen as subjects because the differentiation is well established and it seemed ideal to validate a method by obtaining expected results. Smoking is known to cause large and growing number of premature deaths, and also has direct correlation to both lung and oral cancer [20]. The advantage with urine is that, it can be collected noninvasively, is available in large quantities, has a different metabolome from blood, and sample handling is simple since there is no need to remove protein [21]. It is also widely used as a diagnostic tool in clinical practice [5-8]. One of the major considerations is the presence of large amount of salts, which can interfere with MS analysis. Care should be taken to desalt in the sample preparation step. We have used RPC and HILIC coupled with both positive and negative ionization modes for analysis. To overcome the issues in data normalization we have used reserpine (positive ion mode) or taurocholate (negative ion mode) spiked along with the samples and QCs. QCs were prepared by pooling urine samples (n=20) from both smokers and nonsmokers as previously shown [22]. Data analysis was done using SIEVE 2.2 software in a three step process (chromatogram alignment, component extraction and ID generation via ChemSpider by searching both HMDB and KEGG databases). A number of tobacco-related metabolites were correctly picked as differentiating components between the two sample sets.
2. Materials and Methods 6
2.1.
Materials
Reserpine, taurocholic acid (sodium salt), ammonium acetate, acetone and formic acid were procured from Sigma-Aldrich (Bangalore, India). High purity MS grade solvents (methanol and acetonitrile) were obtained from Merck Millipore (Merck Millipore India Pvt. Ltd., Bangalore). Double-distilled water for LC-MS was obtained from our in-house distillation unit. Solid phase cartridges (Strata-X 33µ Polymeric Reversed Phase; 1mL, 30 mg) were obtained from Phenomenex, Inc. (Hyderabad, India). Urine samples (20-40ml) were collected from 10 healthy non-smoker and 10 healthy smoker male volunteers, aged 22-32 years. The smokers used on an average, 10-15 cigarettes per day and all samples were collected during the day time. All samples were stored at -20 °C till further analysis.
2.2.
Sample Preparation
Before analysis, all samples were thawed in ice, vortexed well and centrifuged at 14000 rpm for 10 min. Equal volumes of the supernatant from all 20 samples were pooled to prepare at least 5 QC samples. From the remaining supernatant of each sample and from each of the pooled QC sample, 800 µL was used for RP (positive and negative) and 200 µL was used for HILIC (positive and negative) analysis.
2.2.1. RPC mode Briefly, 800 µL of urine sample was acidified with 1 µL formic acid and centrifuged (10000 rpm, 5 min). The supernatant was cleaned using RP-SPE cartridges. Prior to loading, the SPE cartridge was conditioned with 1mL methanol followed by 1mL of water. The acidified urine was loaded onto the cartridge and allowed to bind to the column with gravity flow. It was then washed with 1 mL of acidified water (0.1% FA) twice. Metabolites were eluted using 1 mL of acetonitrile:methanol:water (40:10:1) mixture. It was then dried under vacuum and reconstituted with 80 µL of 25% methanol and used for both positive and negative ion mode analysis. For RPC positive 10 µL of reserpine (5 µg/mL), and for RPC negative 10 µL of taurocholate (25 µg/mL) were spiked to 40 µL each of the reconstituted sample and 10 µL was injected into the UHPLC-HRMS for analysis.
2.2.2. HILIC mode 7
Briefly, 800 µL of cold acetonitrile was added to 200 µL of each sample or QC, separately. Samples were vortexed and maintained at 4 °C for 20-30 min and centrifuged (14000 rpm for 5 min). From the resulting supernatant of each sample or QC, 200 µL was taken for HILIC positive analysis and 200 µL was taken for HILIC negative analysis. The HILIC positive samples were each spiked with 10 µL of reserpine (5 µg/mL) standard, while the HILIC negative samples were each spiked with 10 µL of taurocholate (25 µg/mL) standard. Both sets of samples were dried under vacuum and reconstituted in 50 µL of 80% acetonitrile and transferred to auto-sampler vials. For the analysis 10 µL was injected into the UHPLC-HRMS system.
2.3.
UHPLC-MS
The mass spectrometer used for the metabolite analysis is a Q-Exactive Orbitrap (Thermo Fisher Scientific, San Jose, CA, USA) equipped with heated electrospray ionization (HESI) source. It also houses an HCD (higher-energy collision dissociation) cell for carrying out MSn experiments. The Q Exactive is coupled to a Dionex UltiMate 3000 UHPLC system (Thermo Fisher Scientific, San Jose, CA, USA). This system is provided with a column oven (set at 40 ºC), an auto-sampler and a thermo-controller (set at 4 °C). It uses an in-line split-loop injection design and is equipped with an external needle wash system (95% methanol) to ensure zero percent carry over problems. A total of four experiments were performed with sample analysis in RPC (positive and negative) and HILIC (positive and negative) modes. MS operating conditions for all four analyses were as follows: Spray voltage, +2500V (-2500V for negative mode); Capillary temperature, 280 °C; Vaporizer temperature, 320 °C, Sheath gas, 30 arbitrary units (40 for negative mode); and, Auxiliary gas, 10 arbitrary units. Injector settings were as follows: 0-2 min: waste, 2-45 min: load, 45-50 min: waste.
2.3.1. RPC mode For the RPC experiments, a Kinetex XB-C18 column (2.1x100 mm, 1.7 µm, 100 Å, Phenomenex India Pvt. Ltd.) was used. The mobile phase for RPC positive analysis was: Solvent A, Water (10mM Ammonium Acetate, 0.1% FA) and Solvent B, Acetonitrile (0.1% FA). The mobile phase for RPC negative analysis was: Solvent A, Water (10mM Ammonium Acetate) and Solvent B, Acetonitrile. The gradient was optimized for both modes to get maximum separation
8
(0 to 2 min: 0.2% B, 2 to 18 min: 0.2 to 20% B, 18 to 33 min: 20 to 60% B, 33 to 38 min: 60 to 95% B, 38 to 43 min: 95% B, 43.1 to 48 min: 0.2% B) at 300 L/min flow.
2.3.2. HILIC mode For the HILIC experiments, a Kinetex HILIC column (2.1x100 mm, 1.7 µm, 100 Å, Phenomenex India Pvt. Ltd.) was employed. The mobile phase for HILIC positive analysis was: Solvent A, water:acetonitrile (1:1, 10mM Ammonium Acetate, 0.1% FA) Acetonitrile
and Solvent B,
(0.1% FA). The mobile phase for HILIC negative analysis was: Solvent A,
water:acetonitrile (17:3, 10mM Ammonium Acetate) and Solvent B, acetonitrile. The following gradient gave optimal resolution and was used in both ionization modes (0-2 min:100% B, 2-15 min:100-90% B, 15-25 min:90-80% B, 25-30 min:80-75% B, 30-35 min:75-20% B, 35-40 min:20-0% B, 40-45 min:0% B, 45-45.1 min:0-100% B, 45.1-50 min:100% B) at 300 L/min flow.
2.4.
Data Acquisition
Before commencing each experiment, the mass measurement was externally calibrated according to the manufacturer‟s instructions. Data was acquired in the FS/DDS (Full scan/ Data-Dependent Scan) mode. The settings for all four LC-MS methods (RP/HILIC in positive/negative) were same, excepting the polarities. The instrument scanned in the mass range 70-1050 m/z and alternated between MS and MS/MS scans. Resolution for the full scan was set to 70000. The AGC (Automated Gain Control) target for the Orbitrap was 1x106 ions and a fill-time (Max ITMaximum Inject Time) of 80 millisec was allowed. Data was acquired in the centroid mode. For the data-dependent scans, resolution was set to 35000. AGC target was 1x105 ions and Max IT allowed was 200 millisec. The Top 5 (most intense) ions were subjected to MS/MS using normalized collision energy (NCE) of 35 eV. The isolation window was set to 4.0 m/z and dynamic exclusion time was kept at 15 sec (depending on the average chromatographic peak width) to afford greater coverage of ions during MS/MS. To aid greater MS/MS coverage, also, an exclusion list was provided with 100 most intense ions from the blank or mobile phase. The exclusion lists were different for each mode due to the differences in polarities and mobile phase composition. The DDS data was also acquired in the centroid mode. Each analytical batch began with two or more blanks, followed by the test samples interspersed with blanks (at least two 9
blanks for every five samples). QC samples were included in the batch randomly to monitor the method and instrument reproducibility.
2.5.
Data Analysis
For data analysis Xcalibur .RAW files from each of the four experiments were processed separately
using
SIEVE
2.2.
Details
of
the
software
can
be
found
http://www.thermoscientific.com/en/product/SIEVE-software-differential-analysis.html.
at Peak
alignment and component detection were performed over the retention time range 2-45 min for HILIC positive and negative experiments, and 2-35 min for RPC positive and negative experiments. First, the blank files were used for performing background subtraction. SIEVE then uses a proprietary algorithm called Recursive Base Peak Framing to put signals from all sample files, above the given intensity cut-off, with nearly the same retention times as well as masses in one frame. First, the most intense spectral data point in a sample is identified and framing begins. Next the remaining spectral data points are evaluated to identify one having the greatest intensity that falls within the frame and corresponds to a sample other than the sample corresponding to the spectral data point identified earlier. This is an iterative process similar to the way datadependent scanning operates in which the intensities of the MS peaks trigger MS/MS events. The width and length of each frame depends on user-defined parameters namely the retention time and m/z values. Depending on the chromatography, the maximum retention time shift allowed for alignment and framing was 1 min for the two HILIC modes and about 0.4 min for the two RP modes. The MZ width setting was 40 ppm for all four modes. The intensity cut-off was set to 2000000 units based on the average base-line intensities in the raw files. The assumed charge state was +1 or -1. Integrated intensities were calculated using the in-built period-peak detection algorithm. The default adducts‟ list for positive or negative mode were used for identifying frames which could be possible adducts of the main components. Finally, the integrated ChemSpider tool was used to match components obtained by performing the above steps, with metabolites in HMDB and KEGG. The mass tolerance for the database search was kept at 10 ppm for all modes except RPC positive (20 ppm). The obtained data was normalized with the intensity of the internal standard (reserpine or taurocholate) and then, filtered by coefficient of variation (CV) for the QC samples (<20%) and normalized p-values for the smoker samples (<0.05). The filtered results were exported to MS Excel. Each of the four lists was further 10
scrutinized manually to remove entries with very close masses (within 10 ppm) and close retention times. Possible isotopic masses (with near identical retention times) were also removed. These lists were finally used for statistical analysis.
2.6.
Statistical Analysis using MetaboAnalyst
The four resulting filtered peak lists, containing component identities (numbers corresponding to a specific m/z-retention time pair), sample identities and their normalized peak intensities were imported separately into MetaboAnalyst 3.0 for multivariate analysis [23]. Each dataset was subjected to log transformation and Pareto-scaling, after missing values (if any) were autoreplaced by a small default value, and then submitted for analysis. Principal component analysis (PCA) was used to visualize clustering of the two sample sets and the QCs. To get an overview of the difference between the two sample sets, a heatmap was generated for each of the four datasets, using Euclidean distance measure and Ward clustering algorithm [23]. Cross-validated partial least squares discriminant analysis (PLS-DA) was used to obtain a list of top 15 discriminating components, ranked according to their Variable Importance in Projection (VIP) scores, with VIP > 1.0 being considered relevant for group discrimination. Components from the top 15 list for HILIC positive data, which also showed a high correlation, were used to construct a Receiver Operating Characteristic (ROC) curve using the linear SVM algorithm. The same components were then evaluated in the ROC tester tool, for their ability to correctly classify unknown samples as smoker (class 1) or non-smoker (class 2). Owing to a limited sample size, the 5 QC samples, and 2 smoker (S5 and S10) and 2 non-smoker (NS5 and NS10) samples from the original datasets were taken as the „test set‟ while the remaining 8 smoker and 8 non-smoker samples were taken as the „training set‟. Final confirmation of a few tobacco metabolites was done using the advanced mass spectral database, mzCloud. For this, the respective mzML files were used as input for the mzCloud search. Precursor ions were matched exactly with that in SIEVE and a mass tolerance of ±10 ppm was allowed. Tolerance for normalized collision energy was set to ±5 eV.
3.
Results and Discussions
3.1.
Metabolite Profiling of Nonsmoker and Smoker Urine
11
In recent years metabolite profiling of biological fluids has become increasingly important for finding biomarkers that can be used to diagnose disease or provide warning at the preclinical stage. As such, there is no single bio-analytical method to analyze all metabolites from a given biological extract. However, combining different separation and ionization methods can maximize the coverage of metabolites. Both RPC and HILIC columns are commonly used for metabolome analysis and are known to cover a wide range of metabolites. We have also applied this approach to analyze the urine metabolome from smokers and nonsmokers. The typical workflow used for sample processing is shown in figure 1. We have standardized the urine volume and metabolite extraction as mentioned in the method. The mobile phase system for each mode was carefully optimized based on available standards (data not shown). We used mobile phase systems with neutral pH conditions for both the negative ion mode experiments and with acidic conditions for both the positive ion mode experiments. We have used RP SPE for desalting samples for RPC mode and acetonitrile (2:8, urine to acetonitrile v/v) for desalting and protein precipitation in samples for HILIC mode. The ratio of the volumes of urine to acetonitrile for HILIC analysis turned out to be critical. A ratio higher than 2:8 seemed to be insufficient for precipitating the salts, which was evident from drastically diminished signals of the internal standards (data not shown). Data acquisition settings for both RPC and HILIC experiments were identical excepting the ionization modes. Signals from the standards (reserpine/taurocholate) were used to track all sample processing steps involved.
3.2.
Internal Standard based Data Analysis
One of the biggest challenges in untargeted metabolite profiling, post data acquisition is normalization of data. Previous studies have shown normalization with intensities of constant background ions, with protein concentration, UV-based normalization, or by estimating a specific metabolite like creatinine concentration in case of urine etc., [24]. In one of the recent studies, normalization based on specific gravity of urine has been shown to improve recovery of meaningful information from urine metabolome analysis [25]. There is, however, no consensus on a particular method of normalization. In this study, we have spiked known amounts of internal standards (reserpine or taurocholate) and used it for normalization. Normally, these two compounds are used to check the ionization efficiency of the mass spectrometer in positive and negative ionization modes, respectively. The typical UHPLC-HRMS chromatograms for the 12
RPC positive and negative modes of nonsmoker, smoker and QC samples, and the corresponding extracted ion chromatograms for reserpine (retention time 28.3 min) and taurocholate (retention time 24.7 min) are shown in figure 2A and 2B. Similarly, the UHPLC-HRMS chromatograms for the HILIC positive and negative modes of nonsmoker, smoker and QC samples, and the corresponding extracted ion chromatograms for reserpine (retention time 12.2 min) and for taurocholate (retention time 12.4 min) are shown in figure 2C and 2D. The added advantage of using the internal standard was that we were able to use it to set the processing parameters in SIEVE. It not only gave us a clearer idea of the retention time shifts that occurred across the samples, but also an estimation of the mass accuracy in the data acquired. In our four datasets, the retention time shifts for the internal standards ranged from 0.4 min for the two RPC experiments to 1 min for the two HILIC experiments. Likewise, the maximum error in the masses of the internal standards across the 20 samples and QCs, for the four experiments ranged from 10 ppm to 30 ppm. These were taken as the minimum values for setting up the alignment and component extraction parameters. This helped in making sure that we obtained only one frame/feature corresponding to the internal standard‟s aligned mass and retention time. By extrapolation, it was reasonable to assume that most, if not all, of the other features would also be unique. The retention time shifts also indicated that the reproducibility of chromatograms in RPC is better compared to HILIC (Figure 3). At the end of the alignment and component extraction process there were 6650 components for RPC positive, 5090 for RPC negative, 4006 for HILIC positive and 3290 for HILIC negative mode. Another advantage of having the internal standard was that the error between its aligned mass from SIEVE and its theoretical mass in the metabolic database could be used to set a reasonable search criterion. Accordingly, a 10 ppm mass tolerance was allowed for search against the HMDB and KEGG databases, in all modes except in RPC positive, where the tolerance was set to 20 ppm. Finally we could use the internal standards for normalization of the intensities. The integrated intensities of the standards before and after normalization for all four modes are shown in figure 3.
3.3.
Data Filtering
A major problem in mass spectrometry based profiling data, difficult to avoid, is presence of false positives in the final output. Before statistical analysis, therefore, it is useful to filter the data to remove features that may not really be signals. This is where having QC samples is 13
useful. It is difficult to apply a stringent CV cut-off for the test samples (smoker or non-smoker) since these are biologically different, having varying amounts of a given metabolite. We have preferred to apply a CV cut-off to the QC samples. Ideally, the QC samples are representative of all the masses present in all the test samples, and each QC has essentially the same composition. Thus, only features that are consistently detected at similar levels in all of the QC samples should be considered as real signals. A 20% cut-off for CV of QC samples was chosen in accordance with the validation guidelines for routine analytical methods. An additional filter of p-value for smoker samples < 0.05 was used to obtain a meaningful input for statistical analysis. In the filtered lists, it was found that there remained certain entries, which were nearly same in mass and retention time within each of the four datasets. In fact, a few of the entries turned out to be isotopic masses. A manual survey of each of the four lists was therefore deemed necessary. Masses within 10 ppm of each other and also within 0.5 min (RPC) and 1 min (HILIC) retention were considered repetitive. Like-wise, masses with a unit difference and with near identical retention times were considered isotopic. After removing these repetitive and/or isotopic entries, we obtained finally 165 unique components in RPC positive, 83 in RPC negative, 177 in HILIC positive and 100 in HILIC negative experiments (Supplementary Table 1). The statistics were performed on these four lists separately. As an additional analysis, we wanted to know the total number of unique masses obtained from the four experiments. Accordingly, masses within 15 ppm of each other in the two ionization modes (of RPC and HILIC) and then amongst the two chromatographic modes were considered as one. After accounting for these, we obtained in total 482 unique masses from four experiments, of which, 275 had suggested ChemSpider IDs and 207 remained unknown. The complete list of components in each of the four modes is shown in supplementary table 2 to 5.
3.4. Hierarchical Clustering and Discriminant Analysis of Metabolite Data We have used the freely-available MetaboAnalyst 3.0 software to analyze the filtered dataset obtained from SIEVE 2.2. The normalized and filtered data was pre-treated by applying gLog transformation and Pareto scaling, which were necessary to make the data points more comparable. This way, two components/metabolites with large differences in intensities are scaled to similar levels to give them both equal weightage while assessing their discriminatory powers. The PCA plots showed separation of the two sample sets in each mode while showing 14
tight clustering of the QCs (data not shown). The same differences were also evident in the PLSDA score plots (Figure 4). Maximum separation was observed in the HILIC and RPC positive datasets. A slight overlap was visible in the HILIC negative set. Clustering of QCs indicated the robustness of each of the four methods. Non-clustering of QCs normally indicates changes in signal intensities, changes in peak shapes or drifts in retention times during the run of an analytical batch. The heatmap for each mode allowed a direct visualization of the components responsible for the differences (Figure 4). As is evident, most of the metabolites are higher in the smoker set compared to the nonsmoker set in all four modes of analyses. In the RPC and HILIC positive modes, a few metabolites are noticeably higher in nonsmokers compared to smokers. However, the fold-change in all these cases was much lower than observed for the metabolites higher in smokers (Supplementary Tables 2 to 5). Also, one of the nonsmokers (No: 5) and one of the smokers (No: 5) were visibly different from the other samples in their respective classes, showing a somewhat opposing pattern. This was confirmed with the random forest test, which also indicated nonsmoker-5 and smoker-5 as being outliers in all four datasets (supplementary Figure 1). This result was unexpected, but could be due to the possibility that non-smoker 5 was in fact a sparse smoker or that smoker 5 had not actually smoked in the recent past. The outliers are indicated with arrows on the PLS-DA score plots (Figure 4). The cross-validated PLS-DA model resulted in a VIP score plot of the top fifteen metabolites differentiating smokers from nonsmokers (Figure 5). In all four modes, the top 15 components were higher in smokers compared to nonsmokers. It has been known that for higher n value, a VIP score of 1.5 is considered to enable discrimination between two phenotypes. In our study of 10 samples for each dataset, the VIP scores of the top 15 metabolites ranged from 1.2 to 2.9 (for all four modes), a higher score signifying higher discriminatory power. The ChemSpider suggestions for many of these metabolites, especially in the HILIC positive and RPC positive modes were nicotinerelated metabolites. Also suggested were some compounds which are possibly tobacco-flavoring agents/additives and poly-aromatic hydrocarbons. The next step should ideally be confirmation of these metabolites with standards. In the absence of standards, we decided to perform a spectral match of the MS/MS spectra of the top 15 metabolites in each mode with the spectral library, mzCloud. This produced high-confidence hits for a few nicotine metabolites, namely, nicotine (98.4), trans-3-hydroxycotinine (97.2), and cotinine (94.3) (Supplementary Figure 2 to 4). The mzCloud database is relatively new and is by no means exhaustive at present. A number of other 15
nicotine metabolites/tobacco components were also suggested but not confirmed. These included Dibenzo[h,rst]pentaphene which is a polyaromatic hydrocarbon; and cotinine-N-oxide, hydroxynicotine and 4-(methylamino)-1-(3-pyridyl)-1-butanone which are all major nicotine metabolites.
3.5. Biomarker Analysis We next wanted to check the predictive efficacy of the significantly discriminating metabolites in the HILIC positive dataset. This set was chosen because it not only had the most number of differential metabolites (177), but also showed the most separation, as observed in the PLS-DA score plots (Figure 4). For this, the area under the ROC curve (AUC) was used, which is a simple statistical tool to assess the predictive value of discriminating metabolites. The ROC curves can be used individually for each candidate biomarker or even for a combination of them. An AUC of 1.0 indicates a perfect predictor while 0.5 indicates a random prediction. In our experiment, the ROC curve was plotted using linear SVM algorithm. Incidentally, each of the top 15 metabolites in our dataset had an AUC value of 1. Using all these for a combined ROC might give a very good AUC value, but it runs the risk of over-fitting the data. Therefore, metabolites or a combination of metabolites need to be selected based on prior knowledge for example based on whether they belong to a common metabolic pathway. As a general strategy, we chose to select those components in the significant top 15, which also showed a high correlation (high positive or negative correlation coefficients), for constructing the ROC curve. We therefore selected components highly correlating (Pearson‟s correlation coefficients ≥ 0.89) with component HP 1408 for the ROC analysis and obtained an AUC value of 1 (Figure 6). Additionally, the ROC tester tool was used to assess the classification power of these 15 metabolites. There are two requirements for validating a classifier namely, model selection and performance estimation. For this, the data is split into two sets: the training set, for building the model and the test set, for assessing its performance. Sometimes, a separate validation set is also affordable in cases where the sample size is large. The suggested ratio of the training: validation: test set is 50:25:25. In keeping with the above guidelines, we split our data into the training set: 8 non-smokers and 8 smokers (64%) and the test set: remaining 9 samples (36%). Also, the model should be chosen in a way that it provides the lowest error rate on the entire population. Therefore, the identified outliers NS5 and S5 were excluded from the training set and instead, 16
included in the test set. NS10 and S10 were included in the test set randomly while all the QCs were included in the test set for lack of availability of any other samples. MetaboAnalyst then uses random sub-sampling to cross-validate the model with the selected biomarkers (15, in our case) and classify the unknown. As a result, samples S5, S10, QC1, QC2, QC3, QC4 and QC5 were correctly classified as smokers (class 1) while NS5 and NS10 were correctly classified as nonsmokers (class 2). Indeed, the probabilities of classification of both S5 and NS5 were slightly lower which is consistent with these being outliers. The QCs being pooled samples eventually have all of these 15 metabolites present, as compared with nonsmoker samples, which would be characterized by their absence; hence their assignment to class 1 was appropriate (Supplementary table 6).
Such a strategy when applied to a larger dataset and when tested in different
populations can lend credence to the selection of the hypothesized biomarkers.
4. Conclusions The objective of our study was to establish a method for profiling metabolites in urine, for justifying the use of internal standards in the component extraction process and data normalization, and outlining a strategy for selecting possible markers that distinguish one sample set from another. The best approach, we observed, was to work with two sample sets that are already known to be different (smokers and nonsmokers in this case) and derive the same conclusion from our experiment. This was necessary because the field of metabolic profiling is still at a nascent stage with continuous up-gradation of software and databases being the order of the day, questioning the reliability of the data produced. Our study reinstates confidence in the data generated by following such an approach. We have developed here a simple workflow for metabolite profiling from urine samples, using spiked standards which act as a control for the chromatography and are also used for data normalization. The UHPLC-HRMS data generated between smoker and nonsmoker urine samples indicated that many metabolites are higher in smokers compared to nonsmokers. As expected, most of these were putatively identified as being tobacco-related. A similar approach can also be applied to other disease models to understand the changes manifested in the metabolome.
17
Acknowledgement We would like to thank the Department of Biotechnology (DBT), Government of India, for financial assistance in establishing the Metabolomics Facility at C-CAMP, National Centre for Biological Science Campus, Bangalore, India.
18
References [1]
R. Goodacre, S. Vaidyanathan, W.B. Dunn, G.G. Harrigan, D.B. Kell, Metabolomics by numbers: acquiring and understanding global metabolite data, Trends Biotechnol. 22 (2004) 245–252.
[2]
C. Schmidt, Metabolomics takes its place as latest up-and-coming “omic” science, Jnci. 96 (2004) 732–734.
[3]
J.K. Nicholson, J.C. Lindon, Metabonomics, 455 (2008) 1054–1056.
[4]
X. Liu, Z. Ser, J.W. Locasale, Development and quantitative evaluation of a high resolution metabolomics technology, Anal. Chem. 86 (2014) 2175-2184.
[5]
B.M. Wittmann, S.M. Stirdivant, M.W. Mitchell, J.E. Wulff, J.E. McDunn, Z. Li, et al., Bladder Cancer Biomarker Discovery Using Global Metabolomic Profiling of Urine, PLoS One. 9 (2014) e115870.
[6]
H. Luan, L. Liu, Z. Tang, M. Zhang, K. Chua, J. Song, Comprehensive urinary metabolomic profiling and identification of potential noninvasive marker for idiopathic Parkinson‟s disease, Nat. Publ. Gr. (2015) 1–11.
[7]
G.F. Giskeødegård, S.K. Davies, V.L. Revell, H. Keun, D.J. Skene, Diurnal rhythms in the human urine metabolome during sleep and total sleep deprivation, Sci. Rep. 5 (2015) 14843.
[8]
X. Wang, a. Zhang, Y. Han, P. Wang, H. Sun, G. Song, et al., Urine Metabolomics Analysis for Biomarker Discovery and Detection of Jaundice Syndrome in Patients With Liver Disease, Mol. Cell. Proteomics. 11 (2012) 370–380.
[9]
Q. Wang, P. Gao, X. Wang, Y. Duan, The early diagnosis and monitoring of squamous cell carcinoma via saliva metabolomics, Sci. Rep. 4 (2014) 6802.
[10]
T. Zhang, D.J. Creek, M.P. Barrett, G. Blackburn, D.G. Watson, Evaluation of Coupling Reversed Phase (RP), Aqueous Normal Phase (ANP) and Hydrophilic Interaction (HILIC) Liquid Chromatography with Orbitrap Mass Spectrometry for Metabolomic Studies of Human Urine, Anal. Chem. 84 (2012) 1994–2001.
[11]
I. Gertsman*, and B.A.B., Jon A. Gangoiti, Validation of a dual LC-HRMS platform for clinical metabolic diagnosis in serum, bridging quantitative analysis and untargeted metabolomics, Metabolomics. 10 (2014) 312–323.
[12]
D.S. Wishart, Current Progress in computational metabolomics. Briefings 19
in Bioinformatics. 8 (2007) 279–293. [13]
D.S. Wishart, D. Tzur, C. Knox, R. Eisner, A.C. Guo, N. Young, et al., HMDB: the Human Metabolome Database, Nucleic Acids Res. 35 (2007) D521–6.
[14]
D.S. Wishart, C. Knox, A.C. Guo, R. Eisner, N. Young, B. Gautam, et al., HMDB: a knowledgebase for the human metabolome, Nucleic Acids Res. 37 (2009) D603–D610.
[15]
M. Kanehisa, From genomics to chemical genomics: new developments in KEGG, Nucleic Acids Res. 34 (2006) D354–D357.
[16]
S.G. Smith CA1, O‟Maille G, Want EJ, Qin C, Trauger SA, Brandon TR, Custodio DE, Abagyan R, METLIN: a metabolite mass spectral database, Ther Drug Monit. 27 (2005) 747–751.
[17]
D.L. Wheeler, T. Barrett, D.A. Benson, S.H. Bryant, K. Canese, V. Chetvernin, et al., Database resources of the National Center for Biotechnology Information, Nucleic Acids Res. 35 (2007) D5–D12.
[18]
R. Taguchi, M. Nishijima, T. Shimizu, Basic Analytical Systems for Lipidomics by Mass Spectrometry in Japan, 432 (2007) 185–211.
[19]
J. Zierer, C. Menni, G. Kastenmu, T.D. Spector, Integration of „omics‟ data in aging research: from biomarkers to systems biology, Aging Cell, 14 (2015) 933-944.
[20]
How Tobacco Smoke Causes Disease: The Biology and Behavioral Basis for SmokingAttributable Disease: A Report of the Surgeon General. Atlanta (GA): Centers for Disease Control and Prevention (US); 2010. 9, A Vision for the Future. Available from: http://www.ncbi.nlm.nih.gov/books/NBK53009/
[21]
K. a. Veselkov, L.K. Vingara, P. Masson, S.L. Robinette, E. Want, J. V. Li, et al., Optimized preprocessing of ultra-performance liquid chromatography/mass spectrometry urinary metabolic profiles for improved information recovery, Anal. Chem. 83 (2011) 5864–5872.
[22]
E.J. Want, I.D. Wilson, H. Gika, G. Theodoridis, R.S. Plumb, J. Shockcor, et al., Global metabolic profiling procedures for urine using UPLC–MS, Nature Protocols. 5 (2010) 1005 – 1018.
[23]
J. Xia, I. V. Sinelnikov, B. Han, D.S. Wishart, MetaboAnalyst 3.0-making metabolomics more meaningful, Nucleic Acids Res. (2015) 1–7.
20
[24]
T. Zhang, D.G. Watson, A short review of applications of liquid chromatography mass spectrometry based metabolomics techniques to the analysis of human urine, Analyst. 140 (2015) 2907–15.
[25]
W.M.B. Edmands, P. Ferrari, A. Scalbert, Normalization to Specific Gravity Prior to Analysis Improves Information Recovery from High Resolution Mass Spectrometry Metabolomic Profiles of Human Urine, Anal. Chem. 86 (2014) 10925–10931.
21
Figure Captions
Figure 1: Metabolite profiling workflow to analyze urine samples from nonsmokers and smokers. Samples were prepared for two chromatographic separations and also for two ionization modes. Data analysis was carried out with SIEVE 2.2.
22
23
Figure 2: UHPLC-HRMS chromatograms of RPC and HILIC modes. (A) RPC positive chromatograms of nonsmoker, smoker and QC sample and corresponding extracted ion chromatograms of reserpine. (B) RPC negative chromatograms of nonsmoker, smoker and QC sample and corresponding extracted ion chromatograms of taurocholate. (C) HILIC positive chromatograms of nonsmoker, smoker and QC sample and corresponding extracted ion chromatograms of reserpine. (D) HILIC negative chromatograms of nonsmoker, smoker and QC sample and corresponding extracted ion chromatograms of taurocholate.
24
Figure 3: Normalization based on reserpine and taurocholate standards. The integrated peaks of standards corresponding to (A) RPC positive, (B) RPC negative (C) HILIC positive (D) HILIC negative. Normalization of samples (blanks, nonsmokers, smokers and QCs) before and after showed in the bar diagram.
25
Figure 4: Heatmaps of the urine metabolome profile data and PLS-DA analysis for samples and QCs. Heatmaps of (A) RPC positive, (B) RPC negative, (C) HILIC positive and (D) HILIC negative. Each heatmap represents smoker and nonsmoker groups (columns) and hierarchically clustered metabolites (rows). Color key indicates metabolite expression value, red: highest, green: lowest. The corresponding PLS-DA score plots show clustering of samples and QCs. Shaded regions reflect the 95% confidence intervals. Arrows indicate the outliers in smoker and nonsmoker samples.
26
Figure 5: VIP score plot for top fifteen metabolites from all four modes. VIP score plot for (A) RPC positive, (B) RPC negative, (C) HILIC positive and (D) HILIC negative. Color key indicates metabolite expression value, red: highest, green: lowest. Top 15 metabolites in all four modes are higher in smokers compared to nonsmokers.
27
28
Figure 6: Correlation Chart for Metabolite 1408 and ROC Analysis for HILIC positive. (A) Top 25 metabolites highly correlating (Pearson‟s correlation) with metabolite 1408. Highlights in red indicate metabolites also present in VIP Top 15 list. These were used to construct the (B) ROC curve indicating a highly predictive AUC value of 1.
29