Abstract
The overall survival of patients with lower grade glioma (LGG) that might develop into high-grade malignant glioma shows marked heterogeneity. The currently used clinical evaluation index is not sufficient to predict precise prognostic outcomes accurately. To optimize survival risk stratification and the personalized management of patients with LGG, there is an urgent need to develop an accurate risk prediction model. The TCGA-LGG dataset, downloaded from The Cancer Genome Atlas (TCGA) portal, was used as a training cohort, and the Chinese Glioma Genome Atlas (CGGA) dataset and Rembrandt dataset were used as validation cohorts. The levels of various cancer hallmarks were quantified, which identified glycolysis as the primary overall survival-related risk factor in LGGs. Furthermore, using various bioinformatic and statistical methods, we developed a strong glycolysis-related gene signature to predict prognosis. Gene set enrichment analysis showed that in our model, high-risk glioma correlated with the chemoradiotherapy resistance and poor survival. Moreover, based on established risk model and other clinical features, a decision tree and a nomogram were built, which could serve as useful tools in the diagnosis and treatment of LGGs. This study indicates that the glycolysis-related gene signature could distinguish high-risk and low‐risk patients precisely, and thus can be used as an independent clinical feature.
Introduction
Glioma is the most common intracranial malignant tumor, which accounts for 80% of all intracranial malignancies, with 15,000–17,000 new cases every year in America (). The World Health Organization defined grade II and III gliomas as diffuse lower grade gliomas (LGGs), occupying approximately 30% of intracranial tumors (). LGGs predict a better prognosis than glioblastoma; however, despite their highly heterogeneous natural processes, most LGGs progress to aggressive high-grade secondary gliomas that ultimately cause death (). Despite great advances in treatment for LGGs being achieved, including neurosurgery, radiotherapy, and chemotherapy, the treatment of LGGs remains a challenge (). Recently, our ability to diagnose and prognose LGGs has been enhanced by the discovery of various biomarkers. The classification of CNS tumors was revised by the World Health Organization in 2016, based on morphological and molecular features, which highlighted the importance of molecular testing and signature construction for the diagnosis and prognosis of LGGs ().
Compared with healthy cells, cancer cells have very different metabolic patterns, adjusting metabolism to sustain the biosynthetic demands of tumor proliferation and have a higher resistance to cell death pathways (). Changes in the tumor glucose metabolic environment are often associated with oxygen supply and demand, dysfunction of key enzymes, or mitochondrial dysfunction, frequently resulting in advanced tumor aggressiveness, poor prognosis, and limited efficacy of available therapeutic methods (). Although some studies have shown that tumor glucose metabolism disorders are closely associated with cancer progression and low survival in gliomas, there is no glycolysis-based approach to distinguish high-risk patients (; ).
There has been significant progress in molecular markers; however, they are still not satisfactory and there is still much room for improvement (). The present study aimed to determine whether glycolysis is a primary survival-related risk factor for LGGs, to identify powerful biomarkers, and construct a glycolysis-related gene signature for LGG prognosis. Furthermore, the prognostic efficiency of the model was tested in two additional validation sets. Finally, we established a comprehensive model based on clinicopathological and genetic characteristics to enhance the signature’s accuracy and predictive power, which might be applied to guide the clinical care and patient consultation.
Materials and Methods
Preparation and Processing of Data
We downloaded the transcriptome profiles (HTSeq Fragments Per kilobase of transcript per Million mapped reads [FPKM]) and their associated clinical profiles from TCGAbiolinks (), 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 (). 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) (). 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 (). 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 ().
Screening Candidate Genes and Building the Signature
The RNA-seq training data and hallmark gene sets downloaded from the Molecular Signatures database (MSigDB) () were used by a single-sample gene set enrichment analysis (ssGSEA) algorithm () 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 R2 = 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 (). 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 (). We then constructed a glycolysis-related risk score (GRS) that included normalized gene expression values weighted by their LASSO Cox coefficients as follow:
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis Based on DEGs Between High GRS Score and Low GRS Score
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, |log2FC|>1) using the edgeR () 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 (), enrichplot packages (). Standards of significantly enriched terms were set as follow: q-value < 0.05.
Gene Set Enrichment Analysis
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.
Drug Sensitivity of the Members in GRS Score
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.
Development and Evaluation of Clinical Predictive Models
A decision tree was constructed using recursive partitioning analysis (RPA) in the “rpart” package for risk stratification (). 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.
Bioinformatic and Statistical Analyses
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 (). 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 (). In the pooled cohort, to assess the prognostic value, a meta-analysis (I2 = 97%, random-effect model) was used.
Results
Schematic Diagram of Research Design
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).
FIGURE 1
Glycolysis Is the Primary OS-Related Risk Factor in LGGs
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).
FIGURE 2
Construction of a Prognostic Glycolysis-Related Gene Signature
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).
FIGURE 3
GRS Is a Risk Factor for OS in Each Cohort
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).
FIGURE 4
GRS Score Was Associated With the Clinicopathological Characteristics
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).
FIGURE 5
Pathway Enrichment Analysis and Visualization of Differentially Expressed Gene (DEGs) Between the High and Low GRS Score Groups
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, |log2FC|>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.
GRS Acts as an Indicator of Therapeutic Resistance and Potential Treatments for Patients With High GRS Score Patients
Tumor glycolysis increases therapy resistance; therefore, we determined if the GRS could function as an indicator of therapeutic resistance (). 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).
FIGURE 6
Integrating the Glycolysis Signature With Clinicopathological Features to Improve Risk Stratification and Survival Prediction
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.
FIGURE 7
Discussion
Compared with normal cells, tumor cells have many unique hallmarks, which are of great significance in the precise treatment of tumors. Metabolic change is one of the most intriguing areas of tumor research (). To meet their high metabolism, tumor cells must struggle to get enough nutrients (). More than 90 years ago, Warburg et al. found that tumor cells maintain active glycolysis levels even in the absence of oxygen. This landmark discovery is known as the “Warburg Effect” (). Since then, the study of glucose metabolism in tumors has been in full swing. Positron emission tomography (PET) imaging based on cancer glucose metabolism has been widely used in clinical and plays an important role in the diagnosis and monitoring of tumors (). Therapies targeting glucose metabolism, such as glucose transporter (GLUT) inhibitors, also hold great promise in the treatment of tumors (). The role of glucose metabolism in cancer is becoming increasingly attractive.
As research progressed, the role of glucose metabolism in tumors was found to be related to more than just energy supply. Some intermediates of glucose metabolism, such as lactic acid, may be involved in improving the tumor microenvironment and mediating immune reprogramming (). Clinicians can also be able to choose and tailor treatments based on well-established disease models. Zhang et al. constructed immune infiltrating cells-derived risk signature, which can well describe immune characteristics and predict prognosis in glioma patients (). Therefore, the current research on glucose metabolism is far from enough. In addition, several clinical studies have found that high glucose metabolism is associated with poor prognosis in some tumors. Zhang et al. constructed a prognostic signature of glycolysis-related genes in lung adenocarcinoma (), in which patients in the high-risk group had a worse prognosis. Yu et al. also constructed a prognostic model of glycolysis-related genes in gastric cancer and found that glycolysis was related to prognosis and immune infiltration in gastric cancer patients (). However, for the low-grade glioma, a tumor with a very poor prognosis, no glycolysis-related prognostic model has been established. Therefore, it is necessary to evaluate the prognostic significance and therapeutic guidance value of glycolysis in low-grade glioma.
In this study, using the ssGSEA algorithm and univariate regression model, glycolysis was identified as the primary OS-related risk factor in patients with LGG. Combining algorithm WGCNA, univariate Cox regression and LASSO Cox regression model, we obtained robust prognostic candidates and constructed a glycolysis-related gene signature (GRS). Next, the gene signature’s prognostic value was confirmed by multiple analyses in two other cohorts. The clinical characteristics of patients with LGG were analyzed to determine the correlation between the GRS score and clinicopathological features, which showed that the GRS score was elevated significantly in patients with older age (>40), wild-type IDH1 status, and those who had died. Thereafter, GO, KEGG enrichment analysis based on DEGs between the high and low GRS score groups were carried out to explore what aspects of the high GRS score affected the TME of patients with LGG. In addition, patients with LGG with a high GRS score experienced poorer survival compared with those in the low GRS score group in the adjuvant therapy groups, which might have been caused by therapeutic resistance induced by members of the signature, which suggested that the GRS could also be used as a reliable tool to predict therapeutic resistance in patients with LGG. Finally, a decision tree was constructed to improve the stratification of risk, which was integrated with the clinicopathological characteristics. In the decision tree, the GRS score functioned as the major decisive factor. Meanwhile, the GRS was confirmed as an independent prognostic factor after adjusting for other clinicopathological features using multivariate Cox regression analysis. These results indicated that the GRS is a reliable risk factor for OS in patients with LGG. To increase the predictive ability of OS for LGG in a quantitative manner, a prognostic nomogram containing scale-GRS scores and other clinical features was established to predict the probability of 1-, 3-, 5-years OS. Calibration analysis demonstrated the accurate predictive ability of the nomogram, which was in accord with actual survival. Furthermore, ROC analysis demonstrated that the nomogram model exhibited high accuracy to predict survival in the timeline of follow up.
Some candidates included in the signature, have been reported to be involved in the glucose metabolic process of low-grade glioma cancer. ARL3 cycling between an in active GDP-bound and an active GTP-bound form, is involved in energy metabolic process, which may be associated with glucose metabolic process (; ). Rahman et al. demonstrated that IGFBP2 could induce the increase of glioma invasion and malignancy by activating PTEN and AKT pathways, which is enhanced by HIF-α associated with promoting glycolysis related activities (; ). Li et al. reported that the express of RHBDF1 inhibited the RACK1 induced HIF-α degradation in breast cancer, which aggravated the hypoxic environment of the tumor environment and promoted the activation of the glycolysis pathway (). Based on comprehensive bioinformatic analysis, this study further screened some potentials which have not been reported in glucose metabolic process, which may provide new insight for the further research in glucose metabolic process of glioma.
Meanwhile, glycolysis related signature shows high accuracy in predicting the survival rate in 3 cohorts in our study, furthermore the signature displays strong correlation with clinicopathological features. All the results indicate GRS can serve as a reliable indicator in clinical application. To test the correctness of the model and improve the prediction ability of the model, more independent cohorts should be involved in our study. Go and KEGG analysis reveal fact glycolysis is a tumor hallmark involving gene networks (all the members of GRS) rather than some individual “glycolysis genes”, which is involved in tumor microenvironment remodeling. It gives us new enlightenment for the treatment of glioma which should be considered as a whole from multiple targets. In addition, patients with high GRS score exhibited poorer survival treated with traditional adjuvant therapy methods, which implicated urgent need for new means for treatment. we exhibited the drug sensitivity (IC50) of members of the GRS signature in LGG cells, which is need further biological experiments to validate the effectiveness in vitro and vivo.
At last, this study constructed a robust predictive model which can increase the abilities of risk stratification and survival prediction. However, the clinical information of cases in incorporated datasets were inconsistent. We established a predictive model based on the information available in TCGA training set, not including some important clinical features which can further improve the predictive ability of the model, such as p/19q status, WHO grade. In addition, the further verifications are needed by biological experiments and clinical studies.
Conclusion
The present study has the advantages of incorporating different cohorts from the TCGA and CGGA databases to construct a robust glycolysis signature that can predict patient survival and therapeutic resistance with high accuracy in patients with LGG.
Combining the GRS signature with clinicopathological characteristics allowed us to construct a decision tree and a nomogram, which increased the abilities of risk stratification and survival prediction. The GRS might function as a reliable clinical prediction tool and might aid the future development of therapeutic targets.
Statements
Data availability statement
Publicly available datasets were analyzed in this study. This data can be found here: This study used online resources, which are available from the TCGA database at https://portal.gdc.cancer.gov/ and the CGGA database at http://www.cgga.org.cn/.
Author contributions
WL and CC conceived the study. WL and CL wrote the paper and drew the figures. WL XH, and QY analyzed the data. YY, YT, and BP edited and reviewed the paper. YY, CC, and XH were responsible for funding acquisition. All the authors contributed to the study and approved the submitted version of the article.
Funding
This work was supported by the National Natural Science Foundation of China (grant number 81773354 and 82102971), the Guangdong Basic and Applied Basic Research Foundation, China (grant number 2020A1515111074), and the Guangzhou Science and Technology Plan Project (grant number 202102020034).
Acknowledgments
We would like to thank the Cancer Genome Atlas (TCGA) and the Chinese Glioma Genome Atlas (CGGA) working groups for kindly sharing of data. We also thank R Development Core Team and R community (http://www.r-project.org) for the time and effort they put into to create and optimize R.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.843711/full#supplementary-material
References
1
AbbaszadehZ.ÇeşmeliS.Biray AvcıÇ. (2020). Crucial Players in Glycolysis: Cancer Progress. Gene726, 144158. 10.1016/j.gene.2019.144158
2
AnceyP. B.ContatC.MeylanE. (2018). Glucose Transporters in Cancer - from Tumor Cells to the Tumor Microenvironment. FEBS J.285, 2926–2943. 10.1111/febs.14577
3
BlancheP.DartiguesJ.-F.Jacqmin-GaddaH. (2013). Estimating and Comparing Time-Dependent Areas under Receiver Operating Characteristic Curves for Censored Event Times with Competing Risks. Statist. Med.32, 5381–5397. 10.1002/sim.5958
4
BrägelmannJ.Lorenzo BermejoJ. (2019). A Comparative Analysis of Cell-Type Adjustment Methods for Epigenome-Wide Association Studies Based on Simulated and Real Data Sets. Brief Bioinform20, 2055–2065. 10.1093/bib/bby068
5
BratD. J.BratD. J.VerhaakR. G.AldapeK. D.YungW. K.SalamaS. R.et al (2015). Comprehensive, Integrative Genomic Analysis of Diffuse Lower-Grade Gliomas. N. Engl. J. Med.372, 2481–2498. 10.1056/NEJMoa1402121
6
CarvalhoB. S.IrizarryR. A. (2010). A Framework for Oligonucleotide Microarray Preprocessing. Bioinformatics26, 2363–2367. 10.1093/bioinformatics/btq431
7
ColapricoA.SilvaT. C.OlsenC.GarofanoL.CavaC.GaroliniD.et al (2016). TCGAbiolinks: An R/Bioconductor Package for Integrative Analysis of TCGA Data. Nucleic Acids Res.44, e71. 10.1093/nar/gkv1507
8
DasS. K.BhutiaS. K.AzabB.KegelmanT. P.PeachyL.SanthekadurP. K.et al (2013). MDA-9/Syntenin and IGFBP-2 Promote Angiogenesis in Human Melanoma. Cancer Res.73, 844–854. 10.1158/0008-5472.can-12-1681
9
FriedmanJ.HastieT.TibshiraniR. (2010). Regularization Paths for Generalized Linear Models via Coordinate Descent. J. Stat. Softw.33, 1–22. 10.18637/jss.v033.i01
10
GabrielyG.WheelerM. A.TakenakaM. C.QuintanaF. J. (2017). Role of AHR and HIF-1α in Glioblastoma Metabolism. Trends Endocrinol. Metab.28, 428–436. 10.1016/j.tem.2017.02.009
11
GusevY.BhuvaneshwarK.SongL.ZenklusenJ.-C.FineH.MadhavanS. (2018). The REMBRANDT Study, a Large Collection of Genomic Data from Brain Cancer Patients. Sci. Data5, 180158. 10.1038/sdata.2018.158
12
HamiltonK. E.RekmanJ. F.GunninkL. K.BusscherB. M.ScottJ. L.TidballA. M.et al (2018). Quercetin Inhibits Glucose Transport by Binding to an Exofacial Site on GLUT1. Biochimie151, 107–114. 10.1016/j.biochi.2018.05.012
13
HänzelmannS.CasteloR.GuinneyJ. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. Bmc Bioinformatics14, 7. 10.1186/1471-2105-14-7
14
LangfelderP.HorvathS. (2008). WGCNA: An R Package for Weighted Correlation Network Analysis. Bmc Bioinformatics9, 559. 10.1186/1471-2105-9-559
15
LibertiM. V.LocasaleJ. W. (2016). The Warburg Effect: How Does it Benefit Cancer Cells?Trends Biochem. Sci.41, 211–218. 10.1016/j.tibs.2015.12.001
16
LiberzonA.BirgerC.ThorvaldsdóttirH.GhandiM.MesirovJ. P.TamayoP. (2015). The Molecular Signatures Database (MSigDB) Hallmark Gene Set Collection. Cel Syst.1, 417–425. 10.1016/j.cels.2015.12.004
17
LouisD. N.PerryA.ReifenbergerG.von DeimlingA.Figarella-BrangerD.CaveneeW. K.et al (2016). The 2016 World Health Organization Classification of Tumors of the Central Nervous System: A Summary. Acta Neuropathol.131, 803–820. 10.1007/s00401-016-1545-1
18
MayakondaA.LinD.-C.AssenovY.PlassC.KoefflerH. P. (2018). Maftools: Efficient and Comprehensive Analysis of Somatic Variants in Cancer. Genome Res.28, 1747–1756. 10.1101/gr.239244.118
19
OstromQ. T.GittlemanH.TruittG.BosciaA.KruchkoC.Barnholtz-SloanJ. S. (2018). CBTRUS Statistical Report: Primary Brain and Other Central Nervous System Tumors Diagnosed in the United States in 2011-2015. Neuro Oncol.20, iv1–iv86. 10.1093/neuonc/noy131
20
OuX.LvW. (2020). Metabolic Changes and Interaction of Tumor Cell, Myeloid-Derived Suppressor Cell and T Cell in Hypoxic Microenvironment. Future Oncol.16, 383–393. 10.2217/fon-2019-0692
21
PavlovaN. N.ThompsonC. B. (2016). The Emerging Hallmarks of Cancer Metabolism. Cel Metab.23, 27–47. 10.1016/j.cmet.2015.12.006
22
RahmanM. S.ThomasP. (2011). Characterization of Three IGFBP mRNAs in Atlantic Croaker and Their Regulation during Hypoxic Stress: Potential Mechanisms of Their Upregulation by Hypoxia. Am. J. Physiology-Endocrinology Metab.301, E637–E648. 10.1152/ajpendo.00168.2011
23
ReckzehE. S.WaldmannH. (2020). Development of Glucose Transporter (GLUT) Inhibitors. Eur. J. Org. Chem.2020, 2321–2329. 10.1002/ejoc.201901353
24
RobinsonM. D.McCarthyD. J.SmythG. K. (2010). edgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics26, 139–140. 10.1093/bioinformatics/btp616
25
StricklandM.StollE. A. (2017). Metabolic Reprogramming in Glioma. Front. Cel Dev. Biol.5, 43. 10.3389/fcell.2017.00043
26
TerryT.BethA. (2019). Rpart: Recursive Partitioning and Regression Trees. R Package Version 4.1-15. Available at: https://CRAN.R-project.org/package=rpart.
27
VaupelP.SchmidbergerH.MayerA. (2019). The Warburg Effect: Essential Part of Metabolic Reprogramming and central Contributor to Cancer Progression. Int. J. Radiat. Biol.95, 912–919. 10.1080/09553002.2019.1589653
28
VeltelS.KravchenkoA.IsmailS.WittinghoferA. (2008). Specificity of Arl2/Arl3 Signaling Is Mediated by a Ternary Arl3-Effector-GAP Complex. Febs Lett.582, 2501–2507. 10.1016/j.febslet.2008.05.053
29
WalkerS. M.LimI.LindenbergL.MenaE.ChoykeP. L.TurkbeyB. (2020). Positron Emission Tomography (PET) Radiotracers for Prostate Cancer Imaging. Abdom. Radiol.45, 2165–2175. 10.1007/s00261-020-02427-4
30
WilkersonM. D.HayesD. N. (2010). ConsensusClusterPlus: a Class Discovery Tool with Confidence Assessments and Item Tracking. Bioinformatics26, 1572–1573. 10.1093/bioinformatics/btq170
31
YoulandR. S.SchomasD. A.BrownP. D.ParneyI. F.LaackN. N. I. (2017). Patterns of Care and Treatment Outcomes in Older Adults with Low Grade Glioma: a 50-year Experience. J. Neurooncol.133, 339–346. 10.1007/s11060-017-2439-3
32
YuG. (2021). Enrichplot: Visualization of Functional Enrichment Result. Available at: https://yulab-smu.top/biomedical-knowledge-mining-book/.
33
YuG.WangL.-G.HanY.HeQ.-Y. (2012). clusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS16, 284–287. 10.1089/omi.2011.0118
34
YuS.HuC.CaiL.DuX.LinF.YuQ.et al (2020). Seven-Gene Signature Based on Glycolysis Is Closely Related to the Prognosis and Tumor Immune Infiltration of Patients with Gastric Cancer. Front. Oncol.10, 1778. 10.3389/fonc.2020.01778
35
ZhangC.ChengW.RenX.WangZ.LiuX.LiG.et al (2017). Tumor Purity as an Underlying Key Factor in Glioma. Clin. Cancer Res.23, 6279–6291. 10.1158/1078-0432.ccr-16-2598
36
ZhangH.LiX.LiY.ChenB.ZongZ.ShenL. (2020). An Immune-Related Signature for Predicting the Prognosis of Lower-Grade Gliomas. Front. Immunol.11, 603341. 10.3389/fimmu.2020.603341
37
ZhangL.ZhangZ.YuZ. (2019). Identification of a Novel Glycolysis-Related Gene Signature for Predicting Metastasis and Survival in Patients with Lung Adenocarcinoma. J. Transl Med.17, 423. 10.1186/s12967-019-02173-2
38
ZhangN.ZhangH.WangZ.DaiZ.ZhangX.ChengQ.et al (2021). Immune Infiltrating Cells-Derived Risk Signature Based on Large-Scale Analysis Defines Immune Landscape and Predicts Immunotherapy Responses in Glioma Tumor Microenvironment. Front. Immunol.12, 691811. 10.3389/fimmu.2021.691811
39
ZhouC.CunninghamL.MarcusA. I.LiY.KahnR. A. (2006). Arl2 and Arl3 Regulate Different Microtubule-Dependent Processes. Mol. Biol. Cel17, 2476–2487. 10.1091/mbc.e05-10-0929
40
ZhouZ.LiuF.ZhangZ.-S.ShuF.ZhengY.FuL.et al (2014). Human Rhomboid Family-1 Suppresses Oxygen-Independent Degradation of Hypoxia-Inducible Factor-1α in Breast Cancer. Cancer Res.74, 2719–2730. 10.1158/0008-5472.can-13-1027
Summary
Keywords
low-grade glioma, glycolysis, signature, prognosis, risk model
Citation
Liu W, Liu C, Chen C, Huang X, Yi Q, Tian Y, Peng B and Yuan Y (2022) Construction and Verification of a Glycolysis-Associated Gene Signature for the Prediction of Overall Survival in Low Grade Glioma. Front. Genet. 13:843711. doi: 10.3389/fgene.2022.843711
Received
26 December 2021
Accepted
07 March 2022
Published
23 March 2022
Volume
13 - 2022
Edited by
Li Xing, University of Saskatchewan, Canada
Reviewed by
Longbo Zhang, Yale University, United States
Mingjie Wang, Shanghai Jiao Tong University, China
Updates
Copyright
© 2022 Liu, Liu, Chen, Huang, Yi, Tian, Peng and Yuan.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yunhong Tian, tianyunhong2020@gzhmu.edu.cn; Biao Peng, 2747617668@qq.com; Yawei Yuan, yuanyawei@gzhmu.edu.cn
† These authors have contributed equally to this work
This article was submitted to Computational Genomics, a section of the journal Frontiers in Genetics
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.