We downloaded the transcriptome profiles (HTSeq Fragments Per kilobase of transcript per Million mapped reads [FPKM]) and their associated clinical profiles from TCGAbiolinks (Colaprico et al., 2016), which was used as The Cancer Genome Atlas (TCGA) training cohort (n = 501) to establish the risk score signature. The IDH1 (encoding isocitrate dehydrogenase 1) statuses were calculated using maftools (Mayakonda et al., 2018). The first validation cohort comprised RNA sequencing (RNA-seq) data from the Chinese Glioma Genome Atlas (CGGA) database (n = 552). Gene expression was also normalized and calculated using the FPKM method. Rembrandt microarray dataset was downloaded from CGGA database and used as the second validation cohort (n = 123) (Gusev et al., 2018). The 13,609 common genes in 3 independent cohorts are used for the following analysis.

The batch effect of RNA-seq data in the CGGA database, which comprised two independent cohorts, was removed using SVA packages (Brägelmann and Lorenzo Bermejo, 2019). Patients with grade II and III glioma with complete sex, age, IDH1 status, and survival information were used in the following study. The patients with grade II and III glioma in the Rembrandt dataset with complete survival information were used for the follow-up study. For all the included RNA-seq and microarray data, normalization and log2 transformation were performed. The RNA-seq data were adopted in FPKM methods, the microarray data were normalization using RMA methods by oligo package (Carvalho and Irizarry, 2010).

The RNA-seq training data and hallmark gene sets downloaded from the Molecular Signatures database (MSigDB) (Liberzon et al., 2015) were used by a single-sample gene set enrichment analysis (ssGSEA) algorithm (Hänzelmann et al., 2013) as the basis to calculate the quantified score of each cancer hallmark. Univariate Cox analysis was used to calculate the significance of different cancer hallmarks in patients with LGG. With a network type of unsigned and a soft threshold of β = 5 (scale free R ² = 0.8997766), the expression values of protein coding genes in the LGG samples were subjected to weighted gene co-expression network analysis (WGCNA) to construct a scale-free co-expression network, which was used to screen those genes that were most correlated with glycolysis, based on their ssGSEA scores (Langfelder and Horvath, 2008). Subsequently, 19 modules were identified by setting the merged threshold function at 0.25, the green module was identified the genes of significantly related module. The correlation between the gene expression profiles and module eigengenes was measured using module membership (MM), and the correlation between the glycolysis ssGSEA score and individual genes represented the gene significance (GS). We screened 407 extracted candidates from the “glycolysis module” using a cutoff p value of GS of <0.0001 and the univariate Cox regression calculated p value of <0.01. Then, a least absolute shrinkage and selection operator (Lasso) Cox regression algorithm were employed to identify the most significant prognostic markers (Friedman et al., 2010). We then constructed a glycolysis-related risk score (GRS) that included normalized gene expression values weighted by their LASSO Cox coefficients as follow: G R S = ∑ i = 1 n c o e f f i c i e n t ( m R N A i ) ∗ e x p r ( m R N A i )

Based on median value of GRS score, we divided the patients in training cohort into high-risk group and low-risk group, we obtain 1827 differentially expressed genes (p < 0.05, |log₂FC|>1) using the edgeR (Robinson et al., 2010) package of R between two group. Based on the differential genes that have been obtained, Gene ontology and Kyoto encyclopedia of genes and genomes enrichment analyses were launched to probe the potential biological functions and signaling pathways by clusterProfiler (Yu et al., 2012), enrichplot packages (Yu, 2021). Standards of significantly enriched terms were set as follow: q-value < 0.05.

We utilized the msigdb. v7.0. entrez.gmt gene sets from the MsigDB database and clusterProfiler packages to explore the therapeutic resistance and possible cellular pathways. We set following standards for significantly enriched terms: I. NOM p-value < 0.05; II. FDR q-value < 0.05.

The drug sensitivity of the members in GRS score were calculated by online webtool GSCALite (http://bioinfo.life.hust.edu.cn/web/GSCALite/). The results from CTRP database were adopt and displayed in our study.

A decision tree was constructed using recursive partitioning analysis (RPA) in the “rpart” package for risk stratification (Terry and Beth, 2019). For the quantitative prediction of the prognosis of patients with LGG, we established a prognostic nomogram, comprising the glycolysis-based risk model and other clinical parameters, which was used to predict the probability of overall survival (OS) for 1, 3, and 5 years. The predictive ability of the model was demonstrated by plotting a calibration curve, in which a curve close to 45° indicates a good predictive ability.

Data processing and graph construction were carried out using R software (version 4.0.3, http://www.r-project.org). Adobe Illustrator 2020 was used to fine tune the graphics. Survival was evaluated using the Kaplan-Meier method. The differences between the high-risk group and low-risk group in the targeted cohorts were compared using a Log-rank test. To evaluate the risk signature’s predictive performance, receiver operating characteristic curve (ROC) and area under the curve (AUC) at 1, 3, and 5 years were calculated in all cohorts using the “timeROC” packages (Blanche et al., 2013). In the absence of a full-scale gene signature expression profile, each cohort was divided into different clusters according to the optimum k value using K-means consensus clustering in the “ConsensusClusterPlus” packages, The optimal cluster numbers were determined by constructing CDF curves respectively (Wilkerson and Hayes, 2010). In the pooled cohort, to assess the prognostic value, a meta-analysis (I² = 97%, random-effect model) was used.

The information for the patients in each cohort who met the inclusion criteria is listed in Supplementary Table S1. First, among the various hallmarks of cancer, we identified glycolysis as the primary OS-related risk factor in patients with LGG (Figure 1A). Then, WGCNA, univariate Cox regression analysis, and the LASSO algorithm were used to screen candidate genes and construct a robust signature, which could be used to predict patient survival reliably (Figure 1B). Next, the glycolysis signature’s prognostic efficiency was assessed in the training set and two independent validation datasets. In addition, to confirm the signature’s prognostic accuracy, a meta-analysis was preformed, and the usability of the signature in clinical practice was assessed by choosing the response to therapy as an important evaluation indicator (Figure 1C). Finally, we built a decision tree for prognostic precision improvement, and constructed a nomogram based on GRS and other clinicopathological indicators to evaluate the level of risk and the individual patients’ survival probability (Figure 1D).

ssGSEA is an analytical algorithm that uses RNA expression data to score different cancer hallmarks in each sample. The score of 50 hallmark gene sets were calculated and 30 gene sets with a p-value < 0.05 were considered to represent significantly enriched pathways (Figure 2A). Compared with other cancer gene sets, such as those representing apical junction, apoptosis, and angiogenesis, the glycolysis-related gene set had the most powerful influence (highest HR) on survival (Figure 2B). Based on their median ssGSEA score, patients (n = 501) could divide into a low score group and a high score group, which exhibited poorer OS compared with those in the low score group (p < 0.05; Figure 2C).

The top 50 percent of differentially expressed genes (6,805 genes), selected on the basis of an analysis of variance and their glycolysis ssGSEA score in the training set, were chosen for WGCNA. The optimal soft-thresholding power was 5, which ensured that the co-expression network was scale free (Supplementary Figure S1A). Specifically, 20 co-expression modules were identified after merging modules with similarities above 0.75 (Figure 3A). To assess the stability of each module identified in the training cohort, we divided the TCGA data into training and testing cohort to conduct module preservation analysis using the module preservation test (nPermutations = 200) in the WGCNA package; modules with Z-score >10 were considered to be strongly preserved (Supplementary Figure S1B). Among the modules with Z-score >10, the green module showed the highest correlation (r = 0.57, p = 3e-44) and was thus considered to correlate the most with glycolysis (Figure 3B). Using a p value <0.0001 as the filtering criterion for GS, univariate Cox regression analysis was performed on hub genes extracted from the green module. Based on a threshold p value < 0.01 for the univariate Cox regression analysis, 407 candidate genes were submitted to the LASSO Cox regression algorithm, which can identify the strongest prognostic markers for prognosis. Ten-fold cross-validation was used to abrogate over-filtering, by selecting an optimal λ value of 0.0684 (Figure 3C). Fifteen genes (ARL3, EMP3, IGFBP2, PTGFRN, ADAMTS3, ARL9, SEMA4G, RYR3, TNFRSF11B, SSFA2, ABCC3, EMILIN3, IGF2BP2, KLHL9, and RHBDF1) retained their individual non-zero LASSO coefficients (Figure 3D). The LASSO coefficients of each gene in the GRS. Thus, we established the GRS formula, among the members of the signature, KLHL9, ARL9 and SEMA4G are protective factors, the others are risk factors (Figure 3E and Supplementary Table S2).

The median GRS risk score in the training set and two validation sets were used to categorize patients with LGG into low- and high-risk groups. In each cohort, patients in the high-GRS score group had worse OS than those with a low GRS score (Figures 4A–C, p < 0.0001). Time-dependent ROC curves were used to evaluate the reliability of the GRS (Figures 4D–F). The AUC values were 0.78, 0.711, and 0.603 for 1-year, 3-years, and 5-years survival, respectively, in the TCGA training set, implying good reliability of the GRS to monitor survival. In the CGGA validation cohort, the AUC values were 0.731, 0.753, and 0.717; and in the Rembrandt cohort, the AUC values were 0.777, 0.842, and 0.749 for 1-year, 3-years and 5-years survival, respectively. Clustering was performed using a k-means algorithm which divided each cohort into different groups based on the best k value of genes in the GRS signature (Figures 4G–I), basing on the cumulative distribution function (CDF) curves (Supplementary Figure S2A: TCGA, Supplementary Figure S2B: CGGA, Supplementary Figure S2C: Rembrandt). The results showed that the OS differed significantly among the k-means derived groups, which suggested that the classification means based on genes in the signature could be used directly for tumor subtyping (Figure 4J-L). Moreover, a meta-analysis of the prognostic value of the GRS in the pooled cohort, which integrated the TCGA training cohort and the two validation cohorts, showed that among all the 1,176 patients in three cohorts, those with a lower GRS had a better prognosis than those with a high GRS (pooled HR = 1.23, 95% confidence interval (CI): 1.09–1.39; Figure 4M).

To identifying the relationship between the GRS score and the clinicopathological characteristics, we explored the corresponding clinical information of LGG cases in the TCGA training cohort and the CGGA testing cohort (The clinical information in the Rembrandt cohort was incomplete). The GRS score was associated significantly with age, survival status, and IDH1 status. The results showed that the GRS score was high in “Dead (survival status)”, “> 40 (Age)”, “IDH1 wild-type”, non-codel (1p19q status) unmethylated (MGMT status) group, the respective p values were below 0.05 (Figures 5A–F), the same trend was observed in the CGGA cohort (Figure 5G-L).

To explore the association of gene expression with GRS score, gene expression in 501 patients with LGG in the training set was compared between the high and low GRS score groups, to identify differentially expressed genes (DEGs). We obtain 1827 differentially expressed genes (p < 0.05, |log₂FC|>1), Gene ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) were performed on the DEGs. The GO analysis showed that a higher GRS score was associated with extracellular matrix organization and extracellular structure organization (Figure 5M). The KEGG analysis indicated that a higher GRS score was related to ECM−receptor interaction, Complement and coagulation cascades (Figure 5M). Taken together, these results suggested that the proteins encoded by the GRS genes have important functions in tumor microenvironment (TME) remodeling in patients with LGG.

Tumor glycolysis increases therapy resistance; therefore, we determined if the GRS could function as an indicator of therapeutic resistance (Hamilton et al., 2018). In addition to being associated with poor survival, a higher GRS score correlated significantly with resistance to various therapies, such as targeted therapy, radiation therapy, and chemotherapy, according to GSEA in the training cohort (Figure 6A). The online tool GSCALite was used to draw landscape plot (bubble heatmap), which demonstrated the relationships between gene signature members and drug responses (Figure 6B). Genes in the signature correlated significantly with the half-maximal inhibitory concentration (IC50) data in LGG cells. The genes ABCC3, RHBDF1, and PTGFRN conferred drug resistance, which were consistent with the results shown in Figure 3E. Next, clinical outcomes and treatment information from the training cohort were used to verify these predictions. Following surgery, the clinical benefit rate (CR, complete remission; PR, partial remission; SD, stable disease) was remarkably lower (p < 0.0001) in the higher GRS group, both in terms of follow-up treatment outcome and primary treatment outcome (Figure 6C). Furthermore, among patients who had been treated with adjuvant therapies (chemotherapy and radiotherapy), patients with a higher GRS score had worse (p < 0.0001) OS than those with a lower GRS score (Figure 6D).

A decision tree was constructed to increase the stratification of risk using the three parameters available age (≥40 or <40), IDH status (wild and mutant), and the GRS (high and low), the results of which showed that all factors (nsplit = 3, xerror is minimum) remained in the decision tree and four different risk subgroups were identified (Figure 7A). The GRS score played an important role in the model. Kaplan–Meier curves showed that OS differed significantly among the four risk subgroups. Patients in the high GRS score, wild-type IDH1 status, advanced age (age >40) subgroups had the highest risk (Figure 7B). The patients in the training cohort with age, sex, IDH1 status, and GRS scores were used for further study. Univariate Cox regression analysis showed that age, IDH1 status, and GRS score were associated significantly with the prognosis of LGG (Figure 7C, p < 0.05). Multivariate Cox regression analysis identified GRS as an independent prognostic factor combined with other clinicopathological factors (Figure 7D, p < 0.2). These three factors (The GRS score was scaled) were constructed into a nomogram to quantify the risk assessment and survival probability for individual patients with LGG (Figure 7E). The forecast curve of the calibration analysis (black line) of the nomogram for 3-years and 5-years survival probability closely resembled that of the ideal performance (the grey line in Figures 7G,H), indicating that the nomogram was highly accurate. Finally, the nomogram’s reliability was evaluated using time-dependent ROC curves (Figure 7F). The AUC values were 0.875, 0.893, 0.799, and 0.762 for 1-, 3-, 5-, and 10-years survival, respectively, indicating a good potential in clinical practice for monitoring survival. We conduct same analysis in the CGGA cohort, the result was exhibited in Supplementary Figure S3, which also suggest high reliability of our model.