Predicting Network Activity from High Throughput Metabolomics
Knowledge of many metabolic pathways has accumulated over the past century ., For instance , glycolysis , citric acid cycle and oxidative phosphorylation fuel cellular processes through the generation of adenosine triphosphate; glycans and cholesterols not only serve as structural blocks but also mediate intercellular communication ., In fact , metabolites pervade every aspect of life 1 , 2 ., Their roles are increasingly appreciated , as advancing tools allow deeper scientific investigations ., The most notable progresses in recent years come from metabolomics and genome-scale metabolic models ., Metabolomics is the emerging field of comprehensive profiling of metabolites ., As metabolites are the direct readout of functional activity , metabolomics fills in a critical gap in the realm of systems biology , complementing genomics and proteomics 3–6 ., The technical platforms of metabolomics are mainly based on mass spectromety and nuclear magnetic resonance 4 , 7 ., Among them , untargeted LC/MS ( liquid chromatography coupled mass spectrometry ) , especially on high resolution spectrometers , produces unparalleled throughput , measuring thousands of metabolite features simultaneously 5 , 8–10 ., On the other hand , genome-scale metabolic models have been largely driven by genomics , as the total list of metabolic enzymes of a species can be derived from its genome sequence 11 , 12 ., The reconstruction of microbial metabolic network models is an established process 13 , 14 ., Intense manual curation , however , was required in the building of two high-quality human models 15 , 16 , which were followed by a number of derivatives 17–20 ., The coverage of these metabolic models greatly exceeds the conventional pathways ., Even though they are a perfect match in theory , metabolomics and genome-scale metabolic models have had little overlap so far ., The use of metabolomics in building metabolic models is rare 21 , due to the scarcity of well annotated metabolomics data ., The application of genome-scale metabolic models to metabolomics data is not common either 22 ., The limited throughput of targeted metabolomics usually does not motivate large scale network analysis ., Untargeted metabolomics cannot move onto pathway and network analysis without knowing the identity of metabolites ., A typical work flow of untargeted metabolomics is illustrated in Figure 1A ., After ionized molecules are scanned in the spectrometer , the spectral peaks are extracted , quantified and aligned into a feature table ., At this point , each feature is identified by a mass-to-charge ratio ( m/z ) and retention time in chromatography , but its chemical identity is not known ., To assign a spectral feature to a bona fide metabolite , it usually involves tandem mass spectrometry to examine the fragmentation pattern of a specific feature , or coelution of isotopically labeled known references - both are inherently low throughput ., Considerable effort is needed to build a spectral library , which is often of limited size and interoperability ., Thus , metabolite identification forms the bottleneck of untargeted metabolomics 23 ., A number of informatics tools have been developed for LC/MS metabolomics , ranging from feature extraction 24–26 , pathway analysis and visualization 27–29 to work flow automation 30–32 ., Yet , whereas pathway and network analysis is concerned , the existing tools require identified metabolites to start with ., Computational prediction of metabolite identity , based on m/z alone , is deemed inadequate as a single m/z feature can match multiple metabolites even with high instrumental accuracy 33 , 34 , and multiple forms of the same metabolite often exist in the mass spectra 35 ., Although automated MS/MS ( tandem mass spectrometry ) search in databases is improving the efficiency of metabolite identification 36 , 37 , this requires additional targeted experiments and relies on extensive databases , where data from different platforms often do not match ., How to bring untargeted metabolomics data to biological interpretation remains a great challenge ., In this paper , we report a novel approach of predicting network activity from untargeted metabolomics without upfront identification of metabolites , thus greatly accelerating the work flow ., This is possible because the collective power in metabolic networks helps resolve the ambiguity in metabolite prediction ., We will describe the computational algorithms , and demonstrate their application to the activation of innate immune cells ., The genome-scale human metabolic network in mummichog is based on KEGG 38 , UCSD Recon1 15 and Edinburgh human metabolic network 16 ., The integration process was described previously 39 ., The organization of metabolic networks has been described as hierarchical and modular 40 ., When we perform a hierarchical clustering on the metabolic reactions in our network , its modular structure is clear ( Figure 2A ) ., This modular organization , as reported previously 41 , often but not always correlates with conventional pathways ( Figure 2B ) ., The module definition in this work is adopted from Newman and Girvan 42 , 43 , where a module is a subnetwork that shows more internal connections than expected randomly in the whole network ., Modules are less biased than pathways , which are defined by human knowledge and conventions , and outgrown by genome-scale metabolic networks ., Activity of modules may exist within and in between pathways ., Deo et al 22 convincingly demonstrated the advantage of unbiased module analysis over pathways ., On the other hand , pathways have built-in human knowledge , which may be more sensitive under certain scenarios ., Pathway analysis and module analysis are rather complementary , and both are included in mummichog ., The reference metabolic network model contains both metabolites and enzymes ., Since metabolomics only measures metabolites , the model is converted to a metabolite centric network for statistical analysis ., Enzymes are only added later in the visualization step to aid user interpretation ., Within the predefined reference metabolic network model , mummichog searches for all modules that can be built on user input data , and compute their activity scores ., This process is repeated many times for the permutation data to estimate the background null distribution ., Finally , the statistical significance of modules based on user data is calculated on the null distribution ., The specific steps are as follows: The basic test for pathway enrichment here is Fishers exact test ( FET ) , which is widely used in transcriptomic analysis ., The concept of FET is , when we select features ( ) from a total of features ( ) , and find of the features present on a pathway of size , the chance of getting in theory can be estimated by enumerating the combinations of , and ., To apply FET to an enrichment test of metabolic features on pathways , we need to understand the additional layer of complexity ., Our metabolic features can be enumerated either in the m/z feature space or in the metabolite ( true compound ) space ., Since metabolic pathways are defined in the metabolite space , either way needs to factor in the many-to-many mapping between m/z features and metabolites ( Figure S1 ) ., This mapping is effectively covered in our permutation procedure , which starts from the m/z features and reruns the mapping every time ., The overall significance of a pathway enrichment is estimated based on a method by Berriz et al 44 , which ranks the p-value from real data among the p-values from permutation data to adjust for type I error ., Finally , a more conservative version of FET , EASE , is adopted to increase the robustness 45 ., The key idea of EASE is to take out one hit from each pathway , thus preferentially penalize pathways with fewer hits ., The specific steps are as follows: Both the module analysis and pathway analysis above serve as a framework to estimate the significance of functional activities ., In return , the predicted metabolites in significant activities are more likely to be real ., Mummichog collects these metabolites , and look up all their isotopic derivatives and adducts in ., A confidence rating system is applied to filter for qualified metabolites ., For instance , if both the single-charged form M+H1+ and the form M ( C13 ) +H1+ are present , this metabolite prediction carries a high confidence ., All the qualified metabolites carry over their connections in the reference metabolic network , and form the “activity network” for this specific experiment ( e . g . Figure 3 ) ., The activity network gears towards a quality and clear view of user data , as modules and pathways can be redundant and fragmented ., It also accentuates the activity in a specific experimental context , in contrast to the generic nature of the reference metabolic network ., We next illustrate the application of these algorithms to a novel set of immune cell activation data , and two published data sets on human urinary samples and yeast mutants ., The innate immunity plays a critical role in regulating the adaptive immunity , and the field was recognized by the 2011 Nobel Prize in Physiology or Medicine 46 ., According to the nature of stimuli , innate immune cells direct different downstream molecular programs , which are still under intense scientific investigation 47 , 48 ., In this study , we examine the metabolome of human monocyte-derived dendritic cells ( moDC ) under the stimulation of yellow fever virus ( YF17D , a vaccine strain ) ., We have shown previously that yellow fever virus activates multiple toll-like receptors , and induces cellular stress responses 49–51 ., This data set is , to our knowledge , the first high throughput metabolomics on any immune cells ( macrophages were previously studied by limited throughput ) ., The cell extracts from our activation experiment were analyzed by LC/MS metabolomics , and yielded 7 , 995 spectral features ( denoted as ) after quality control ., Among them , 601 features were significantly different between the infected samples and both the baseline and time-matched mock controls ( , student t-test; denoted as ) ., Using and , mummichog computes significant pathways and modules and the activity network ., Viral infection induced a massive shift of metabolic programs in moDCs ( pathways in Table S1 , modules in Figure S2 ) ., The predicted activity network is shown in Figure 3A , and we will focus our investigation on a small subnetwork ( Figure 3B ) , which includes the metabolisms of nucleotides , glutathione/glutathione disulfide and arginine/citrulline ., Nucleotides are required for viral replication , and the hijacking of host nucleotide metabolism by virus has been well described 52–54 ., Glutathione is best known as intracellular antioxidant , where it is oxidized to glutathione disulfide ( GSSG ) ., However , our data show that both glutathione and GSSG are depleted in activated moDCs , departing from this conventional wisdom ., The across-the-board depletion is consistent with the down-regulation of genes for glutathione synthesis ( Figure 4B ) ., Our data support the notion that glutathione is released by dendritic cells and conditions the extracellular microenvironment during their interaction with T cells 55–57 ., Arginine is known to be an important regulator in innate immune response 58 , 59 ., Arginine metabolism can lead to two pathways: to ornithine ( catalyzed by arginase ) or to citrulline ( catalyzed by nitric oxide synthase ) ., The decrease of arginine and increase of citrulline suggests the latter pathway , which is the main reaction of producing intracellular nitric oxide ., We indeed detected the inhibition of eNOS and iNOS expression later ( Figure 4C ) , a well documented feedback effect by nitric oxide 60 ., We also performed tandem mass spectrometry on the metabolites in Figure 3B , using authentic chemicals as references ., All the metabolites , except glutamylcysteine and thyroxine , were confirmed ( Figure 5 , Figure S3 ) ., The depletion of arginine and accumulation of citrulline in moDC was also triggered by dengue virus but not by lipopolysaccharide ( LPS , Figure S4 ) ., It is known that iNOS is induced in dendritic cells by LPS but not by virus 47 , 61 ., Our data suggest a different nitric oxide pathway in viral infection , driven by constitutive nitric oxide synthases ., The intracellular nitric oxide has a fast turnover and we did not detect its accumulation by fluoremetric assays ( data not shown ) ., We previously demonstrated that the phosphorylation of EIF2A was induced by YF17D 50 ., An upstream mechanism is now suggested by this metabolomic experiment , as both the production of nitric oxide and depletion of arginine induce the activity of EIF2A kinases 62 ., The nature of metabolomics data often varies by platforms and sample types ., We thus extend our mummichog approach to two published data sets on human urinary samples 63 and on yeast cell extracts 64 ., Both data sets carry metabolite annotation by the original authors , which can be used to evaluate the prediction by mummichog ., The human urinary data contained both formal identification by matching to local library of chemical references and putative identification by combining multiple public resources 63 ., We used mummichog to investigate the gender difference in this data set , and predicted an activity network of 45 metabolites ., Among them , 13 were not found in the original annotation ., For the remaining metabolites , 97% ( 31/32 ) were agreed between mummichog and the original annotation ( Figure 6 ) ., There is an option in mummichog to enforce the presence of M+H+ form ( for positive mode , M−H− for negative mode ) in metabolite prediction ., With this option , 3 out of 44 metabolites were not in the original annotation , and the remaining 41 metabolites were in 100% agreement ., The mummichog algorithms are not tied to a specific metabolic model ., We adopted the yeast metabolic model from BioCyc database 11 for the yeast data 64 , to predict an activity network contrasting mutant and wild type strains ., This data set was only annotated for 101 metabolites through the authors local library ., As a result , the majority of metabolites in the predicted network by mummichog were not found in the original annotation ., Out of the remaining 28 metabolites , 24 ( 86% ) were agreed between mummichog and the original annotation ( Figure 6 ) ., Enforcing the presence of primary ion M−H− ( data collected in negative ion mode ) had little effect to the result , since the original annotation was already biased to metabolites that are ionized easily ., These results show that the prediction by mummichog is robust cross platforms and sample types ., Critical to the success of mummichog is the integration of genome-scale metabolic models ., We have used in this study a recent human metabolic model ., An alternative human model from BioCyc 11 produced comparable results ( Figure S6 ) ., The coverage of the models in all three case studies is shown in Table 1 ., These genome-scale metabolic models are more extensive than conventional pathways , and were shown to capture activities in between pathways 22 ., The pathway organizations differ between the two human models , as the BioCyc model tends to use smaller pathways ., This creates some model dependency in the pathway analysis , but little effect to the “activity network” , as mummichog is more network centric ., The two test cases in Figure 6 also indicate that these models tend to capture more information than conventional annotations ., However , as mentioned earlier , the new data from metabolomic studies are yet to be integrated into these genome-scale metabolic models ., For example , a number of metabolites in metabolomics databases 36 , 65 , 66 are not in any of these metabolic models ., In general , the features from a high resolution profiling experiment by far exceed the current annotations in metabolite databases ., This leads to a de facto filtering when data are run on mummichog ( similiar situation in database searches ) ., Meanwhile , the features that can be mapped to the current metabolic model are more likely to be biologically relevant ., This “filtering” is pertinent to the metabolic model , not to mummichog algorithms - mummichog still has to choose the true metabolites from multiple possible candidates ( Figure S1B ) ., It will be an important future direction to advance metabolic modeling with the chemical data ., We also expect the metabolic models to improve on lipid annotation , physiological context and tissue specificity ., As lessons learned from transcriptomics , pathway and network analysis not only provides functional context , but also the robustness to counteract noises at individual feature level , which are commonly seen in omics experiments ., Similarly , the prediction on activity by mummichog is tolerant to errors at individual feature level ., In the moDC data , we chose by a cutoff value ., When we vary this cutoff from to , the program returned networks of a stable set of metabolites ( Figure S7 ) ., The module finding procedure in the program was designed to extensively sample subnetwork structures ., Among the modules will be many variations , but the subsequenct “activity network” will collapse on stable results ., In deed , we tested an alternative algorithm of modularization 67 , and it returned almost identical predicted networks , in spite of moderately different intermediate modules ( Figure S8 ) ., In theory , there are merits to incorporate a statistical matrix from the feature selection step prior to mummichogs analysis and mass flow balance of metabolic reactions 22 , 68 ., While these are appealing directions for future research , the current version of mummichog confers some practical robustness , such as tolerance to technological noise and biological sampling limitation ., For example , mass balance is almost impossible within serum or urine samples , because the reactions producing these metabolites are likely to occur in other tissues ., The number of overlap metabolites is used in the enrichment calculation in both module analysis and pathway analysis ., Sometimes , a single m/z feature may match to several metabolites in the same module/pathway , inflating the overlap number ., Thus , mummichog always compares the number of overlap metabolites and the number of corresponding m/z features , and uses the smaller number for enrichment calculation , since the smaller number is more likely to be true ., The size of each metabolic pathway is defined by the number of metabolites in the pathway ., mummichog uses only the metabolites that can be matched in to define a pathway size , because this reflects the analytical coverage of the experiment and is confined by the same coverage ., Overall , mummichog uses the whole feature list from an experiment for resampling , therefore the computation of statistical significances effectively circumvents analytical biases ., In spite of the fantastic progress in mass spectrometry , these are the early days of metabolomics ., Effective computational integration of resources , the combination of cheminformatics and bioinformatics , will greatly benefit the field 69 , 70 ., As data accumulate , further method refinement will become possible ., Mummichog presents a practical solution of one-step functional analysis , bypassing the bottleneck of upfront metabolite identification ., It trades off some sensitivity in the conventional approach for the much accelerated work flow of high throughput LC/MS metabolomics ., Mummichog is not designed to replace tandem mass spectrometry in metabolite identification ., It is the biological activity not metabolites per se that mummichog predicts ., Even with some errors on individual metabolites , as long as the biology is pinpointed to a subnetwork structure , investigators can focus on a handful of validations , steering away from the lengthy conventional work flow ., In conclusion , we have demonstrated that mummichog can successfully predict functional activity directly from a spectral feature table ., This benefits from the convergence of genome-scale metabolic models and metabolomics ., Mummichog will continue to improve as the metabolic network models evolve ., We expect this method to greatly accelerate the application of high throughput metabolomics ., The mummichog software is available at http://atcg . googlecode . com ., Human peripheral blood mononuclear cells ( PBMCs ) were isolated from Buffy coats by separation over a Lymphoprep gradient ., CD14+ monocytes were isolated from the PBMCs with MACS beads ( Miltenyi Biotec , Auburn , CA ) and cultured for 7 days with 20 ng/ml GM-CSF and 40 ng/ml IL-4 ( Peprotech , Rocky Hill , NJ ) ., MoDCs were then harvested , washed twice and resuspended in serum-free medium ., MoDCs ( ) were stimulated in triplicate in 48-well plates in a 200 µL volume with Yellow Fever virus ( M . O . I . of 1 ) , or mock infected ., After 2 hrs , 800 µL of 10% FBS-RPMI was added to all wells ., MoDCs were harvested at 6 hr or 24 hr after infection and centrifuged ., Supernatants were aspirated , and dry cell pellets were frozen at −80°C ., Supernatants of moDC cultures were assayed for the concentration of IL-6 and TNF using ELISA kits ( BD , San Diego , CA ) ., Three biological replicates were used for LC/MS and QPCR ., Full scan LC/MS ( m/z range 85–2000 ) was performed essentially as previously described 8 ., Cell extracts or supernatants were treated with acetonitrile ( 2∶1 , v/v ) and centrifuged at 14 , 000× g for 5 min at 4°C to remove proteins ., Samples were maintained at 4°C in an autosampler until injection ., A Thermo Orbitrap-Velos mass spectrometer ( Thermo Fisher , San Diego , CA ) coupled with anion exchange chromatography was used for data collection , via positive-ion electrospray ionization ( ESI ) ., Metabolites of interest were identified by tandem mass spectrometry on a LTQ-FTMS , where the biological sample , biological sample spiked in with authentic chemical and authentic chemical reference were run sequentially ., The and were done in the ion trap of the LTQ-FTMS , with an isolation width of 1 amu and a normalized collision energy of 35 eV ., The LC/MS data were processed with apLCMS program 25 for feature extraction and quantification ., Significant features were also verified by inspecting the raw data ( Figure S5 ) ., Features were removed if their intensity is below 10 , 000 in every sample class ., Missing intensity values were imputed to 1000 ., The intensities were log2 transformed ., Low quality features were further filtered out if their averaged in-class coefficient of variation was greater than 0 . 2 ., Averaged ion intensity over three machine replicates was used for subsequent analysis ., These 7 , 995 features constituted the reference list ., No normalization was used because total ion counts in all samples were very similar ., Students t-test was used to compare infected samples ( at 6 hr ) versus mock infections ( at 6 hr ) , and infected samples ( at 6 hr ) versus baseline controls ( at 0 hr ) ., Features with in both tests were included in the significant list ., The feature table , , and predictions are given in Dataset S1 ., For gene expression quantification , mRNA was extracted by RNeasy Mini Kit ( Qiagen ) according to manufacturers protocol , where the cell lysate was homogenized by QIAshredder spin columns ., Reverse transcription was performed with SuperScript III reverse transcriptase and oligo-dT ( Invitrogen ) according to manufacturers recommendation ., Real-time PCR was performed on a MyiQ Icycler ( BioRad ) , using SYBR Green SuperMix ( Quanta Biosciences ) ., The PCR protocol used 95°C 3 mins; 40 cycles of 95°C 30 seconds and 60°C for 60 seconds ., The amplicons were verified by melting curves ., Quantafication was performed by the method , normalized by microglobulin levels ., The primer sequences are given in Table S2 ., Data on human urinary samples 63 were retrieved from MetaboLights server 71 ., The positive ion feature table for study “439020” contained 14 , 720 features ., A feature is only included if its ion intensity is above 100 , 000 in 5 or more samples ., This leaves 11 , 086 features , which consist for this study ., Data were normalized by total ion counts ., We next compared the metabolite difference between females ( 8 samples of low testosterone glucuronide level ) and males ( 11 samples of high testosterone glucuronide level ) ., is consisted of 524 features ( by student t-test ) ., The original authors annotated 3 , 689 metabolite features , and their annotation was used to compare with the prediction by mummichog ., The yeast data 64 were downloaded from MAVEN website 32 in mzXML format ., Feature extraction was performed in MAVEN through two approaches: targeted approach and untargeted approach ., The targeted approach used chemical library from the same lab and produced 177 annotated features , which corresponded to 101 metabolites ., The untargeted approach extracted 6318 features without annotation ., After the same processing procedure as in our moDC data , contained 5707 features ., We thus used mummichog to predict on the untargeted data , and compared the result to MAVEN annotation ., The consisted of 426 features that were significantly different between the prototrophic wild type and the auxotrophic mutant ( by student t-test ) ., The yeast metabolic model was compiled from BioCyc data 11 .
The functional interpretation of high throughput metabolomics by mass spectrometry is hindered by the identification of metabolites , a tedious and challenging task ., We present a set of computational algorithms which , by leveraging the collective power of metabolic pathways and networks , predict functional activity directly from spectral feature tables without a priori identification of metabolites ., The algorithms were experimentally validated on the activation of innate immune cells .
Mass spectrometry based untargeted metabolomics can now profile several thousand of metabolites simultaneously ., However , these metabolites have to be identified before any biological meaning can be drawn from the data ., Metabolite identification is a challenging and low throughput process , therefore becomes the bottleneck of the filed ., We report here a novel approach to predict biological activity directly from mass spectrometry data without a priori identification of metabolites ., By unifying network analysis and metabolite prediction under the same computational framework , the organization of metabolic networks and pathways helps resolve the ambiguity in metabolite prediction to a large extent ., We validated our algorithms on a set of activation experiment of innate immune cells ., The predicted activities were confirmed by both gene expression and metabolite identification ., This method shall greatly accelerate the application of high throughput metabolomics , as the tedious task of identifying hundreds of metabolites upfront can be shifted to a handful of validation experiments after our computational prediction .
systems biology, metabolic networks, biology, computational biology