Identification of a Metabolism-Related Risk Signature Associated With Clinical Prognosis in Glioblastoma Using Integrated Bioinformatic Analysis

Altered metabolism of glucose, lipid and glutamine is a prominent hallmark of cancer cells. Currently, cell heterogeneity is believed to be the main cause of poor prognosis of glioblastoma (GBM) and is closely related to relapse caused by therapy resistance. However, the comprehensive model of genes related to glucose-, lipid- and glutamine-metabolism associated with the prognosis of GBM remains unclear, and the metabolic heterogeneity of GBM still needs to be further explored. Based on the expression profiles of 1,395 metabolism-related genes in three datasets of TCGA/CGGA/GSE, consistent cluster analysis revealed that GBM had three different metabolic status and prognostic clusters. Combining univariate Cox regression analysis and LASSO-penalized Cox regression machine learning methods, we identified a 17-metabolism-related genes risk signature associated with GBM prognosis. Kaplan-Meier analysis found that obtained signature could differentiate the prognosis of high- and low-risk patients in three datasets. Moreover, the multivariate Cox regression analysis and receiver operating characteristic curves indicated that the signature was an independent prognostic factor for GBM and had a strong predictive power. The above results were further validated in the CGGA and GSE13041 datasets, and consistent results were obtained. Gene set enrichment analysis (GSEA) suggested glycolysis gluconeogenesis and oxidative phosphorylation were significantly enriched in high- and low-risk GBM. Lastly Connectivity Map screened 54 potential compounds specific to different subgroups of GBM patients. Our study identified a novel metabolism-related gene signature, in addition the existence of three different metabolic status and two opposite biological processes in GBM were recognized, which revealed the metabolic heterogeneity of GBM. Robust metabolic subtypes and powerful risk prognostic models contributed a new perspective to the metabolic exploration of GBM.

Altered metabolism of glucose, lipid and glutamine is a prominent hallmark of cancer cells. Currently, cell heterogeneity is believed to be the main cause of poor prognosis of glioblastoma (GBM) and is closely related to relapse caused by therapy resistance. However, the comprehensive model of genes related to glucose-, lipid-and glutamine-metabolism associated with the prognosis of GBM remains unclear, and the metabolic heterogeneity of GBM still needs to be further explored. Based on the expression profiles of 1,395 metabolism-related genes in three datasets of TCGA/CGGA/GSE, consistent cluster analysis revealed that GBM had three different metabolic status and prognostic clusters. Combining univariate Cox regression analysis and LASSO-penalized Cox regression machine learning methods, we identified a 17-metabolism-related genes risk signature associated with GBM prognosis. Kaplan-Meier analysis found that obtained signature could differentiate the prognosis of high-and low-risk patients in three datasets. Moreover, the multivariate Cox regression analysis and receiver operating characteristic curves indicated that the signature was an independent prognostic factor for GBM and had a strong predictive power. The above results were further validated in the CGGA and GSE13041 datasets, and consistent results were obtained. Gene set enrichment analysis (GSEA) suggested glycolysis gluconeogenesis and oxidative phosphorylation were significantly enriched in high-and low-risk GBM. Lastly Connectivity Map screened 54 potential compounds specific to different subgroups of GBM patients. Our study identified a novel metabolism-related gene signature, in addition the existence of three different metabolic status and two opposite biological processes in GBM were recognized, which revealed the metabolic heterogeneity of GBM. Robust metabolic subtypes and powerful risk prognostic models contributed a new perspective to the metabolic exploration of GBM.
Keywords: glioblastoma, metabolism, prognosis, signature, heterogeneity INTRODUCTION Glioblastoma (GBM) is the most common, most aggressive and worst prognosis glioma in adults (1), accounting for about 55% of gliomas, with median survival of only 14-16 months (1,2). The diffusivity and invasiveness of GBM itself and the inter-/intra-tumor heterogeneity lead to GBM therapy resistance and high recurrence (3)(4)(5)(6). Therefore, despite the standard treatment protocol for GBM such as surgical resection, radiotherapy and chemotherapy, the prognosis of GBM still remains dismal, and the 5-year survival rate is only about 5% (2,7,8). Novel molecular markers or therapeutic targets are urgently needed to improve the prognosis of GBM.
Metabolic reprogramming is one of the hallmarks of cancer cells, and there is growing evidence that metabolic dysregulation plays an important role in the growth, proliferation, angiogenesis, and invasion of cancer cells (9)(10)(11). Civita et al. (12) revealed the landscape of GBM heterogeneity using laser capture microdissection and RNA-seq analysis, which showed dysregulation of metabolic pathways, providing direct evidence for metabolic alterations in GBM. The metabolism of glucose, lipid, and glutamine in cancer cells is altered (13). According to Warburg's basic research, cancer cells obtain energy mainly through the glycolytic pathway rather than the oxidative phosphorylation (OxPhos) pathway, so abnormal glycolytic metabolism is one of the basic characteristics of malignant cells (14). Recent studies have found that lipid metabolism reprogramming plays a crucial role in membrane synthesis, energetic production and signal transduction in the progression of cancer cells (15). Nuclear magnetic resonance (NMR) spectroscopy revealed that unsaturated fatty acids, cholesterol esters and phosphatidylcholine are only present in the GBM (16,17). At present, the biological phenotype and molecular mechanism of lipid composition change leading to glioma need further study. In addition, studies on cancer cell metabolism have provided evidence that tumor-specific activation of signaling pathways, such as the upregulation of the oncogene Myc, can regulate glutamine uptake and its metabolism through glutaminolysis to provide the cancer cell with a replacement of energy source (18). Therefore, the in-depth exploration of three major metabolites of GBM may provide an important theoretical basis for the development of new treatment (10,19). However, there has been no comprehensive analysis of the three-major metabolism-related genes and their prognostic value in GBM.
In the present study, we comprehensively analyzed GBM mRNA sequencing data from three public datasets, the Cancer Genome Atlas (TCGA), The Chinese Glioma Genome Atlas (CGGA), and GSE13041, to explore the metabolic status of GBM patients. Through cluster analysis, the patients can be divided into 3 stable clusters according to the gene expression profile, and the prognosis and molecular characteristics of the 3 clusters are significantly different. In addition, we further screened potential specific therapeutic compounds for each cluster. More importantly, we identified a metabolism-related risk signature to assess the prognosis of patients with GBM in the TCGA datasets which can be served as an independent predictor closely related to the prognosis of GBM patients. And the high-and low-risk groups had distinctly different biological processes. The above results were further validated in CGGA and GSE13041 datasets. In summary, robust prognostic risk models and subtypes contribute to a better understanding of the molecular pathogenesis of GBM.

Data Collection
Whole genome mRNA expression sequencing data and corresponding clinical information [histology, subtype, gender, age, isocitrate dehydrogenase1 (IDH1) mutational status, methylguanine methyltransferase (MGMT) promoter status, glioma cytosine-phosphate-guanine island methylator phenotype (G-CIMP) status and survival information] of 165 GBM patients were downloaded from TCGA (PanCancer Atlas) (20) datasets as training set. Similarly, GSE13041 (21) and CGGA RNA expression data and clinical information were obtained as validation sets. Of which, mRNA sequencing data from CGGA contains two datasets, mRNAseq_693 and mRNAseq_325, whose platforms are Illumina HiSeq and Illumina HiSeq 2000 or 2500, respectively. Therefore, we removed the batch effect from these two datasets and normalized them to obtain an integrated data of 216 GBM patients. The characteristics of GBM patients from these three datasets were summarized in Supplementary Table 1. Five hundred and thirty two Glucose-, 1034 lipid-, and 13 glutamine-metabolism related genes were downloaded from Molecular Signature Database v7.0 (MSigDB) (http://www.broad.mit.edu/gsea/msigdb/) (22). The detailed metabolism-related genes were listed in the Supplementary Table 2.

Consensus Clustering
We took the expression profile of metabolism-related genes for consistent clustering by using "Cancersubtype" R package (http://cran.r-project.org). The euclidean distance was applied to calculate the similarity distance between samples, and K-means methods was utilized for clustering. By performing resampling analysis, 80% of the samples were sampled for 100 times. The optimal number of clusters was determined by the cumulative distribution function (CDF) and was validated in the CGGA and GSE13041 datasets. And principal component analysis (PCA) was carried out using the R package "princomp" to validate the molecular subtype.

Gene Signature Identification
Univariate Cox regression were performed on the 169 GBM patients in TCGA datasets to select the optimal prognostic gene set with R package "glmnet." After getting the corresponding hazard ratio (HR) and p-value of each gene, the genes with p < 0.05 were selected as seed genes for Cox LASSO regression with 10-fold cross-validation (CV). Risk score for each patient of the TCGA training set was calculated with the linear combinational of the signature gene expression (expr) weighted by their regression coefficients (23).
Risk score = expr gene1 ×β gene1 + expr gene2 ×β gene2 Frontiers in Oncology | www.frontiersin.org In the above equation, "N" was the total number of key genes, "expr" represented the expression value of geneN, and β represented the selected gene coefficient from LASSO analysis. Patients in the training datasets was then categorized into high and low risk score groups according to the median risk score. The risk score for each patient in the validation datasets were also calculated based on the same risk formula. The multivariate Cox regression analysis was conducted to determine whether the risk score was an independent predictor for GMB patients with R package "survival." The differences in overall survival (OS) between high-risk and low-risk score in the training and validation datasets were estimated using the Kaplan-Meier (KM) method. Receiver operating characteristic (ROC) analysis was performed to evaluate the accuracy of the risk model with R package "pROC." Gene Set Enrichment Analysis (GSEA) GSEA was performed using Gene Set Enrichment Analysis v3 software downloaded from the Broad Institute (http://www. broadinstitute.org/gsea/index.jsp) (24). The mRNA expression profile of GBM samples from the TCGA/CGGA/GSE13041 datasets were analyzed by GSEA. For GSEA, risk score was selected as a binary variable divided into high-and low-risk by a criterion of whether the score was greater than the median value. The collection of annotated gene sets of c2.cp.kegg.v7.0.symbols.gmt in MSigDB was chosen as the reference gene sets in GSEA software, the NOM p < 0.05 and false discovery rate (FDR) < 0.25 was set as the cutoff.

Connectivity Map (CMap) Analysis
CMap is a systematic, data-driven program for detecting correlations among genes, compounds, and biological conditions. We queried the recently updated CMap to screen potential compounds that might target metabolismrelated pathways. A list of differential expression genes (DEGs) among the 3 clusters in TCGA was obtained using the "lmFit" function of the R package "limma" with default parameters (25), and the top 373 genes (148 upregulated and 225 downregulated) were selected to uploaded into the CMap database. Compounds with an absolute value of enrichment ≥ 0.7 and p < 0.05 were selected as potential therapeutic drugs for GBM.

Statistical Analysis
One-way ANOVA was performed to compare the differences of risk score between/among subtypes/clusters. The one-way ANOVA method and Tukey's test was applied to identify the DEGs between the high-and low-risk groups (q < 0.05, |log2FC| > 2). KM curve with log-rank test was used to assess the OS differences among/between different groups. Univariate and multivariate Cox regression analyses were conducted to assess the independent prognostic factors. The statistical analyses were conducted using R software version 3.5.1 (R Core Team, R Foundation for Statistical Computing, Vienna, Austria), p < 0.05 was regarded as statistically significant.

Molecular Cluster Identification and Validation
The intersection of the three datasets and metabolism genes was extracted, and the overlapping genes were removed. The gene expression profiles of 1,395 metabolism-related genes (Supplementary Table 2) were exploited to identify the GBM clusters in TCGA cohort. All GBM samples were grouped into k (k = 2, 3, 4, 5, 6, 7, 8, 9) different subtypes using "Cancersubtype" R package. According to the CDF curves of the consensus score, we selected the k = 3 as the optimal division (Figures 1A-C and Supplementary Figure 1). KM analysis showed that the OS of the 3 clusters were significantly different ( Figure 1D, p < 0.05), and PCA revealed that the 3 clusters could be separated from each other ( Figure 1E). Figure 1F showed the heatmap of these 3 clusters defined by the top 100 variable expression genes. The above results indicated that there were significant metabolic phenotypes among the 3 clusters. In order to validate the stability of molecular subtypes, we further selected CGGA and GSE13041 datasets for clustering. The clustering results of molecular subtypes in CGGA and GSE13041 datasets were consistent with those in TCGA, and the relevant results were shown in Supplementary Figures 2-5, respectively. Therefore, we identified three stable clusters of GBM based on the expression of metabolism-related genes.

Identification of Metabolism-Related Genes Signature for Prognostic Prediction
The 1,395 putative metabolism-related genes were exploited to conducting univariate cox regression analysis. Firstly, we identified 29 significantly metabolism-related genes associated with the survival of GBM with p < 0.05 (Supplementary Table 2). We further performed LASSO Cox regression algorithm with cross-validation (Figures 2A,B), after 1,000-time iterations, a 17-gene risk signature was constructed ( Table 1) and the risk score for each patient was calculated with their expression level and regression coefficient. To comprehensively investigate the relationship between risk score and patients' survival, we further stratified patients into highand low-risk groups based on the median risk score ( Figure 2C). And as shown in Figure 2C, patients of GBM in TCGA cohort with high risk score have more death cases when compare to low risk score. KM curves analysis result revealed that patients in high risk group had a shorter survival time than in the low risk group ( Figure 2D). To validate this gene set, we also calculated patients' risk scores of CGGA and GSE13041 cohorts with same regression coefficient. And as expected, consensus result was also obtained in the validation datasets (Supplementary Figure 6).

The 17-Gene Risk Signature Shows Strong Prognostic Power
Univariate and multivariate Cox regression analyses were performed to determine prognostic factors for survival in TCGA/CGGA/GSE13041 patients with GBM. As shown in Table 2, the 17-gene risk signature was independently correlated Frontiers in Oncology | www.frontiersin.org with OS (p < 0.05), and the genes signature could also be served as an independent prognostic factor in the CGGA and GSE13041 cohorts (p < 0.05).
In order to further clarify the significance of risk signature in GBM, we detected the correlation between risk score and major clinical features. The risk score distribution of patients with different IDH1 status, MGMT promotor methylation status and G-CIMP status was significantly different in TCGA cohort. The risk scores were lower in GBM patients with IDH1-mutant type, MGMT promoter methylated and G-CIMP subgroups than the IDH1-wild type (wt), MGMT promoter unmethylated and non G-CIMP ones of TCGA cohort, respectively (Figures 3A-C). Consistent results could also be observed in the CGGA cohort, the risk scores of IDH1-mutant type and MGMT promoter methylated subgroups were lower than the IDH1wt and MGMT promoter unmethylated ones, respectively (Figures 3D,E). However, there was no significant difference in risk score between the 1p/19q codel and the non-codel groups, which may be due to the small sample size of the 1p/19q codel group ( Figure 3F). In addition, no statistically significant differences were observed between risk scores for different age stratifications and for different molecular subtypes (Supplementary Figure 7). We further compared the differences in risk scores among three different clusters we identified above in TCGA/CGGA/GSE13041. Interestingly, as shown in the Figures 3G-I, the cluster with the worst prognosis had the highest risk score (cluster2 in TCGA and GSE13041, and cluster1 in CGGA), while the cluster with the best prognosis had the lowest risk score (cluster1 in TCGA and GSE13041, cluster3 in CGGA).
We further validated the prognostic predictive power of the risk signature we identified in different clinical feature stratified groups, such as IDH1-wt/-mutant, MGMT promoter methylated/unmethylated, G-CIMP/non G-CIMP, and 1p/19q codel/non-codel cohorts. In the TCGA and CGGA cohorts, KM analysis showed that cases with high-risk score had shorter OS than the low-risk ones in most stratified patients (Figure 4). However, this result was not observed in the TCGA_IDH1mutant group, TCGA_G-CIMP group and CGGA_ 1p/19q-codel group due to the small sample size of cases.
By performing the ROC analysis in the training datasets and validation datasets, we next evaluate the accuracy of the risk signature. The result showed that the AUC value of 1-, 2-, and 3-year for the TCGA datasets were 0.710, 0.783, and 0.873, respectively ( Figure 5A). And in the CGGA/GSE validation sets, a similar strong prognostic power was obtained (Figures 5B,C). In addition, we also compared the accuracy of the clinical features and risk score for the survival prediction   of GBM and discovered that our risk signature is the optimal (Figures 5D-F). These results all confirmed that the metabolism risk signature we constructed had a strong prognostic prediction power.

The High-and Low-Risk Groups Present Different Biologic Processes
In order to explore the difference in biological process between the high-and the low-risk group in the TCGA cohort, GSEA was performed to compare the gene expression of patients in the two groups. PCA showed that in the three public cohorts of TCGA/CGGA/GSE13041 cases in the high-risk and low-risk groups were significantly distributed in two regions (Supplementary Figure 8). As shown in the Figure 6, GSEA indicated that glycolysis gluconeogenesis, WNT signaling pathway, pathways in cancer, MAPK signaling pathway, ERBB signaling pathway, adherens junction, focal adhesion, ECM receptor interaction and tight junction were significantly enriched in high-risk patients, while low-risk cases showed enrichment of OxPhos. The significantly pathways were selected based on the screening criteria of nominal p < 0.05 and false discovery rate (FDR) < 0.25. The results indicated that there were significant differences in biological processes between the two groups of GBM with different metabolic risk levels.

CMap Analysis Identifies Novel Candidate Compounds Targeting the GBM Clusters
To identify potential compounds capable of targeting the pathways associated with metabolism-related genes, we queried the CMap database using the mRNA expression signatures by applying differential expression analysis to GBM subgroups samples. The 54 compounds that were able to repress the above gene expression profile of GBM are shown in Figure 7. CMap mode of action (MoA) analysis of the 54 compounds revealed 44 mechanisms of action shared by the above compounds. Three compounds (fludroxycortide, fluorometholone, and hydrocortisone) shared the MoA of glucocorticoid receptor agonist, and two compounds (physostigmine and skimmianine) shared the MoA of acetylcholinesterase inhibitor. Moreover, a total of 14 compounds shared the following 7 mechanisms: adrenergic receptor agonists, inhibitors of bacterial cell wall synthesis, carbonic anhydrase inhibitors, DNA synthesis inhibitor, dopamine receptor antagonist, and phosphodiesterase inhibitors.

DISCUSSION
In the present study, through the comprehensive analysis of genes involved in glucose, lipid, and glutamine metabolism of GBM patients in three datasets, the risk signature closely related to the prognosis of GBM was constructed, and the heterogeneity of GBM metabolism was further revealed. Firstly, leveraging a large cohort of GBM profile, 3 clusters of GBM with different clinical features were successfully obtained through the cluster analysis. Subsequent analysis revealed that the 3 clusters had significantly different risk scores. In addition, we identified a 17 metabolism-related gene risk signature as a significant independent predictor for the prognosis of GBM by univariate cox regression analysis and LASSO analysis. GSEA suggested glycolysis gluconeogenesis and OxPhos were significantly enriched in high-and low-risk GBM. Moreover, we further demonstrated the robustness of molecular subtype and the predictive power of 17 metabolism-related gene risk signature in two validation datasets, CGGA and GSE13041, and achieved consistent results. Finally, we used CMap database to screen compounds that may target metabolism-related genes, and it is hoped that targeted therapy can be performed on GBM clusters with different metabolic status.
Considering that univariate Cox model has insufficient dimensional data on variable selection, we first performed univariate Cox model to obtain the genes related to overall survival and applied Cox LASSO regression to improve the performance index for predicting prognosis (26). None of the 17 genes showed high coefficients in the Cox model, but the cumulative effect of the 17 genes signature on OS obtained the optimal survival prediction. Interestingly, in the 17-metabolism gene risk signature we constructed, it has been reported that FOXO3 acetylation plays a central role in the regulation of glycolytic metabolism in glioblastoma, and the survival of GBM patients with FOXO3 acetylation is shorter (27). In addition, FOXO3 is related to GBM temozolomide (TMZ) resistance, and the phosphorylated AKT/FOXO3 axis regulates the expression of long non-coding RNA related to TMZ resistance GBM cells (28). Another gene, PIK3R1, is part of the RTK/ RAS/(3)K signaling pathway, which is mutated in many cancers and plays a key role in the proliferation, differentiation, and survival of cancer cells. Nearly 90% of normal GBMs showed different degrees of PIK3R1 changes, leading to abnormal activation of RTK/RAS/PI(3)K signal cascade (29). PIK3R1 has been found to promote the transformation of malignant astrocytes into glioma-like state (30). Stratifin (SFN, 14-3-3 sigma), as an oncogene related to cell proliferation, facilitates the development and progression of a variety of cancers (31,32) including gliomas (33) and has the potential to be a new therapeutic target. In addition, APOD has been identified to be associated with astrocytoma progression (34). And SH3GLB1, as the autophagy-related gene, is associated with glioma prognosis. Knockdown of SH3GLB1 inhibits glioma cell proliferation, migration and invasion, and improves sensitivity to temozolomide (35). In addition, some of the genes in the 17-metabolism risk signature are involved in the development of a variety of cancers. NR1H4, also known as farnesoid X receptor (Fxr), is the gene with the highest positive coefficient and hazard ratio in the metabolism-related risk signature. Previous studies have shown that NR1H4 plays an important role in the development of colon cancer by regulating the stability of c-Myc (36). KLF15 has been reported to inhibit cell growth in lung adenocarcinoma (37), gastric cancer (38), and colorectal cancer (39), and can be used to predict prognosis. High ADRA2A expression was associated with poor overall survival for breast (40) and bladder cancer (41). As a genetic marker, RNASEL has been linked to lethal prostate cancer (42). And ESRRB (or ERRβ) is a negative regulator of cell cycle and may be a therapeutic target for breast cancer (43). Inhibition of HSPH1 downregulates the expression of Bcl-6 and c-Myc and hampers the growth of human aggressive B cell non-Hodgkin lymphoma (44). ACADS, as one of the key metabolic genes related to the metabolic response involved in carcinogenesis, is regulated by DNA methylation and can be used as a potential methylation biomarker related to the proliferation and metastasis of hepatocellular carcinoma (45). RUFY1, named RUN, and FYVE domain containing 1, is a member of RUFY family (46). RUFY1 regulates the transport of integrins and participates in the migration of NIH-3T3 fibroblasts (47). Evidence shows that RUFY1, as a tumor promoter gene, plays an important role in the development of gastric cancer. Targeting PODXL/RUFY1 complex may improve the prognosis of gastric cancer (GC) and provide new treatment opportunities for GC patients (48). Studies have shown that PCSK1 is the gene most significantly associated with poor response to concurrent chemoradiotherapy (CCRT) in rectal cancer. Therefore, overexpression of PCSK1 is one of the risks of poor CCRT response and prognosis in rectal cancer patients (49). And PCSK1 is also over expressed in pure fibrolamellar hepatocellular carcinoma as one of neuroendocrine genes (50). As for ALAS gene, in non-small-cell lung cancer, nonerythrocyte ALAS1 gene transcription levels and ALAS1 protein levels were significantly elevated in cancer cells, while ALAS2 transcription levels were increased nearly 5-fold in HCC4017 cells (51). However, 2 of the 17-metabolism risk genes, including SPTSSA and ARSF have not been studied in cancers. SPTSSA, also known as serine palmityl transferase small subunit A (SPTSSA), encodes A subunit of SPT that is a rate-limiting enzyme in the sphingolipid biosynthesis pathway (52). As a component of eukaryotic membranes, SLs has a variety of critical functions in the growth and development of embryos and the maintenance of normal physiology (53). In a clear cell renal cell carcinoma (ccRCC) study, bioinformatics analysis revealed that 10 genes, including SPTSSA, significantly coregulated ccRCC with SPTLC1, and the low expression of SPTLC1 was significantly correlated with the disease progression and poor prognosis of ccRCC (54). In addition, another gene ARSF, which is located in Xp22.3 with ARSD and ARSE, has significant homology and highly similar intron/exon structure. Meanwhile, the splicing sites of ARSF, ARSD, ARSE, and ARSC are all at the same position. The biological role of ARSF may therefore be masked by the other three genes (55). The biological role of 17 metabolismrelated genes in glioblastoma remains to be further explored.
The Warburg effect has been widely described in GBM and other tumor types. According to the Warburg hypothesis (14), cancers are partly caused by impaired mitochondrial function and OxPhos, which is characterized by cancer cells producing most of their energy through glucose fermentation (such as aerobic glycolysis), with limited OxPhos capacity (56). This metabolic reprogramming is thought to be an adaptive mechanism for the rapid growth of tumor cells to meet their increasing energy needs. However, the intrinsic cellular heterogeneity of GBM raises the question of whether the survival and proliferation of different cell subsets is limited to glucose fermentation or other metabolic pathways. Recent studies suggest that the residual activity of mitochondrial function in GBM cells can still provide OxPhos for cancer cells (57)(58)(59). According to Deleyrolle LP's study (6), subgroups of cells with different metabolic requirements exist in GBM, in which fast-cycling cells utilize aerobic glycolysis, while slow-cycling cells preferentially utilize mitochondrial OxPhos to obtain metabolism energy. But some studies hold the opposite view that GBM tissue are unable to obtain significant energy from OxPhos. Because ultrastructural and biochemical evidence suggests that GBM cells exhibit defects in the number, structure and function of mitochondria, thus incompatible with OxPhos energy production (60)(61)(62)(63). In addition, large numbers of mitochondrial DNA (mtDNA) mutations have been found in 13 cancers including GBM that compromise OxPhos function (64). Therefore, emerging evidence indicates that cancer cells including GBM obtain energy through the glutaminolysis pathway using mitochondrial substrate-level phosphorylation (mSLP) as an alternative to OxPhos (61). Moreover, some studies based on the above theory, ketogenic metabolic therapy, as an alternative standard of care, has the potential to improve outcomes for GBM patients and other malignant brain cancers, and has yielded impressive results in clinical practice for GBM treatment (65)(66)(67)(68)(69). Here, functional analysis in our study indicated that GBMs are metabolic heterogeneity and the biological process differences between the high-and the low-risk group mainly focused on the metabolic patterns and signaling pathways. Glycolysis gluconeogenesis and OxPhos were significantly enriched in high-and low-risk GBMs. Whether OxPhos is caused by mSLP remains to be confirmed by future studies. In addition, our functional enrichment analysis identified metabolic signaling pathways associated with high-risk GBM, including WNT signaling pathway, MAPK signaling pathway, ERBB signaling pathway and pathways in cancer. Of which, abnormal Wnt signaling is recognized to drive metabolic alterations such as glycolysis, lipogenesis and glutaminolysis, which are critical to the survival of cancer stem cells (70). The pathogenesis of GBM involves multiple levels of WNT signal pathway, including tumorigenesis, stem cell maintenance, invasion, and drug resistance. Inhibition of WNT signal transduction is expected to be a new direction of GBM therapy (70,71).
Interestingly, our results suggest that some conventional antipsychotics among drugs mentioned above, including amitriptyline, fluoxetine, and imipramine, may also play a surprising role in the treatment of GBM, which is consistent with previous studies and offers new hope for patients with GBM (12,75,76,80,83,84,104). In addition, Leite et al. (105) demonstrated that clomipramine (tricyclic antidepressant drugs such as imipramine) has an impact on GBM growth and has no toxicity in normal cells (astrocytes and microglia). Reaserches have further elucidated the potential mechanism of conventional antidepressants therapy for GBM. By targeting the respiratory chain complex III and changing membrane potential, the tricyclics antidepressants mentioned above induce mitochondria-mediated apoptosis of malignant glioma cells, activate the intrinsic pathway of cytochrome-C release and caspase-3 dependent apoptosis process, and finally results in glioma cell death (106)(107)(108)(109)(110)(111)(112). It is worth mentioning that previous studies have also indicated that the anticoagulant ticlopidine synergized with imipramine to induce the death of human glioma cell lines and primary mouse glioma cells (98,99). In addition, some of the other drugs in our findings mentioned above have been reported as a combination of glioblastoma chemotherapy, such as rollipram, a promising immunotherapeutic adjuvant that can potentiate the effect of bevacizumab on GBM (89,90). The antiangiogenic effect of rofecoxib makes GBM chemotherapy more effective (96). Tacrolimus endows the GBM stem-like cells chemosensitive to the MRP1 drug substrate (97). Of course, it is still controversial whether some of the drugs identified by our study are effective against GBM. Such as vorinostat combined with bevacizumab significantly inhibited angiogenesis of GBM stem cells in vitro (100), but no significant benefit was observed in patients with GBM in clinical trials (113). However, our study provides evidence to explore the underlying mechanism of these drugs in the treatment of GBM, and it may be a new direction to investigate whether these drugs have the potential to treat GBM from a metabolic perspective.
However, some limitations of our present study should be acknowledged. First, GSE13041, one of the external validation datasets which was selected in this study lacked some clinical information, so the prognostic value of signatures needed to be validated by more external datasets, and multicenter prospective studies are needed to assess the feasibility of risk gene signatures in the future. Second, the deep molecular mechanisms behind our findings need to be further elucidated in experimental studies to facilitate our understanding of the functional roles and clinical applications of metabolism-related genes and potential molecular compounds. In particular, the issue of whether GBM tumor cells can utilize Oxphos (or mSLP) to generate the required energy remains to be further explored due to the lack of support from relevant datasets.

CONCLUSION
Our study identified a novel metabolism-related gene signature, in addition the existence of three different metabolic status and two opposite biological processes in GBM were recognized, which revealed the metabolic heterogeneity of GBM. Robust metabolic subtypes and powerful risk prognostic models contributed a new perspective to the metabolic exploration of GBM.

AUTHOR CONTRIBUTIONS
GL and HX conceived and designed the study. RZ put forward constructive suggestions to the study. ZH and CW analyzed the results, completed the visualization of the image, and contributed equally to this work. ZH completed the manuscript. All authors participated in the preparation approved the final manuscript.

ACKNOWLEDGMENTS
ZH sincerely thanks his family for their help and support, especially his parents, wife and son for providing him with a quiet and undisturbed working environment. At the same time, due to the excellent performance of Apple MacBook Pro laptop, the whole data analysis process was very pleasant and enjoyable.