Oncotarget

Research Papers:

The integrative metabolomic-transcriptomic landscape of glioblastome multiforme

PDF |  HTML  |  How to cite  |  Order a Reprint

Oncotarget. 2017; 8:49178-49190. https://doi.org/10.18632/oncotarget.16544

Metrics: PDF 2032 views  |   HTML 2709 views  |   ?  

Dieter Henrik Heiland _, Jakob Wörner, Jan Gerrit Haaker, Daniel Delev, Nils Pompe, Bianca Mercas, Pamela Franco, Annette Gäbelein, Sabrina Heynckes, Dietmar Pfeifer, Stefan Weber, Irina Mader and Oliver Schnell

Abstract

Dieter Henrik Heiland1,5, Jakob Wörner2, Jan Gerrit Haaker1,5, Daniel Delev1,5, Nils Pompe2, Bianca Mercas1,5, Pamela Franco1,5, Annette Gäbelein1,5, Sabrina Heynckes1,5, Dietmar Pfeifer3,5, Stefan Weber2, Irina Mader4,5 and Oliver Schnell1,5

1Department of Neurosurgery, Medical Center, University of Freiburg, Freiburg, Germany

2Institute of Physical Chemistry, Faculty of Chemistry and Pharmacy, University of Freiburg, Freiburg, Germany

3Department of Hematology, Oncology and Stem Cell Transplantation, Medical Center, University of Freiburg, Freiburg, Germany

4Department of Neuroradiology, Medical Center, University of Freiburg, Freiburg, Germany

5Faculty of Medicine, University of Freiburg, Freiburg, Germany

Correspondence to:

Dieter Henrik Heiland, email: dieter.henrik.heiland@uniklinik-freiburg.de

Keywords: metabolomics, transcriptomics, network analysis, glioblastoma multiforme, WGCNA

Received: January 02, 2017    Accepted: February 23, 2017    Published: March 24, 2017

ABSTRACT

The purpose of this study was to map the landscape of metabolic-transcriptional alterations in glioblastoma multiforme. Omic-datasets were acquired by metabolic profiling (1D-NMR spectroscopy n=33 Patient) and transcriptomic profiling (n=48 Patients). Both datasets were analyzed by integrative network modeling. The computed model concluded in four different metabolic-transcriptomic signatures containing: oligodendrocytic differentiation, cell-cycle functions, immune response and hypoxia. These clusters were found being distinguished by individual metabolism and distinct transcriptional programs. The study highlighted the association between metabolism and hallmarks of oncogenic signaling such as cell-cycle alterations, immune escape mechanism and other cancer pathway alterations. In conclusion, this study showed the strong influence of metabolic alterations in the wide scope of oncogenic transcriptional alterations.


INTRODUCTION

Glioblastoma multiforme (GBM) is the most common primary malignant brain tumor in adults, with an annual incidence of 3–4 cases per 100 000 people in Europe [1, 2] and the United states [3]. In spite of the best available treatment, the prognosis for patients with GBM is poor, with a median survival of not more than 14–16 months [4-8]. During the last decade, numbers of novelties in cancer development, genetic and metabolic alterations were discovered in the context of glioma [9]. Most notably, altered metabolism was described as a hallmark of cancer development and frequently occurs among different cancer types [10]. Furthermore, tumorigenesis is driven by the reprogramming of cellular metabolism, which derives directly or indirectly from genetic or epigenetic alterations [10]. The recently published revised WHO classification of brain tumors accounts for genetic alterations highlighting in gliomas especially the IDH mutation [11]. In glioma, the “Glioma-CpG island methylator phenotype“(G-CIMP) has been found to arise as a consequence of tumor-associated metabolic alterations [12, 13]. Usually, isocitric acid is transformed into α-ketoglutaric acid by IDH1/2. In case of a catalytic located IDH1/2 mutation, isocitric acid is converted into the oncometabolite 2-hydroxyglutaric acid (2-HG) [12, 13]. Accumulation of 2-HG reshapes the tumor methylome and constitutes the G-CIMP [12, 13]. Interestingly, patients affected by the IDH1/2 mutation reveal a better clinical course compared to non-mutated patients [14]. In particular, the IDH1/2 mutation and its metabolic changes pinpointed the strong coherence between metabolic and genetic alterations and highlighted the relevance of tumor-metabolism in glioblastoma multiforme.

Metabolic/genetic profiling in glioma

The first observations on metabolic alterations in gliomas were made in 1H-magnetic resonance spectroscopy (H-MRS) some decades ago [15, 16]. Aided by array-based molecular methods and improved analytic approaches, new studies in cancer cell metabolism were expanded to the understanding of the mechanisms and functional consequences of tumor-associated metabolic alterations. Nuclear magnetic resonance (NMR) and 1H-MRS are established tools for metabolic profiling in-vivo or ex-vivo. A recent study analyzed nine glioma-cell lines by NMR spectroscopy [17]. The metabolic profiles were compared to publicity available gene expression data of these cell-lines. Four distinct profiles were detected, which were classified according to specific metabolic alterations [17]. In another study, the metabolic/proteomic signature of glioma-cell-lines was investigated. Eight metabolites were detected by NMR and compared to RPPA (reverse-phase protein array) proteomic analysis and expression data. The cell-lines were then clustered into “full-stem” and “restricted stem” subgroups, based on their stem-cell capacity [18]. Moreover, in a study performed by Chinnaiyan et al., 2012, where metabolic analysis was performed by liquid chromatography, altered anaerobic metabolism in the mesenchymal subgroup described by Verhaak and Phillips [19, 20] was shown. Pantel et al., 2014 analyzed tumor heterogeneity by single-cell RNA sequencing and investigated 4 subgroups of transcriptomic profiles, which co-existed within the same tumor [21]. The profiles were named by dominant biological function as oligodendrocytic differentiation, cell-cycle, immunoresponse and hypoxia [21]. Therefore, environmental conditions as hypoxia with its metabolic alteration may influence tumor heterogeneity and associated expression profiles in GBM and other cancers [22-24].

The purpose of this study was to integrate metabolomic and transcriptomic data by comprehensive network-based modeling. This analytic approach attempts to improve the knowledge about interacting regulation mechanism between metabolic and transcriptional alteration in glioblastoma multiforme.

RESULTS

Workflow and overview of combined methods

This study contained a pipeline including several steps from, tissue sampling guided by neuronavigation, genetic and metabolic engineering up to bioinformatic analysis. In short, tissue samples were subjected to a combined genetic and spectroscopic analysis. Metabolites from 1H NRM spectroscopy and genetic expression data were normalized. Classical analysis (Consensus Cluster (k=5), unsupervised clustering and survival Analysis) of transcriptomic and metabolic data was performed separately on both datasets. For data integration, a network analysis using a topological overlapping measurement identified genetic modules [25]. Associations of each gene and the traits of interest (metabolites) were quantified by defining Gene Significance (GS) as the correlation between each gene and metabolite. Additionally, module membership (MM) for each gene was measured by correlation of module eigengene and the gene expression profile. The intramodule connectivity (kME) was calculated using the GS and MM measures. Therefore, genes were defined by its correlation with specific metabolites as well as its importance (ranked within its module membership). The kME allowed an identification of those genetic modules, which were associated to specific intracellular metabolism. All steps of the workflow were illustrated in Figure 1.

The figure reveals the workflow and data processing of the “in-house” pipeline.

Figure 1: The figure reveals the workflow and data processing of the “in-house” pipeline. This semi-automated analysis served as a robust method for integrative analysis of metabolic and genetic data.

Separate metabolic and transcriptomic profiling of glioblastoma multiforme

56 patients with neuronavigation-guided intraoperative sampling, available fresh frozen tissue samples and a histological confirmed glioblastoma multiforme were included in this study. RNA preparation and transcriptomic profiling was performed in 48 cases (low RNA quality in 8 cases), a total number of 38 tissue samples were analyzed by NMR (low tissue quality in 8 cases, low metabolite quality in 10 cases) and 4 samples not achieved post-processing quality criteria (fitting not successful).

Raw spectra (Figure 2, upper panel) were processed by a Bayesian-algorithm implemented in the R-software packed “Batman” (Figure 2, lower panel, detailed description in the method part). An unsupervised hierarchical clustering of the normalized metabolic intensity values of all tissue samples identified three clusters with distinct metabolic profiles (Figure 3C). Cluster yellow was associated to the proneural subgroup and revealed a significant better clinical course in comparison to the other cluster groups. Cluster green showed an association to the mesenchymal expression group and showed a poor clinical course. In comparison to the separate metabolic analysis the transcriptomic profiles of 48 patients were explored. A consensus clustering showed 3 main subgroups, which were confirmed by unsupervised hierarchical clustering (Figure 3A3B). The expression subgroups of each patient were predicted by a random forest algorithm, which based on a training dataset (TCGA data set of 484 patients) with implemented subgroups by Verhaak (Figure 3B bottom panel). Three subgroups were found with distinct transcriptomic profiles. The proneural (PN) and neural (N) group was not separately grouped weather classical (CL) and mesenchymal (MES) samples grouped into distinct cluster. A survival statistics revealed no strong differences between all subgroups, only the PN subgroup revealed a significant improved survival compared to the MES subgroup. This analysis showed a clear separation of transcriptome and metabolic profiles into the described subgroups of Verhaak. To quantify these correlations an integrative analysis of both omic-datasets was performed.

High resolutionNMR-Spectra of one patient including a raw spectra (upper panel) and fitted curves (lower panel).

Figure 2: High resolutionNMR-Spectra of one patient including a raw spectra (upper panel) and fitted curves (lower panel).

Figure 3:

Figure 3: (A) Consensus cluster of expression data (n=48) revealed 4 cluster groups, which were summarized in an unsupervised clustering (B). The bar below the heatmap indicated the expression subgroups (identified by random forest analysis). (C) Unsupervised hierarchical cluster of normalized metabolite values. Bars below the heatmap describe the expression subgroup of each patient. (D-E) Survival analysis of all clustergroups (derived from transcriptome (D) and metaboliom (E)) shows a significantly different overall survival with a more favorable outcome for the proneural subgroup of proneural-associated metabolic cluster.

Combined integrative network-based analysis of metabolic and transcriptome profiling

For further characterizing of each metabolite, a gene expression network with integrated metabolic data was build. Genes were summarized in 79 expression modules. Module-containing genes were ranked by their intramodule connectivity and correlated to each metabolite. These correlation coefficients were analyzed by an unsupervised hierarchical clustering (Figure 4A). Highly correlated metabolites and expression modules revealed four specific subgroups colored in green (Cluster I), yellow (Cluster II), violet (Cluster III) and red (Cluster IV).

Figure 4:

Figure 4: (A) Unsupervised hierarchical cluster of correlation coefficients (kME and normalized metabolite values). High correlations was colored in red, low correlation in blue. (B) Gene set enrichment analysis plots summarized enrichment scores of indicated biological functions. Enrichment Scores (ES) and p-values of all expressions modules (derived from WGCNA) were illustrated in a volcano plot. On the y-axis, the negative logarithm of GSEA p-values were presented, the x-axis contained ES-values. The size of each point indicated the level of gene set enrichment.

Cluster I (green): proneural/oligo-related cluster

Cluster I contained 10 metabolites and 16 expression modules. These expression modules were analyzed by gene set enrichment analysis (GSEA, detailed description in the method part). The results of the GSEA were summarized in a volcanoplot (Figure 3B). This plot illustrated each enrichment-score and the corresponded negative log of FDR p-values of each module within cluster 1. Most expression modules that were presented in cluster I significantly enriched proneural genes. In addition, oligodendrocytic genes were also highly enriched in most of the expression modules summarized in cluster I (Figure 4B). Creatine, was identified as “key-metabolite” defined as highest total kME correlation in each cluster group, which was significantly associated with all expression modules contained in cluster I.

Cluster II (yellow): immune-related cluster

Cluster II comprised 5 metabolites and 8 expression modules. A GSEA was used for further functional characterization. An analysis of the known glioblastoma multiforme expression subgroups showed a stronger overall enrichment for the mesenchymal gene set. Interestingly, immune response and activation of the immune system were highly enriched in cluster (Figure 4B). The key-metabolites of this cluster were phosphocholine and choline being significantly correlated to all related expression modules.

Cluster III (violet): cell-cycle-related cluster

Cluster III included 9 metabolites and 11 expression modules. Four expression modules showed a strong association with the proneural subgroup, while others enriched mesenchymal genes. Further analysis of other biological functions revealed a strong connection with cell-cycle functions including regulation of Phase M2, DNA repair mechanism and regulation of mitosis (Figure 4B). The metabolite glycine was identified as the key-metabolite of this cluster. It was significantly correlated to all related expression modules.

Cluster IV (red): hypoxia/mesenchymal-related cluster

Cluster IV contained 7 metabolites and 6 expression modules. A strong enrichment of mesenchymal genes was found by GSEA. Only one expression module showed a non-significant enrichment of proneural genes. Key-metabolite of this cluster was lactate, which is a known agent of anaerobic metabolism (Figure 4B). Interestingly, expression modules of cluster IV highly enriched more genes that were up-regulated in hypoxic conditions than those being down-regulated.

All connections between metabolite and correlated genes were given in a connectivity-weighted network and illustrated in Figure 5. This scheme represents the relationship between the metabolites and the genetic modules. It is apparent from this that some metabolites, such as myo-inositol and lactate, also have connections to further clusters than only to those identified primarily.

Figure 5:

Figure 5: (A) Integrative network of metabolite and expression modules (WGCNA). Size and color indicated the importance of each gene/metabolite in the network. Hub-genes/metabolites were colored in red.

Integrative interpretation of clusters

Individual expression modules of each identified cluster (1-4) were separately analyzed by gene-set-variation analysis (GSVA) to find highly enriched metabolic-pathways (KEEG) (Figure 6A). Cluster 1 and 3 dominantly expressed proneural genes. These clusters showed an overlap of enriched metabolic pathways (Figure 6A), including arginine, proline-, glutathione-, and amino-sugar-metabolism. Cluster 3, which contained cell-cycle functions, exclusively enriched glycerolipid- and sphingolipid metabolism. Cluster 2 and 4 highly enriched mesenchymal genes and were characterized as immune response- or hypoxia-related. These clusters showed an overlap of metabolic pathways (pathways in line 22 – 25 on Figure 6A), spotlighted by the pyruvate metabolism. Most notably, Cluster 2 contained expression modules, which enriched dominantly genes of the immune response and exclusively enriched the tryptophan metabolism.

Figure 6:

Figure 6: (A) KEEG metabolism-pathway of Cluster 1-4 was illustrated. The exclusive and overlapping enriched pathways were marked. The color code at the bottom showed the dominant expression subgroup and indicated the functional subgroup. (B) Map of metabolic differences on single-metabolite level and expression of enzymes belonging to the mapped metabolic pathways. *** p<0.001, ** p<0.01, *p<0.05.

A map of metabolites associated with the proneural subtype (creatine) including the creatine-, arginine-, glycine-, serine- and pyruvate metabolism as illustrated in Figure 6B. Enzymes of the creatine degradation were significantly stronger expressed in mesenchymal samples, while enzymes of the arginine-proline metabolism were up-regulated in proneural tumors (Figure 6B).

DISCUSSION

Gene expression and cell metabolism are frequently altered in glioblastoma and support cell proliferation, epigenetic alterations and tumor aggressiveness. In 2010, Verhaak proposed four glioblastoma multiforme subclasses by expression profiling [19]. These subgroups were characterized by specific gene sets, which were exclusively expressed in each subclass. Patel et al. (2014) showed a strong variance of the Verhaak expression subgroups based on single-cell sequencing of up to 100 cells within the same tumor [21]. Intratumoral heterogeneity within those subgroups could be demonstrated, and four transcriptional modules were additionally identified on single cell lines within the same tumor (Hypoxia, Immune, Oligo, Cell Cycle). These metabolic profiles reflect the final downstream product of gene transcription [26]. Therefore the purpose of this study was to combine genetic and metabolomic information to identify the transcriptomic-metabolomic landscape of glioblastoma multiforme, which mirrors the up- and downstream products of cellular regulation mechanisms.

In the first cluster analysis containing only spectral intensities, a significant association of metabolism to the Verhaak expression subgroups was found (Figure 3C). This connection was to be expected and has already been described in the literature [27]. In the next step, a network model was developed that integrates functional expression modules and correlated metabolites of GBM. This model resulted in four different cluster-groups, illustrated in Figure 4. Each cluster was interpreted in the context of its dominant underlying biological functions and contained: oligodendrocytic differentiation (Cluster 1), cell-cycle (Cluster 3), immune response (Cluster 2) and hypoxia (Cluster 4). Cluster 2 (immune response) and 4 (hypoxia) highly enriched mesenchymal genes. They showed an overlap of their enrichment of metabolic-pathway signatures, which highlight the pyruvate metabolism. Pyruvat plays an important role by connecting anaerobic and aerobic metabolism [28]. The immune response related cluster additionally enriched exclusively the tryptophane metabolism. A study by Moffett et al., 2003 described the major role of the tryptophane metabolism in several immune functions including the immune escape mechanism [29]. The key metabolite of the immune response-related cluster was choline. Although choline-containing compounds are generally discussed to be involved within membrane [30] turnover, other authors discussed the choline-containing compound to be related to immune response in crap [31].

The key metabolite of the hypoxia-related cluster was L-lactic acid. As generally known, lactate acid is used for energy metabolism in hypoxic environment [10, 27]. Hypoxic conditions were also known for a proneural-to-mesenchymal transition, which revealed strong coherence between metabolic environment and transcriptional regulation mechanism [32]. The metabolic and transcriptional differences of the pyruvate and lactate metabolism are presented in Figure 6. Our results point to the fact that metabolic alterations may drive the appearance of specific glioblastoma multiforme expression subclasses as a specific adaption upon microenvironmental conditions.

Clusters 1 and 3 dominantly enriched proneural genes and therefore belonged to the proneural subtype. Furthermore, cluster 1 showed a strong enrichment of oligodendrocytic differentiation. The key-metabolite was creatine, which is in line with the literature [33], where creatine was found to be highly associated with increased expression of proneural genes in a MRS in-vivo study [33]. Additionally, proneural samples were typically characterized by an up-regulation of oligodendrocytic genes as OLIG1/2 [19]. The cluster 3 contained cell-cycle functions and showed an exclusive enrichment of the glycerolipid and sphingolipid pathway. This metabolic alteration is well known and highly associated with G0/G1- cell-cycle functions and other parts of cell-cycle regulation [34]. Its key-metbolite was glycine, which has several unspecific functions and could metabolized in several substrates (KEGG). A specific connection between glycine and cell cycle function has not been reported so far.

The study by Patel et al., 2014 analyzed single-cell transcriptomic profiles of glioblastoma multiforme within the same tumor. All profiles were summarized in four clusters containing: oligodendrocytic differentiation, cell-cycle, immune response and hypoxia. These results were totally in-line with our reported findings. In fact, tumor metabolism is highly inhomogeneous within the same tumor and our study was limited by only single biopsies, taken from the contrast-enriched part. In contrast to limited analysis by metabolomic/transcriptomic data alone, an integrative approach of metabolomic-transcriptomic features allowed to capture all subgroups being described by Patel et al. This is attributed to the fact, that integrative approach mirrors the intracellular functional processes better than transcriptomic profiling alone.

All these findings confirm a strong link between metabolic alterations and specific expression subclasses. Known metabolic alterations (IDH1/2 mutation) are able to reshape the tumor methylome and may also affect specific traits of tumor gene expression. For the first time, our study described metabolic alterations, which were related to specific expression subtypes of glioblastoma multiforme. Therefore, our results point to the fact that metabolic alterations may drive the appearance of specific glioblastoma multiforme expression subclasses as a specific adaption upon microenvironmental conditions.

Finally, these findings highlight the strong mutual dependence of transcriptional oncogenic alterations and glioblastoma metabolism. The global network, which described the landscape of transcriptional/metabolomics alterations (Figure 5) revealed lots of new findings beyond the reported metabolic connections. These GBM specific alterations could be used to explore further new therapeutic strategies.

Limitations of the study

This study has some limitations. First, the small sample size and the disregarded heterogeneity could have lead to false-positive associations or confounder effects. Second, samples were guided by intraoperative neuronavigation without correction for brain shift. H&E stainings confirmed the occurrence of tumor in each biopsy. Additionally, tissue samples also contained non-tumor cells (1-10%), which could have an impact on global metabolomics/transcriptomic profiling. Conservative statistical methods with corrections for multiple testing at each level of analysis were applied. Only corrected p-values (Bonferroni, Benjamini-Hochberg) were reported for the sake of robustness. Nevertheless, these findings have to be confirmed in a larger cohort of patients.

Conclusion

This study displayed the landscape of metabolic-transcriptomic alterations in glioblastoma multiforme. Our results point to the fact that metabolic alterations may drive the appearance of specific glioblastoma multiforme expression subclasses. Additionally, we highlighted the association between metabolism and hallmarks of tumorigenesis such as cell-cycle dysfunctions or immune escape mechanisms. In conclusion, the strong influence of metabolic alterations in the wide scope of oncogenic transcriptional alterations.

MATERIALS AND METHODS

Workflow and study design

An illustration of the following workflow is shown in Figure 1. A detailed description of each analytic step is presented below.

Patients

For this prospective study we included 48 patients with primary glioblastoma multiforme WHO grade IV (without known lower-grade leasion in the patients history), who underwent surgery at the Department of Neurosurgery of the Medical Center, University of Freiburg between 2012 and 2016. The local ethics committee of the University of Freiburg approved data evaluation, imaging procedures and experimental design (protocol 100020/09 and 5565/15). The methods were carried out in accordance with the approved guidelines. Written informed consent was obtained from all patients.

Tissue collection and histology

Tumor tissue was sampled from contrast enhancing regions identified by intraoperative neuronavigation (Cranial Map Neuronavigation Cart 2, Stryker, Freiburg, Germany) during tumor resection. The tissue was snap-frozen in liquid nitrogen immediately after resection and processed for further genetic/metabolic analysis. Tissue samples were fixed using 4% phosphate buffered formaldehyde and paraffin-embedded with standard procedures. H&E staining was performed on 4 µm paraffin sections using standard protocols. Immunohistochemistry was applied using an autostainer (Dako) after heat-induced epitope retrieval in citrate buffer. IDH1 mutation was assessed by immunohistochemistry using an anti-IDH1-R123H antibody (1:20, Dianova).

1H-NMR spectroscopy

Frozen tissue samples (n=48) were extracted by methanol-water (M/W) extraction. The M/W extraction was performed as described in Beckonert et al [35]. Extracts were disintegrated by sonication. Half of the extract was used to extract DNA and determine its total concentration (for normalization step-1, as described below). Only 38 samples achieved the quality criteria. The hydrophilic-phase was separated, lyophilized and resolved in deuteriumoxide (D2O). 600μl suspension was transferred to NMR-tubes for further NMR procedures. 1H-NMR was performed at the Institute of Physical Chemistry of the University of Freiburg. 1D-NMR spectra were performed on a Bruker Avance III HDX 600-MHz FT-NMR spectrometer (Rheinstetten, Germany), equipped with the probe: PABBO BB/19F-1H/D Z-GRD. Each single spectrum within the spectra was recorded with 2 dummy scans and 32 scans with 64k points in the time domain. The sweep width was set to 16.02 ppm with an offset of 4.691 ppm. This resulted in an acquisition time of 3.4 seconds for each scan with a dwell time of 52 microseconds. The relaxation delay was set to 2 seconds, so that the total acquisition time of each spectrum was 3 minutes and 5 seconds. Water-suppressed 1H NMR spectra were acquired by using a zgesgp sequence [36]. The FID was Fourier transformed and automatically phase corrected without any further zero filling or apodization. The spectra were manually calibrated by setting the peak of L-lactate acid at 1.310 ppm. All acquisition and processing of the spectra was performed with TopSpin 3.2 patchlevel 6. A total number of 29 spectra achieved high quality, which is necessary for complex post-processing.

Post-processing and evaluation

Raw data was analyzed by “batman”, a R-software based tool for metabolite detection in complex spectra. The batman software fits a predefined list of metabolites by a Bayesian approach. A detailed description of the batman algorithm was given by Hao et al [37]. Four spectra did not archive quality criteria and were excluded from further analysis. Metabolite (spectral) intensities were used as further input for the integrative analysis.

Data normalization

Data was normalized by a two-step normalization-approach. First, raw intensity values were normalized by the given DNA concentration of each patient to balance different tumor mass of each extract. Second, normalization was performed by a median-based normalization-algorithm.

Genome-wide expression analysis

RNA of 48 patients was prepared using the RNAeasy kit (Qiagen). An amount of 1.5 µg RNA was obtained for expression arrays analysis. 8 samples did not achieve RIN>8.5 and could not used for further transcriptomic profiling. Arrays were performed by human gene ST 2.0 chip (Affymetrix). Raw data were processed, normalized and controlled by R software and the “affy” R-package. Different expression analysis and statistical testing (pairwise t-test) were performed by limma R-package.

Identification of expression-subgroups

A consensus cluster analysis was performed by the implemented R package “ConsensusClusterPlus”, with a kmax of 10. The samples were dominantly split into 3 Clusters, which was also revealed by an unsupervised hierarchical clustering (Figure 3A3B). Additionally, a random forest algorithm was used to identify the Verhaak expression subgroups for each sample (R-package: “randomForest”) [19]. Expression data of 484 TCGA samples were processed at level-3 and analyzed in the forest model as training data. The model showed a specificity of 94.3% and a sensitivity of 89.3%. Expression subgroups were predicted in the expression data of 48 patients based on the trained algorithm [38].

Weighted gene co-expression network analysis and gene set enrichment analysis

The WGCNA analysis is a robust tool for integrative network analysis and was used in recent studies [39-41]. It is based on a scaled-topology-free based network approach and uses the topological overlapping measurement to identify corresponding modules as shown in Figure 1. These modules were analyzed by their eigengene correlation to each metabolite. The correlation of the intramodule connectivity (kME) and metabolites was used as input for a “Cluster of Clusters analysis”. This analysis integrates expression modules and metabolites, which present equal correlation values (kME and metabolite intensity values). A detailed description of WGCNA is given in Heiland et al. 2016 [33].

Gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA)

A permutation-based pre-ranked Gene Set Enrichment Analysis (GSEA) was applied to each module to verify its biological functions and pathways [42]. The predefined gene sets of the Molecular Signature Database v5.1 were taken. Enrichment score was calculated by the rank order according to normalized intramodule connectivity of each gene in the expression module [42]. For significant enrichment, p-values were adjusted by FDR. Gene Set Variation Analysis (GSVA) was performed by the GSVA package implemented in R-software. The analysis based on a non-parametric unsupervised approach, which transformed a classic gene matrix (gene-by-sample) into a gene set by sample matrix resulted in an enrichment-score for each sample and pathway [43].

Statistic analysis

Normalized intensity values were clustered in an unsupervised hierarchical clustering. Cluster analysis was performed in r-software and affiliated packages. The Kaplan–Meier method was used to provide median point estimates and time-specific rates. The Hazard-Ratio (HR) was calculated by Cox-Regressions tests.

ACKNOWLEDGMENT

We thank the MagRes, University of Freiburg, Germany for coordinating the NMR usage.

CONFLICTS OF INTEREST

The authors declare no conflicts of interest.

REFERENCES

1. Crocetti E, Trama A, Stiller C, Caldarella A, Soffietti R, Jaal J, Weber DC, Ricardi U, Slowinski J, Brandes A, RARECARE working group. Epidemiology of glial and non-glial brain tumours in Europe. Eur J Cancer. 2012; 48: 1532–42. doi: 10.1016/j.ejca.2011.12.013.

2. Ostrom QT, Bauchet L, Davis FG, Deltour I, Fisher JL, Langer CE, Pekmezci M, Schwartzbaum JA, Turner MC, Walsh KM, Wrensch MR, Barnholtz-Sloan JS. The epidemiology of glioma in adults: a &quot;state of the science&quot; review. Neuro Oncol. 2014; 16: 896–913. doi: 10.1093/neuonc/nou087.

3. Ostrom QT, Gittleman H, Liao P, Rouse C, Chen Y, Dowling J, Wolinsky Y, Kruchko C, Barnholtz-Sloan J. CBTRUS Statistical Report: Primary Brain and Central Nervous System Tumors Diagnosed in the United States in 2007-2011. Neuro Oncol. 2014; 16 Suppl 4: iv1-iv63. doi: 10.1093/neuonc/nou223.

4. Gilbert MR, Dignam JJ, Armstrong TS, Wefel JS, Blumenthal DT, Vogelbaum MA, Colman H, Chakravarti A, Pugh S, Won M, Jeraj R, Brown PD, Jaeckle KA, et al. A randomized trial of bevacizumab for newly diagnosed glioblastoma. N Engl J Med. 2014; 370: 699–708. doi: 10.1056/NEJMoa1308573.

5. Chinot OL, Wick W, Cloughesy T. Bevacizumab for newly diagnosed glioblastoma. N Engl J Med. 2014; 370: 2049. Available from http://www.ncbi.nlm.nih.gov/pubmed/24860870.

6. Taal W, Oosterkamp HM, Walenkamp AM, Dubbink HJ, Beerepoot LV, Hanse MC, Buter J, Honkoop AH, Boerman D, de Vos FY, Dinjens WN, Enting RH, Taphoorn MJ, et al. Single-agent bevacizumab or lomustine versus a combination of bevacizumab plus lomustine in patients with recurrent glioblastoma (BELOB trial): a randomised controlled phase 2 trial. Lancet Oncol. 2014; 15: 943–53. doi: 10.1016/S1470-2045(14)70314-6.

7. Vredenburgh JJ, Desjardins A, Herndon JE, Marcello J, Reardon DA, Quinn JA, Rich JN, Sathornsumetee S, Gururangan S, Sampson J, Wagner M, Bailey L, Bigner DD, et al. Bevacizumab plus irinotecan in recurrent glioblastoma multiforme. J Clin Oncol. 2007; 25: 4722–9. doi: 10.1200/JCO.2007.12.2440.

8. Friedman HS, Prados MD, Wen PY, Mikkelsen T, Schiff D, Abrey LE, Yung WKA, Paleologos N, Nicholas MK, Jensen R, Vredenburgh J, Huang J, Zheng M, et al. Bevacizumab alone and in combination with irinotecan in recurrent glioblastoma. J Clin Oncol. 2009; 27: 4733–40. doi: 10.1200/JCO.2008.19.8721.

9. Van Meir EG, Hadjipanayis CG, Norden AD, Shu H-K, Wen PY, Olson JJ. Exciting New Advances in Neuro-Oncology. CA Cancer J Clin. 2010; 60: 166–93. doi: 10.3322/caac.20069.

10. Pavlova NN, Thompson CB. The Emerging Hallmarks of Cancer Metabolism. Cell Metab. Elsevier Inc. 2016; 23: 27–47. doi: 10.1016/j.cmet.2015.12.006.

11. Louis DN, Perry A, Reifenberger G, von Deimling A, Figarella-Branger D, Cavenee WK, Ohgaki H, Wiestler OD, Kleihues P, Ellison DW. The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary. Acta Neuropathol. 2016; 131: 1–18. doi: 10.1007/s00401-016-1545-1.

12. Turcan S, Rohle D, Goenka A, Walsh LA, Fang F, Yilmaz E, Campos C, Fabius AW, Lu C, Ward PS, Thompson CB, Kaufman A, Guryanova O, et al. IDH1 mutation is sufficient to establish the glioma hypermethylator phenotype. Nature. 2012; 483: 479–83. doi: 10.1038/nature10866.

13. Noushmehr H, Weisenberger DJ, Diefes K, Phillips HS, Pujara K, Berman BP, Pan F, Pelloski CE, Sulman EP, Bhat KP, Verhaak RGW, Hoadley KA, Hayes DN, et al. Identification of a CpG island methylator phenotype that defines a distinct subgroup of glioma. Cancer Cell. 2010; 17: 510–22. doi: 10.1016/j.ccr.2010.03.017.

14. Brat DJ, Verhaak RG, Aldape KD, Yung WKA, Salama SR, Cooper LAD, Rheinbay E, Miller CR, Vitucci M, Morozova O, Robertson AG, Noushmehr H, Laird PW, et al. Comprehensive, Integrative Genomic Analysis of Diffuse Lower-Grade Gliomas. N Engl J Med. 2015; 372: 2481–98. doi: 10.1056/NEJMoa1402121.

15. Lehnhardt FG, Röhn G, Ernestus RI, Grüne M, Hoehn M. 1H-and 31P-MR spectroscopy of primary and recurrent human brain tumors in vitro: Malignancy-characteristic profiles of water soluble and lipophilic spectral components. NMR Biomed. 2001; 14: 307–17. doi: 10.1002/nbm.708.

16. Hagberg G, Burlina AP, Mader I, Roser W, Radue EW, Seelig J. In vivo proton MR spectroscopy of human gliomas: definition of metabolic coordinates for multi-dimensional classification. Magn Reson Med. 1995; 34: 242–52.

17. Cuperlovic-Culf M, Ferguson D, Culf A, Morin P, Touaibia M. 1H NMR metabolomics analysis of glioblastoma subtypes: Correlation between metabolomics and gene expression characteristics. J Biol Chem. 2012; 287: 20164–75. doi: 10.1074/jbc.M111.337196.

18. Marziali G, Signore M, Buccarelli M, Grande S, Palma A, Biffoni M, Rosi A, D’Alessandris QG, Martini M, Larocca LM, De Maria R, Pallini R, Ricci-Vitiani L. Metabolic/Proteomic Signature Defines Two Glioblastoma Subtypes With Different Clinical Outcome. Sci Rep. 2016; 6: 21557. doi: 10.1038/srep21557.

19. Verhaak RG, Hoadley KA, Purdom E, Wang V, Qi Y, Wilkerson MD, Miller CR, Ding L, Golub T, Mesirov JP, Alexe G, Lawrence M, O’Kelly M, et al. Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell. 2010; 17: 98–110. doi: 10.1016/j.ccr.2009.12.020.

20. Phillips HS, Kharbanda S, Chen R, Forrest WF, Soriano RH, Wu TD, Misra A, Nigro JM, Colman H, Soroceanu L, Williams PM, Modrusan Z, Feuerstein BG, et al. Molecular subclasses of high-grade glioma predict prognosis, delineate a pattern of disease progression, and resemble stages in neurogenesis. Cancer Cell. 2006; 9: 157–73. doi: 10.1016/j.ccr.2006.02.019.

21. Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, Cahill DP, Nahed BV, Curry WT, Martuza RL, Louis DN, Rozenblatt-Rosen O, Suvà ML, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science (80-). 2014; 344: 1396–401. doi: 10.1126/science.1254257.

22. Liang C, Qin Y, Zhang B, Ji S, Shi S, Xu W, Liu J, Xiang J, Liang D, Hu Q, Ni Q, Xu J, Yu X. Metabolic plasticity in heterogeneous pancreatic ductal adenocarcinoma. Biochim Biophys Acta. 2016; 1866: 177–88. doi: 10.1016/j.bbcan.2016.09.001.

23. Garner JM, Ellison DW, Finkelstein D, Ganguly D, Du Z, Sims M, Yang CH, Interiano RB, Davidoff AM, Pfeffer LM. Molecular heterogeneity in a patient-derived glioblastoma xenoline is regulated by different cancer stem cell populations. PLoS One. 2015; 10: e0125838. doi: 10.1371/journal.pone.0125838.

24. Haukaas TH, Euceda LR, Giskeødegård GF, Lamichhane S, Krohn M, Jernström S, Aure MR, Lingjærde OC, Schlichting E, Garred Ø, Due EU, Mills GB, Sahlberg KK, et al. Metabolic clusters of breast cancer in relation to gene- and protein expression subtypes. Cancer Metab. 2016; 4: 12. doi: 10.1186/s40170-016-0152-x.

25. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008; 9: 559. doi: 10.1186/1471-2105-9-559.

26. Horgan RP, Kenny LC. SAC review “Omic” technologies: proteomics and metabolomics. Obstet Gynaecol. 2011; 13: 189–95. doi: 10.1576/toag.13.3.189.27672.

27. Chinnaiyan P, Kensicki E, Bloom G, Prabhu A, Sarcar B, Kahali S, Eschrich S, Qu X, Forsyth P, Gillies R. The metabolomic signature of malignant glioma reflects accelerated anabolic metabolism. Cancer Res. 2012; 72: 5878–88. doi: 10.1158/0008-5472.CAN-12-1572-T.

28. Oppermann H, Ding Y, Sharma J, Berndt Paetz M, Meixensberger J, Gaunitz F, Birkemeyer C. Metabolic response of glioblastoma cells associated with glucose withdrawal and pyruvate substitution as revealed by GC-MS. Nutr Metab (Lond). 2016; 13:1–11. doi: 10.1186/s12986-016-0131-9.

29. Moffett JR, Namboodiri MA. Tryptophan and the immune response. Immunol Cell Biol. 2003; 81: 247–65. doi: 10.1046/j.1440-1711.2003.t01-1-01177.x.

30. Kinoshita Y, Yokota A. Absolute concentrations of metabolites in human brain tumors using in vitro proton magnetic resonance spectroscopy. NMR Biomed. 1997; 10: 2–12.

31. Wu P, Jiang J, Liu Y, Hu K, Jiang WD, Li SH, Feng L, Zhou XQ. Dietary choline modulates immune responses, and gene expressions of TOR and eIF4E-binding protein2 in immune organs of juvenile Jian carp (Cyprinus carpio var. Jian). Fish Shellfish Immunol. 2013; 35: 697–706. doi: 10.1016/j.fsi.2013.05.030.

32. Talasila KM, Røsland GV, Hagland HR, Eskilsson E, Flønes IH, Fritah S, Azuaje F, Atai N, Harter PN, Mittelbronn M, Andersen M, Joseph JV, Hossain JA, et al. The angiogenic switch leads to a metabolic shift in human glioblastoma. Neuro Oncol. 2017; 19:383-393. doi: 10.1093/neuonc/now175.

33. Heiland DH, Mader I, Schlosser P, Pfeifer D, Carro MS, Lange T, Schwarzwald R, Vasilikos I, Urbach H, Weyerbrock A. Integrative Network-based Analysis of Magnetic Resonance Spectroscopy and Genome Wide Expression in Glioblastoma multiforme. Sci Rep. 2016; 6: 29052. doi: 10.1038/srep29052.

34. Ogretmen B, Hannun YA. Biologically active sphingolipids in cancer pathogenesis and treatment. Nat Rev Cancer. 2004; 4: 604–16. doi: 10.1038/nrc1411.

35. Beckonert O, Keun HC, Ebbels TMD, Bundy J, Holmes E, Lindon JC, Nicholson JK. Metabolic profiling, metabolomic and metabonomic procedures for NMR spectroscopy of urine, plasma, serum and tissue extracts. Nat Protoc. 2007; 2: 2692–703. doi: 10.1038/nprot.2007.376.

36. Hwang TL, Shaka AJ. Water Suppression That Works. Excitation Sculpting Using Arbitrary Wave-Forms and Pulsed-Field Gradients. Journal of Magnetic Resonance, Series A. 1995; 13: 275–9. doi: 10.1006/jmra.1995.1047.

37. Hao J, Liebeke M, Astle W, De Iorio M, Bundy JG, Ebbels TM. Bayesian deconvolution and quantification of metabolites in complex 1D NMR spectra using BATMAN. Nat Protoc. 2014; 9: 1416–27. doi: 10.1038/nprot.2014.090.

38. Liaw A, Wiener M. Classification and Regression by randomForest. R news. 2002; 2: 18–22. doi: 10.1177/154405910408300516.

39. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008; 9: 559. doi: 10.1186/1471-2105-9-559.

40. Holtman IR, Raj DD, Miller JA, Schaafsma W, Yin Z, Brouwer N, Wes PD, Möller T, Orre M, Kamphuis W, Hol EM, Boddeke EW, Eggen BJ. Induction of a common microglia gene expression signature by aging and neurodegenerative conditions: a co-expression meta-analysis. Acta Neuropathol Commun. 2015; 3: 31. doi: 10.1186/s40478-015-0203-5.

41. Iancu OD, Colville A, Oberbeck D, Darakjian P, McWeeney SK, Hitzemann R. Cosplicing network analysis of mammalian brain RNA-Seq data utilizing WGCNA and Mantel correlations. Front Genet. 2015; 6: 174. doi: 10.3389/fgene.2015.00174.

42. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005; 102: 15545–50. doi: 10.1073/pnas.0506580102.

43. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14: 7. doi: 10.1186/1471-2105-14-7.


Creative Commons License All site content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 License.
PII: 16544