
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)13405-6
10.1016/j.heliyon.2024.e37374
e37374
Research Article
Identification and validation of a novel robust glioblastoma prognosis model based on bioinformatics
Zhang Le ac1
Yan Xiaoling d1
Wang Yahong e1
Wang Qin Wangqin-0209@163.com
ac⁎⁎⁎
Yan Hua yanhua20042007@sina.com
af⁎⁎
Yan Yan yanyandoctor_1982@tju.edu.cn
abc⁎
a Clinical College of Neurology, Neurosurgery and Neurorehabilitation, Tianjin Medical University, Tianjin, 300070, China
b Medical Engineering and Translational Medicine Research Institute, Tianjin University of Medicine, Tianjin, 300072, China
c Department of Clinical Laboratory, Tianjin Huanhu Hospital, Tianjin, 300350, China
d Department of Pathology, Tianjin Huanhu Hospital, Tianjin, 300350, China
e Nero-oncology Center, Tianjin Huanhu Hospital, Tianjin, 300350, China
f Tianjin Key Laboratory of Cerebral Vascular and Neurodegenerative Diseases, Tianjin Neurosurgical Institute, Tianjin Huanhu Hospital, Tianjin, 300350, China
⁎ Corresponding author. Clinical College of Neurology, Neurosurgery and Neurorehabilitation, Tianjin Medical University, Tianjin, 300070, China. yanyandoctor_1982@tju.edu.cn
⁎⁎ Corresponding author. Clinical College of Neurology, Neurosurgery and Neurorehabilitation, Tianjin Medical University, Tianjin, 300070, China. yanhua20042007@sina.com
⁎⁎⁎ Corresponding author. Clinical College of Neurology, Neurosurgery and Neurorehabilitation, Tianjin Medical University, Tianjin, 300070, China. Wangqin-0209@163.com
1 These authors contributed equally to this work and share first authorship.

03 9 2024
30 9 2024
03 9 2024
10 18 e3737428 5 2024
6 8 2024
2 9 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
Background

Glioblastoma (GBM) is a very common primary malignant tumor of the central nervous system (CNS). Aging, macrophage, autophagy, and methylation related genes are hypothesized to be crucial to its pathogenesis. In this study, we aimed to explore the role of these genes in the prognosis of GBM.

Methods

The RNA sequence (RNA-seq) and clinical information were downloaded from The Cancer Genome Atlas database (TCGA) and the Chinese Glioma Genome Atlas database (CGGA). We performed univariate and least absolute shrinkage and selection operator (LASSO) multivariate Cox regression analysis to identify risk signatures related to overall survival (OS). We further developed a nomogram to predict individual outcomes. In addition, the immune microenvironment was analyzed by CIBERSORT.

Results

256 differentially expressed genes (DEGs) were obtained based on aging, macrophage, autophagy, and methylation related genes between GBM samples and normal tissues in TCGA-GBM cohort. We identified five optimal risk signatures with prognostic values in TCGA-GBM cohort and established a prognostic risk score model. The validity of the model was verified in the CGGA cohort and Huanhu cohort. Finally, we constructed a nomogram for clinical application by combining age, O(6)-methylguanine-DNA methyltransferase (MGMT) promoter methylation status, and risk score. Activated NK cells and resting mast cells were highly expressed and memory B cells, plasma cells, resting NK cells, M1 macrophages, and neutrophils exhibited low expression in the high-risk score group. GBM patients with a low-risk score had a higher Tumor Immune Dysfunction and Exclusion (TIDE) score. The risk score of hot tumors was higher than that of the cold tumors. Additionally, 29 genes involved in glucose and lipid metabolism were highly expressed with a high-risk score. 31 metabolism-related pathways were significantly different between high-risk and low-risk groups.

Conclusions

We constructed and validated a novel prognostic model for GBM. Aging, macrophage, autophagy, and methylation related genes may serve as prognostic and therapeutic biomarkers. The model developed may assist in guiding treatment for GBM patients. Our research had great significance in accurately predicting the prognosis of GBM and may offer reference for immunotherapy decision for GBM patients.

Highlights

• A novel prognostic model for GBM was constructed and validated.

• Aging, macrophage, autophagy, and methylation related genes were integrated to construct a comprehensive prognostic model.

• Immune infiltration and metabolic reprogramming were used to validate the prognostic model.

Keywords

Glioblastoma
Risk score
Nomogram
OS
Immune infiltration
Hot/cold tumor
Metabolic reprogramming
Abbreviations

GBM glioblastoma

CNS central nervous system

LASSO least absolute shrinkage and selection operator

TME tumor microenvironment

RNA-seq RNA sequence

ICIs immune checkpoint inhibitors

MAP1LC3A microtubule-associated protein 1 light chain 3 alpha

MPO myeloperoxidase

GATA4 GATA binding protein 4

S100A4 S100 calcium-binding protein A4

SPP1 secreted phosphoprotein 1

DEGs differentially expressed genes

PPI protein-protein interaction

AUC area under the curve

HRs hazard ratios

TIDE tumor immune dysfunction and exclusion

PCA principle component analysis

SSE squared errors

ECM extracellular matrix

ARGs autophagy-related genes

ROC receiver operating characteristic

CI Confidence Intervals

TCGA The Cancer Genome Atlas database

CGGA Chinese Glioma Genome Atlas database

IDH isocitrate dehydrogenase (IDH)

ATRX alpha-thalassemia/mental retardation syndrome X-linked

TERT telomerase reverse transcriptase

MGMT O(6)-methylguanine-DNA methyltransferase

KEGG Kyoto Encyclopedia of Genes and Genomes
==== Body
pmc1 Introduction

Glioma is a common primary tumor of the CNS that originates from neuroglial progenitor cells. It accounts for 81 % of CNS malignancies [1]. According to the fifth edition of the WHO Classification of Tumors in the CNS, glioma can be classified into four grades: I, II, III, and IV. The higher the grade the tumor is, the worse the prognosis. GBM, which belongs to WHO grade IV tumor, is highly malignant and indicative of poor prognosis [2].The survival time of GBM patients is approximately 15 months [3]. Given the high recurrence rate, surgical removal does not improve survival. Studies have shown that 77,000 new GBM cases were diagnosed in the United States and Europe, but the five-year survival rate was <3 % [4]. In recent years, aging, macrophage, autophagy and methylation are reported to associate with the development of glioma [[5], [6], [7], [8]]. Moreover, these four signatures are closely related to the prognosis of GBM. The prognostic model of GBM based on aging [9], autophagy [10], macrophage [11], methylation [12] can predict the prognosis of GBM patients, respectively. Thus, we hypothesized that the combination of these signatures could better predict the prognosis of GBM. We attempted to construct an OS prognosis model based on these signatures. Although many tumor prognosis models related to these gene signatures were built based on the development of public databases, we integrated these four risk factors for the first time to construct a more comprehensive prognosis model of GBM.

Today, RNA-seq is widely performed to observe the differences of cells at the gene level. Several studies have investigated molecular biomarkers for GBM using RNA-seq technology. Additionally, numerous prognostic models have been constructed for glioma analysis using large-scale genome-sequencing technologies. Jian Zhou et al. developed and verified a glioma prognostic model based on the immune microenvironment [13]. Jia Hao Bao et al. found a novel cuproptosis-related gene signature in lower-grade glioma and verified the expression of hub genes by quantitative real-time PCR and western blotting [14]. Houshi Xu et al. reported that ITGB2 could serve as a prognostic maker and a predictive indicator for immunological therapy in glioma [15]. Additional studies have indicated that a certain single signature does not completely represent the entire feature set of a certain tumor. However, the combination of multiple genes may serve as excellent molecular biomarkers [16,17]. Here, we investigated a novel prognostic model in the pathogenesis of GBM by combining multiple signatures.

In recent years, the tumor microenvironment (TME) has been identified as an important regulator in the promotion of cancer progression in both primary and metastatic brain malignancies. It is the surrounding microenvironment of tumor cells, including blood vessels, immune cells, fibroblasts, bone marrow-derived inflammatory cells, various signaling molecules, and extracellular matrix (ECM). Several studies have shown that a tumor immune microenvironment that contains several types of immune cells can affect the efficacy of immune checkpoint inhibitors (ICIs) [18,19]. Thus, we explored the immunological mechanism of the TME in GBM with the aim of improving patient outcomes.

Over the years, more and more studies showed that metabolic reprogramming is the metabolic changes that cells undergo in response to various stressors, and it is one of the characteristics of tumors. This mechanism not only allows cells to resist external stresses, but also imparts new functions on them. Metabolic reprogramming is a common phenomenon in various tumors, involving metabolic pathways such as glucose metabolism, lipid metabolism, amino acid metabolism, and other pathways, which are closely related to the occurrence and development of tumors [20].

In this study, we acquired the RNA-seq transcriptome data and clinical information of GBM patients from TCGA and CGGA. First, we combined aging, macrophage, autophagy and methylation signatures and identified 256 overlapping genes. Then, the candidate hub genes were selected by univariate and LASSO multivariate Cox regression analysis. Five risk genes were determined: Microtubule-associated protein 1 light chain 3 alpha (MAP1LC3A), myeloperoxidase (MPO), GATA binding protein 4 (GATA4), S100 calcium-binding protein A4 (S100A4) and secreted phosphoprotein 1 (SPP1), respectively. Finally, a novel prognostic nomogram including the risk signatures, age, and MGMT promoter methylation status was constructed. CGGA data and the Huanhu cohort were used as external validation datasets to verify the predictive performance of this signature. Moreover, we used CIBERSORT and the LM22 feature matrix to estimate the differences in immune infiltration in high-risk and low-risk score GBM patients. Although several studies have reported on various aspects of GBM, the prognostic biomarkers have not been elucidated clearly to the best of our knowledge. In view of this, we fully considered the various factors related to the prognosis of GBM and combined them well in order to construct a new and more effective prognosis model of GBM.

2 Materials and methods

2.1 Data collection and processing

We downloaded the mRNA expression profile data and clinical information of 161 GBM samples and 5 normal samples from TCGA (https://tcga-data.nci.nih.gov/tcga/). Among these GBM samples, 152 patients with complete survival information were used for the following analyses. In addition, the mRNA expression profile data and clinical information of 301 GBM samples were obtained from CGGA (http://www.cgga.org.cn/). Among these GBM samples, 285 patients with complete survival information were selected as the validation cohort. All the data was obtained using an Agilent-014850 Whole Human Genome Microarray 4 × 44K G4112F platform and used in accordance with TCGA and CGGA publication requirements. The details of the samples are presented in Tables S1 and S2. The aging, macrophage, autophagy, and methylation related genes were obtained from the Gene Set Enrichment Analysis (GSEA, https://www.gsea-msigdb.org/gsea) database. If there were duplicate data for the same patient, the average RNA expression data was selected for use.

2.2 Patients in the Huanhu cohort and sample collection

Sample collection and patient information were approved by the Ethics Committee of Tianjin Huanhu Hospital (Ethical Approval No.2024-047, 2024/02/28) and written informed consent was provided by all the patients who participated in this study. Between 2020 and 2023, 50 patients received surgery and were diagnosed as GBM by at least two experienced pathologists. All samples were stained with hematoxylin and eosin (H&E) to observe the morphology of tumor cells. Then, samples were fixed using formalin and embedded in paraffin. The protein expression levels of the risk genes were examined by immunohistochemistry (IHC). All the patients were successfully followed up by phone calls. Table S3 shows the detailed information of all patients.

2.3 IHC

Paraffin-embedded surgically removed GBM samples from patients who received surgery at Huanhu Hospital were sectioned in 3 μm intervals and mounted on glass slides. Then, the slides were placed in an oven at 60 °C for 1 h to deparaffinize the samples and then rehydrated. Antigen retrieval was performed in a pressure cooker in 10 mmol/L Tris-citrate buffer (pH: 6.0). To reduce endogenous peroxidase activity, the slides were incubated with 3 % hydrogen peroxide at room temperature for 10 min. PBS was used to wash for 3 times, after which samples were incubated with goat serum at room temperature for 30 min. Slides were incubated with primary antibodies overnight at 4 °C. All slides were incubated with the secondary antibody (ab288151, 1:200; Abcam) for 30 min at room temperature after washing with PBS 3 times. All slides were washed with PBS and then incubated with 3.3′-diaminobenzidine (DAB) solution for 5 min. Subsequently, slides were washed with water before counterstaining with hematoxylin for 5 min at room temperature. At least two pathologists evaluated the results of staining images. The primary antibodies used were: Anti-MAP1LC3A antibody (ab52768, 1:100; Abcam), anti-MPO antibody (ab208670, 1:200; Abcam), anti-GATA antibody (ab199428, 1:200; Abcam), anti-S100A4 antibody (ab197896, 1:200; Abcam) and anti-SPP1 antibody (ab214050, 1:200; Abcam). Images were captured using an upright fluorescence microscope (Olympus BX63, Japan). The IHC images were taken at 400 × magnification.

2.4 Analysis of the DEGs

A total of 152 GBM samples and 5 normal controls in TCGA-GBM cohort were compared to identify DEGs. DEGs were analyzed using the empirical Bayesian approach using the R package limma (version 4.1.0) [18]. The screening criteria were |log2-fold change (FC)|>2.0 and FDR<0.05 [21]. Benjamini and Hochberg developed the adjusted p-value method to control the FDR [22]. All the DEGs in TCGA-GBM cohort were visualized using a volcano plot.

2.5 Enrichment analysis

GSEA (http://software.broadinstitute.org/gsea/index.jsp) was performed using the R package Cluster Profiler [23]. GSEA was used for identifying the immunological signaling pathways and corresponding immune genes involved in high-risk and low-risk groups in TCGA-GBM cohort [24]. GSVA was used for identifying the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways between high-risk and low-risk groups in TCGA-GBM cohort. The enriched gene sets were considered significantly enriched when the P value < 0.05.

2.6 Protein-protein interaction (PPI) network

The String database (https://string-db.org/, version 11.0) is one of the most popular databases for studying protein interactions [25]. A total of 2031 organisms including 9.6 million proteins were installed in the database. A total of 184 million interactions were connected in these proteins. It was used to analyze and predict the functional association and the interaction of proteins. These protein interactions included both direct physical effects and indirect functional correlations. We first used Cytoscape (version 3.7.2) [26] to construct a PPI network for preliminary screening. Then, selected genes were analyzed using the cytohubba plugin in Cytoscape based on the Maximum neighborhood component (MNC) algorithm for further screening. The top 100 genes with the highest connectivity were screened. Finally, a PPI network of the DEGs was constructed.

2.7 Construction and validation of a GBM prognostic risk score model

Based on the clinical survival information provided by the TCGA database, univariate Cox regression analysis was performed to screen genes significantly related to the prognosis of GBM patients in TCGA training cohort. A total of 12 genes were obtained with P < 0.001 as the threshold. Then, to prevent overfitting of the models, LASSO Cox regression analysis was performed to further screen the optimal candidate prognostic related genes using the R package glmnet (http://CRAN.R-project.org/package=glmnet) [27]. In the LASSO Cox analysis, 1000 times subsamples were repeated in TCGA cohort and the genes that we chose were repeated ＞900 times. And five aging, autophagy, macrophage and methylation related genes (MAP1LC3A, MPO, GATA4, S100A4, SPP1) were obtained. Based on the LASSO prognostic coefficient of the optimal genes and the expression levels of the target genes in the database, the risk score of each sample was calculated using the following formula:Riskscore=∑i=1nCoefi*Xi

where Coefi is the risk coefficient calculated using the LASSO Cox model and Xi is the gene expression value.

The prognosis of the GBM patients was predicted using the risk score in TCGA-GBM cohort. All patients were classified into low-risk and high-risk groups according to the median value of the risk score. K-M survival curve using the survival package was performed to evaluate the OS in the two groups. The survival differences of the two groups were evaluated using a two-sided log-rank test.

The performance of the risk score algorithm was validated in the CGGA and Huanhu cohort to verify the predictive accuracy of the model constructed in TCGA training cohort. The K-M survival analysis was employed to verity the validity of the risk model in CGGA and Huanhu validarion cohort.

2.8 Construction and evaluation of the nomogram

Multivariate Cox regression analysis was used to investigate whether the risk score could be regarded as an independent predictor for survival. We included age, sex, isocitrate dehydrogenase (IDH) mutation, alpha-thalassemia/mental retardation syndrome X-linked (ATRX) mutation, telomerase reverse transcriptase (TERT) promoter mutation, and MGMT promoter methylation for independent prognostic analysis. We generated a systematic nomogram based on the clinical parameters and risk score by using the R package rms (https://cran.r-project. org/web/packages/rms/) for predicting the 1-, 3-, and 5-year survival probability. Thus, we combined bioinformatics analysis with the clinical characteristics. The discriminative capability of the nomogram was evaluated using the calibration plots [28]. It was used to assess the difference between the predicted outcomes and the actual outcomes. When the calibration plot is close to the 45° line, it represents best prediction. Furthermore, receiver operating characteristic (ROC) curves were generated to determine the accuracy of the prognostic models. The value of area under the curve (AUC) provides a quantitative measure of the predictive power. The Confidence Intervals (CI) around the AUC values was calculated to provide a measure of the uncertainty associated with the prediction. Whether the risk score model can be used as an independent and reliable prognostic indicator was finally determined by calculating the AUC value and CI. Hazard ratios (HRs) and 95 % CIs were provided where appropriate. P < 0.05 was considered statistically significant.

2.9 Estimation of the immune cell infiltration levels and the immunotherapeutic response between high-risk and low-risk groups

CIBERSORT algorithm was utilized to analyze the immune cell types in TME in detail. CIBERSORT (http://cibersort.stanford.edu/) was used for the deconvolution of the expression matrix of human immune cell subtypes based on the linear support vector regression principle [29]. A total of 22 immune cells and 547 barcode gene signatures were set as the default provision. The sum of all estimated proportions of immune cell types in each sample was equal to one for each GBM sample. The Tumor Immune Dysfunction and Exclusion (TIDE) (http://tide.dfci.harvard.edu/) algorithm was a computational method to predict immunotherapy efficacy. It was used to reveal tumor immune evasion mechanisms [30]. Signatures of T cell dysfunction and exclusion predicted cancer immunotherapy response and determined whether the immunotherapy was effective for individuals. The lower the TIDE score was, the better the immunotherapy response was. The ConsensusClusterPlus R package was used to screen hot/cold tumors [31]. Metabolic reprogramming genes were downloaded from Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/) [32].

2.10 Statistical analysis

Statistical analysis was performed using R software version 3.5.2 and GraphPad Prism version 8 (GraphPad Software, Inc.) For grouping similar mRNA expression profiles among the GBM samples, Principal Component Analysis (PCA) and K-means cluster plots were performed. Datasets were divided into K clusters using the K-means clustering technique. The empirical Bayesian approach of limma was applied to identify the DEGs. The two-sided log-rank test was applied to compare the groups for the K-M survival curves. Moreover, we used a χ2 test to compare sex ratios, IDH mutation ratios, ATRX mutation ratios, TERT promoter mutation ratios, and MGMT promoter methylation ratios between the GBM patients and healthy controls. A Wilcoxon rank sum test was used to compare the differences in immune cell infiltration, TIDE score, hot/cold tumors, as well as metabolic reprogramming in the two risk groups. For all analyses, P < 0.05 was considered statistically significant.

3 Results

3.1 Identification of differentially expressed genes between GBM and normal tissues

It is well known that DEGs are often associated with the occurrence and development of disease, including tumors. Therefore, we downloaded the mRNA expression profile data of 152 GBM patients and the corresponding clinical information from TCGA. All patients had complete survival information. First, we analyzed the DEGs between GBM and normal tissues and obtained a total of 3750 up-regulated and 4528 down-regulated genes. All the DEGs are shown in the volcano plot in Fig. 1A and listed in Table S4). Then, to investigate the functions of aging, macrophage, autophagy and methylation in GBM, we obtained the related genes from the GSEA database. The intersection between the above genes and DEGs was obtained, and a total of 256 overlapping genes were obtained (Table S5) for further study. Finally, to analyze the difference of prognosis after clustering and further explore whether aging, macrophage, autophagy and methylation related genes have potential correlation with the prognosis of GBM, we performed K-means cluster analysis of the cancer samples based on the expression levels of these 256 genes. K-means is an unsupervised clustering technique. All data could be separated into K clusters. We selected the K value according to the sum of the squared errors (SSE). Minimized SSE is widely used as a simple and direct method to obtain the value of K. Here, we used the K-means clustering algorithm to classify the genes into 3 groups (Fig. 1B and C). Meanwhile, another method by Gap statistical method was used to calculate the value of K. The result showed that when K = 3, it is the most significant, which is consistent with our previous result (Fig. 1D). Furthermore, PCA analysis was performed to verify the accuracy of the divided clusters (Fig. 1E). In summary, we screened out DEGs between GBM and normal tissues and samples were divided into 3 clusters for the following study.Fig. 1 Identification and classification of DEGs in GBM patients. (A) Volcano plot showing the DEGs of the 152 GBM and 5 normal samples from TCGA. The horizontal axis represents the log2(fold change [FC]) and the vertical axis represents the -log10 of P value. The red dots indicate the upregulated genes, and the blue dots indicate the downregulated genes. The black dots had no statistical significance. (B and C): The 256 overlapping genes were classified into 3 clusters. (D): Another method (Gap statistic method) further confirmed the 256 genes were appropriate to classify into 3 clusters. (E) PCA analysis was performed to separate the overlapping genes into 3 clusters.

Fig. 1

3.2 Comparison of the characteristics of patients in the 3 clusters

Since the samples were divided into 3 clusters, we wondered if there were any differences of signatures among them. Through statistical analysis, we found that there was no significant difference in the age of patients among the 3 clusters (Fig. 2A). However, sex, IDH mutation, ATRX mutation, TERT promoter mutation, and MGMT promoter methylation did differ significantly among the 3 clusters (Fig. 2B–F). In addition, Kaplan-Meier OS analysis showed that there were significant differences in the overall survival of the samples, and cluster 3 had the worst prognosis (Fig. 2G). These results indicated that GBM samples could be classified into different prognostic types based on these 256 genes. Aging, macrophage, autophagy, and methylation related genes were potentially correlated with the prognosis of GBM.Fig. 2 Comparison of the characteristics of the 3 clusters for prognosis. (A-F) Age, sex, IDH mutation, ATRX mutation, TERT promoter mutation, and MGMT promoter methylation were compared across the 3 clusters. (G) The overall survival of the 3 clusters was compared using K-M curves; cluster 3 had the worst prognosis. **P＜0.01, ***P＜0.001.

Fig. 2

3.3 Construction and validation of the prognostic risk score signatures

To explore the potential interactions among these 256 genes, a PPI network was constructed (Fig. S1A). An algorithm based on MNC was employed to select the top 100 genes with the highest connectivity using the cytohubba plugin in Cytoscape (Fig. S1B, Table S6). Next, we performed univariate Cox regression analysis to screen the top 100 genes related to OS. The HR of each gene was calculated; with P < 0.05 as the threshold, 12 genes were finally obtained (Fig. 3A). Among these, a gene called DUSP6 was not shown on the PPI map because it was not related to any other gene. These 12 genes had an HR＞1 and were associated with a poor prognosis in GBM patients (Fig. 3B). Subsequently, to ensure the stability of the model, LASSO Cox regression analysis was performed for further screening. We determined the optimal five genes with the minimum lambda values (Fig. 3C and D). The five prognostic-related genes were MAP1LC3A, MPO, GATA4, S100A4, and SPP1, respectively. Finally, the normalized expression of each gene was weighted based on the regression coefficient of the LASSO Cox regression analysis. We constructed a risk score model to predict the survival probability of GBM patients (risk score = 0.052998210, normalized expression level of MAP1LC3A + 0.008082386 × normalized expression level of MPO + 0.004006090 × normalized expression level of GATA4 + 0.003168030 × normalized expression level of S100A4 + 0.014389293 × normalized expression level of SPP1). We calculated the risk score of each patient according to the formula. The median risk score was used as a cut-off value to stratify patients of TCGA training cohort and CGGA validation cohort into two groups: the high-risk group and the low-risk group. The K-M curves for OS showing both in TCGA training cohort (Fig. 3E) and CGGA validation cohort (Fig. 3F) were significantly different. The TCGA-GBM cohort included 76 high-risk score patients and 76 low-risk score patients. The high-risk group showed a 2.05-fold higher risk (95 % confidence interval (CI): 1.42–2.96, P＜0.0001). The CGGA cohort included 143 high-risk score patients and 142 low-risk score patients. The high-risk group showed a 1.78-fold higher risk (95 % CI: 1.33–2.38, P＜0.0001).Fig. 3 Construction and validation of a GBM prognostic model. (A) The PPI network of the screened genes. (B) Univariate Cox regression analysis of the top 12 risk signatures among the 100 genes. (C and D) LASSO Cox regression model. (E and F) K-M survival curve in TCGA training cohort (F) and the CGGA validation cohort (G). (G) K-M survival curve in the Huanhu validation cohort. (H) Representative images of immunohistochemical staining in low-risk and high-risk patients. *, P＜0.05; ****, P＜0.0001.

Fig. 3

More importantly, to further validate the role of the five-gene risk signature, we examined the prognostic value of the risk score model using an independent cohort recruited from Tianjin Huanhu Hospital. This cohort included 50 GBM patients with full survival information. The protein levels of the five genes were detected by IHC. Patients were separated into high-risk and low-risk groups according to the same formula and the cut-off value calculated in TCGA-GBM cohort. The OS of patients in the high-risk group and low-risk group differed significantly; the high-risk group exhibited a poorer OS (HR: 2.109; 95 % CI: 1.085–4.099; P = 0.0278) (Fig. 3G). Representative IHC images of MAP1LC3A, MPO, GATA4, S100A4, and SPP1 are shown in Fig. 3H. All the results suggested that the risk score calculated using MAP1LC3A, MPO, GATA4, S100A4, and SPP1 was accurate in predicting the prognosis of patients with GBM.

3.4 Construction and evaluation of a nomogram

To determine whether risk score was an independent prognostic indicator, we performed multivariate Cox regression including age, sex, IDH mutation, ATRX mutation, TERT promoter mutation, MGMT promoter methylation, and risk score in TCGA-GBM (Fig. 4A). The results showed that the lack of methylation of the MGMT promoter (HR = 1.61, P = 0.037), age (HR = 1.02, P = 0.008), and risk score (HR = 355.28, 95 % CI: 2.02–4.98, P＜0.001) were significantly correlated with OS. Samples with a high-risk score had a greater risk of death. These results highlighted that risk score was independent of other clinical factors in TCGA-GBM cohort with promising predictive value.Fig. 4 Independent prognostic factor analysis and a nomogram predicting the 1-, 3-, and 5-year survival probability of GBM patients. (A) Multivariate Cox regression was used to analyze 7 factors in TCGA-GBM cohort. The status of MGMT promoter methylation, the risk score, and age were associated with OS (***, P＜0.001; **, P＜0.01; *, P＜0.05). (B) A prognostic nomogram was constructed to predict the OS at 1-, 3-, and 5-years using the data from TCGA-GBM. (C-E) Calibration plots of the nomogram for predicting survival probability at 1-, 3-, and 5-years in TCGA-GBM cohort. (F-G) The ROC curve showing the AUC at 1-, 3- and 5- years in TCGA and CGGA database.

Fig. 4

Nomograms are a visible prognostic model. By integrating different prognostic variables, nomograms can be used to predict individual survival probability. It is possible to accurately determine the outcomes of GBM patients using these concise scoring systems. Therefore, nomograms are a popular prognostic model and are widely used in oncological settings for predicting prognosis. In the present study, a nomogram that integrated the above three independent prognostic factors was constructed (Fig. 4B). To estimate the OS of individuals, a vertical line was drawn upward to the point axis, and the sum of the points of the 3 variables was determined. Then, the 3 independent variables were added and the score on the total point axis was marked. By doing this, we calculated the 1-, 3-, and 5-year OS by drawing a straight line downward from the total point axis to the outcome axis. For example, a 20-year-old (0 points) GBM patient with a risk score (100 points) and methylated MGMT promoter (0 points) got a total of 100 points. The probability of 1-, 3-, and 5-year survival probability for this patient were 60 %, 10 %, and below 1 %, respectively. The risk score ranged from 0 to 100, which occupied the largest proportion of the risk points compared with the other two factors. This result was consistent with the multivariate Cox regression analysis above.

The prediction accuracy of the nomogram was evaluated by calibration curves. The calibration curve showed the bias between the predicted and the actual survival rate. A 45° ideal curve was used as a reference. The closer the actual curve is to the ideal curve, the better the effectiveness of the model. The 1-, 3-, and 5-year calibration curves were in good agreement with the ideal curve (Fig. 4C–E). Furthermore, ROC analysis showed that the 1-, 3-, and 5-year AUC in TCGA and CGGA datasets were 0.786, 0.661, 0.711 and 0.655, 0.648, 0.647, respectively (Fig. 4F and G). These results suggested that the nomogram model could effectively predict the survival of GBM patients.

3.5 Estimation of the immune landscape between the high-risk and low-risk GBM patients

CIBERSORT is the most popular machine-learning approach for immune cell infiltration analysis at present. It is a computational method with the aim of solving the immune cell composition of tumor samples and helpful for high throughput characterization analysis. LM22 is a signature matrix file containing 547 genes. These genes could distinguish 22 human hematopoietic phenotypes in vivo or in vitro cells, including NK cells, seven T cell types, naïve and memory B cells, plasma cells, and myeloid subsets. Here, to estimate the immune infiltration differences of immune cells between the high-risk and low-risk groups in TCGA-GBM cohort, we used CIBERSORT combined with the LM22 signature matrix. We estimated the immune cell infiltration, and the changes in the proportion of tumor-infiltrating immune cells in different patients may represent the intrinsic features of individual differences (Fig. 5A). There were significant differences in the seven types of immune cells between the high-risk and low-risk groups. These differentially expressed cells were memory B cells, plasma cells, resting NK cells, activated NK cells, M1 macrophages, resting mast cells, and neutrophils, respectively. GBM patients with high-risk scores had a significantly higher proportion of activated NK cells and resting mast cells, and a significantly lower proportion of memory B cells, plasma cells, resting NK cells, M1 macrophages, and neutrophils than the low-risk score patients (Fig. 5B). These results may represent a potential pathway between the high-risk and low-risk groups of prognostic differences. In addition, we conducted GSEA functional enrichment analysis between the high-risk and low-risk groups. The results revealed that there were 194 significant signaling pathways enriched in the high-risk group. The top three greatly enriched pathways were the MAPK signaling pathway, calcium signaling pathway, and cytokine−cytokine receptor interaction pathway (Fig. 5C; Table S7). In summary, all these results indicated that abnormal immune cell infiltration may be a prognostic indicator.Fig. 5 Immune cell proportions in the TME and GSEA analysis of the high-risk and low-risk GBM patients. (A) The CIBERSORT algorithm was utilized to analyze the immune cell infiltration levels. The heap map shows the relative proportion of 22 immune cell types in the GBM patients. (B) A total of 22 immune cell types were compared between the high-risk and low-risk groups using boxplots. (C) Significant enrichment of the GSEA pathways between high-risk and low-risk groups based on TCGA-GBM cohort.

Fig. 5

3.6 TIDE signature and identification of cold/hot tumors for immunotherapy

ICI treatment for cancer can bring long-term benefits for patients. However, only a small proportion of patients exhibit a favorable response. To further estimate the immunotherapy response in TCGA-GBM cohort, we employed a computational algorithm named the TIDE. The TIDE score may assist oncologists in selecting patients who are more likely to respond favorably to ICIs. We compared the differences in immune response between the high-risk and low-risk score groups (Fig. 6A). The results showed that there was a significant difference in the immune response in the two groups (Fig. 6B). It was also found that the TIDE score was higher in the low-risk score group (Fig. 6C). The TIDE score and risk score were negatively correlated with each other.Fig. 6 TIDE and cold/hot tumor analysis. (A) TIDE analysis was shown using a waterfall plot. The red and blue lines represent the tumors that would or would not respond to immunotherapy, respectively. (B and C): Comparison of the TIDE score between the high-risk and low-risk groups. (D and E): TCGA-GBM samples were classified into 2 clusters based on “hot” tumor expression levels. (F) Heat map showing the expression levels of 12 gene marker of cold and hot tumors in TCGA-GBM cohort. (G) Risk score of cold tumors and hot tumors in TCGA-GBM cohort.

Fig. 6

Depending on whether there was T cell infiltration, tumors were divided into cold tumors and hot tumors. Immunotherapy was not effective in cold tumors which lacked T-cell infiltration. Conversely, hot tumors with T-cell infiltration were associated with improved immunotherapy treatment outcomes. According to the hot tumor markers, the TCGA-GBM samples could cluster into two categories (Fig. 6D and E). The 12 hot tumor markers were CXCL9, CXCL10, CXCL11, CXCR3, CXCR4, CCL5, CD8A, CD8B, CD4, CD3E, CD274, and PDCD1, respectively. The expression levels of these 12 markers in hot tumors were higher than those in cold tumors (Fig. 6F). Meanwhile, the risk score of hot tumors was higher than of cold tumors (Fig. 6G). These results demonstrated that immunotherapy may be helpful for the proportion of GBM patients with a high TIDE score and high expression levels of these hot tumor markers.

4 Differences in metabolic reprogramming between high-risk and low-risk groups

Metabolic reprogramming is important characteristic for tumor, then we further explore the possible variation of metabolic reprogramming in our study. Here, 29 genes involved in glucose and lipid metabolism pathways were selected, and the differences between the high-risk and low-risk groups were analyzed (Fig. 7A). These 29 genes were higher in the high-risk group than in the low-risk group. Furthermore, in the enrichment analysis of KEGG, 186 metabolic pathways were analyzed (Table S8). The result indicated that 31 pathways were statistically significant between high-risk and low-risk groups (Fig. 7B–Table S9). These results suggested that the risk score model was meaningful in metabolic reprogramming. Further analysis was required to explore the potential mechanism.Fig. 7 Metabolic reprogramming differences between high-risk and low-risk groups. (A) Heat map showing the 29 differentially expressed genes related to glucose and lipid metabolism pathways in the high-risk and low-risk groups. (B)) KEGG pathway enrichment results.

Fig. 7

5 Discussion

5.1 The background of this study

Nowadays, gliomas are the most common primary malignancy of the CNS, and GBM is the most aggressive type of glioma. GBM patients have a poorer prognosis and shorter OS after surgery, chemotherapy, and radiotherapy due to their invasiveness and high likelihood of recurrence. Current therapeutic strategies have failed to improve OS in the past 40 years [33]. Thus, there is an urgent need to identify a novel approach to GBM management. Aging, macrophage, autophagy, and methylation are not independent factors, instead they work together to influence a tumor. Aging is an important risk factor for glioma development. The Nlrp3 inflammasome controls the process of brain aging, where it serves as a key link between aging and glioma progression [34]. Epigenetic aging has been reported as a valuable marker for tumor prognosis. The epigenetic aging pattern of the DNA methylation model is closely related to OS and recurrence in glioma [35]. Hypoxia is a microenvironment feature in glioma. Macrophages are one of the primary cells residing in the microenvironment. They can polarize into classically activated macrophages (M1-like type) and alternatively activated macrophages (M2-like type) in a hypoxic microenvironment. M2-like type macrophages play an important role in the development and progression of tumors. Autophagy promotes cell survival by providing essential nutrition during hypoxia. By enhancing autophagy, hypoxic glioma promotes M2-like macrophage polarization [36]. Autophagy-related genes (ARGs) associated with methylation were identified in TCGA. A total of six ARGs were related to methylation levels and prognosis in low-grade glioma, and the methylation regulation of ARGs played an important role in the immune microenvironment [37]. This study analyzed the DEGs in GBM samples and healthy controls in TCGA. A total of 256 overlapping genes were identified in the DEGs of aging, macrophage, autophagy, and methylation related genes. The genes were classified into 3 clusters based on K-means analysis. IDH mutation, ATRX mutation, TERT promoter mutation, and MGMT promoter methylation exerted significant differences in the 3 clusters. Cluster 3 had the worst OS compared with the other two clusters.

5.2 Identification and validation of the five genes

Firstly, using a series of univariate and multivariate regression analyses, the value of a five-gene risk signature was verified. The five genes were MAP1LC3A, MPO, GATA4, S100A4, and SPP1, respectively. The risk scores were calculated based on the expression levels of these five genes. Patients were divided into high-risk and low-risk group based on the median of risk score. Next, two other datasets, CGGA and the Huanhu cohort from our hospital were used for validation. MAP1LC3A, a protein maker for autophagy, includes two subtypes: MAP1A and MAP1B. These two proteins were all associated with microtubules. MAP1LC3A was reported to be involved in the carcinogenesis of multiple types of tumors, such as esophageal squamous carcinoma, glioma, gastric cancer, and osteosarcoma, amongst others [38]. A total of 43 % of the GBM samples exhibited upregulated expression of MAP1LC3A, and this upregulation was correlated with GBM carcinogenesis [39]. Previous studies have shown that MAP1LC3A is a risk signature of GBM. NRG1, ITGA3, and MAP1LC3A were selected to construct a prognostic risk score model [40]. MPO is a highly oxidative and inflammatory-associated enzyme, that is primarily expressed in inflammatory monocytes and neutrophils. H2O2 is a substrate of MPO, and several oxidative species are generated following its interaction. MPO serves both a pro- and anti-tumorigenic effect. MPO + macrophages are associated with long term survival of brain tumors based on single-cell technologies [41]. The GATA transcription factors include six family members. GATA4 is a member of the zinc finger transcription factor family of proteins. It serves as a negative OS marker in GBM patients by regulating the expression levels of the anti-apoptotic proteins, Bcl-2 and Bcl-xL. The methylation of GATA4 promoter is a survival prognostic marker in GBM [42]. S100A4 is an upstream regulator of EMT-associated factors. It exerts the function on tumor cell motility and metastasis by binding to calcium. Higher-grade tumors show higher expression levels of S100A4 than lower-grade tumors [43]. S100A can also act as a regulator of myeloid and immune suppressive T cells. The knockdown of S100A resulted in a failure of reprogramming of immune cells and thus improved OS significantly, and it has served as an immunotherapeutic target in GBM patients [44]. SPP1, also known as osteopontin (OPN), is an acidic secreted phosphoprotein with several functions. It is involved in the development of several tumors, including breast carcinoma, glioma, hepatocellular carcinoma, ovarian carcinoma, and cervical cancer, where it exhibits upregulated expression [45]. Using these five genes in GBM, we generated a risk score model.

5.3 Construction of a nomogram

Next, we analyzed seven factors including age, sex, risk score, IDH mutation, ATRX mutation, TERT promoter mutation, and MGMT promoter methylation to verify whether the risk score was an independent prognostic factor. The results showed that age, risk score, and MGMT promoter methylation were all significantly associated with OS. IDH mutation, ATRX mutation, TERT promoter mutation, and MGMT promoter methylation were all closely associated with prognosis. GBM patients with IDH mutations had a relatively favorable survival compared with GBM patients with wild-type IDH [46,47]. TERT promoter mutation was a characteristic of GBM and may thus be used as a diagnostic factor [48]. In addition, it was also found that the MGMT gene promoter methylation was a strong prognostic factor [49]. There were some correlations between IDH, ATRX, and TERT promoter mutations in glioma. IDH mutation has been reported not only in glioma but also in other types of tumors. One of the characteristics between primary and recurrent tumors is IDH mutation. The majority of IDH-mutant GBM patients showed ATRX loss. While GBM patients with wild-type IDH exhibited ATRX expression. The presence of mutant IDH also increased histone lysine methylation by activating the TERT promoter [50].

A nomogram consisting of age, risk score, and MGMT promoter methylation status was established. To the best of our knowledge, this nomogram is the first to combine aging, macrophage, autophagy, and methylation-related genes for prognostic prediction. Calibration curves and ROC analysis indicated that the nomogram was excellent for prognostic prediction. The nomogram could assist in predicting the survival probability of an individual. We found that the established model could accurately predict the survival of GBM patients in the TCGA and CGGA dataset in 1-, 3-, and 5-year, which is useful and meaningful in clinical application.

5.4 Analysis of TME and metabolic reprogramming

TME is a complex environment in which tumor cells can survive and develop, and the immune cells and their regulation in the microenvironment play an important role in the occurrence and development of tumors. Tumor-infiltrating leukocytes (TILs), which are vital components of the TME has been found to be closely related to the prognosis of tumors [51]. The assessment and valuation of TILs are important as their functions vary in different tumor types [52,53]. ICIs are important immunotherapeutic tools. The most well-known ICIs were developed based on programmed cell death protein 1/programmed cell death 1 ligand 1 (PD-1/PD-L1) inhibitors, which block the PD-1/PD-L1 signaling pathway between tumor cells and T cells, resulting in the activation of infiltrating T cells and the attack of tumor cells. For several types of late-stage tumors, ICI treatment is the first choice and has provided renewed hope for numerous patients [54]. However, only a few patients benefit from this treatment. Thus, researchers are attempting to promote the transformation of “cold tumors” into “hot tumors” to enhance the efficacy of ICIs. Tumors can be classified into immunologically inactive tumors (cold) tumors and active (hot) tumors based on immune cell infiltration in the TME. Cold tumors are insensitive to therapies such as chemotherapy and immunotherapy and thus have a poor prognosis. However, hot tumors are characterized by the infiltration of immune cells and thus have a better prognosis when treated with ICIs [55]. In the present study, GBM patients with high-risk scores had a significantly higher proportion of activated NK cells and resting mast cells, and a significantly lower proportion of memory B cells, plasma cells, resting NK cells, M1 macrophages, and neutrophils than the low-risk score patients. GBM samples were divided into cold and hot tumors based on the immune infiltration score. The risk score in hot tumors was higher than in cold tumors. Our results showed that the TME changed between high-risk and low-risk groups.

Metabolic reprogramming is a characteristic of tumors. Glucose, lipid, amino, and lactate metabolism exert a significant impact on the TME, migration, invasion, and angiogenesis. Metabolic reprogramming is closely related to the TME. The tumor cell-derived metabolites affect the composition and distribution of immune cells in the TME, ultimately resulting in immune disorder and tumor progression [56]. Tumor cells alter their metabolic mechanisms to promote growth and survival. These altered metabolisms have common characteristics: including increased uptake of glucose and fermentation of glucose to lactate. This phenomenon is known as the Warburg Effect. The Warburg Effect has been documented for over 90 years and is one of the most extensively studied phenomena in almost all kinds of tumors [57]. Here, 29 genes associated with glucose metabolism and lipid metabolism were selected, and it was found that their expression levels were higher in the high-risk group than in the low-risk group. Moreover, 31 metabolism-related pathways were statistically significant between high-risk and low-risk groups. Thus, the exact function and mechanisms of these genes or pathways, alone and in combination, is deserved of further study. In our work, GBM samples were divided into two classes by the value of risk score. Based on this classification, immune infiltration and metabolic reprogramming showed significantly difference in a lot of genes and pathways. Our work indicated that TME and metabolic reprogramming had impact on the prognosis of GBM. These results are inspiring and encourage us in further exploring the mechanism of GBM. The workflow of this study is shown in Fig. 8.Fig. 8 The workflow of this study.

Fig. 8

5.5 Strengths and limitations

There are some strengths in our work. First of all, the novelty aims at the comprehensive consideration of aging, autophagy, macrophage and methylation on GBM. The risk score calculated using MAP1LC3A, MPO, GATA4, S100A4, and SPP1 was accurate in predicting the prognosis of patients with GBM in TCGA, CGGA and our own database. And we successfully constructed a prognosis model combined risk score and clinical parameters. What's more, immune infiltration and metabolic reprogramming were further analyzed between high-risk and low-risk groups.

However, the present study also has some limitations. First, the data available from the two public databases and the Huanhu cohort were limited. The value of the risk score for clinical application requires further validation in expanded clinical trials to confirm the results of the present study. The public database makes it difficult to collect complete clinical information for each patient. Second, the cut-off value was based on the median risk score in this study. Thus, an optimal cut-off value should be determined instead of the median of the risk score. Third, our novel risk score model did not take into consideration when GBM patients received different treatment therapies. Fourth, it is better to further explore the function of these five risk genes by in vitro and in vivo experiments, such as CCK-8, colony formation, transwell migration and invasion assays, and animal experiments. Lastly, given the retrospective design, potential bias was associated with unbalanced clinicopathological features which cannot be ignored. Prospective, large-scale cohort studies are also needed to validate the prognostic value of MAP1LC3A, MPO, GATA4, S100A4, and SPP1that have been screened out in our study.

6 Conclusion

In summary, we developed a novel five-gene risk score model for predicting the prognosis of GBM in a TCGA cohort. CGGA data and the Huanhu cohort were used for validation. Next, we established a novel prognostic nomogram by age, risk score, and MGMT promoter methylation status. We analyzed the immune landscape and metabolic reprogramming in GBM. According to this classification of GBM, the difference in immune infiltration and metabolic reprogramming is meaningful, which indicates that our prognostic model has been constructed successfully. Our results may be beneficial for predicting the prognosis of glioma and evaluating the effect of immunotherapy. Our study provides novel insights into the prognosis and immunotherapy for GBM.

Funding

This study was supported by the 10.13039/501100001809 National Natural Science Foundation of China (grant number: 81972349 and 82272759 ).

Data availability statement

All data generated or analyzed during this study are available from the corresponding author upon reasonable request. The data associated with this study has not been deposited into a publicly available repository.

CRediT authorship contribution statement

Le Zhang: Writing – review & editing, Software, Methodology. Xiaoling Yan: Investigation, Formal analysis. Yahong Wang: Validation, Resources. Qin Wang: Writing – original draft, Software. Hua Yan: Visualization, Funding acquisition. Yan Yan: Supervision, Funding acquisition, Conceptualization.

Declaration of competing interest

The authors declare no competing interests.

Appendix A Supplementary data

The following are the Supplementary data to this article:Multimedia component 1

Multimedia component 1

figs1 figs1

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.heliyon.2024.e37374.
==== Refs
References

1 Xu S. Tang L. Li X. Fan F. Liu Z. Immunotherapy for glioma: current management and future application Cancer Lett. 476 2020 1 12 10.1016/j.canlet.2020.02.002 32044356
2 Louis D.N. Perry A. Wesseling P. Brat D.J. Cree I.A. Figarella-Branger D. Hawkins C. Ng H.K. Pfister S.M. Reifenberger G. Soffietti R. von Deimling A. Ellison D.W. The 2021 WHO classification of tumors of the central nervous system: a summary Neuro Oncol. 23 2021 1231 1251 10.1093/neuonc/noab106 34185076
3 Thakkar J.P. Dolecek T.A. Horbinski C. Ostrom Q.T. Lightner D.D. Barnholtz-Sloan J.S. Villano J.L. Epidemiologic and molecular prognostic review of glioblastoma Cancer Epidemiol. Biomarkers Prev. 23 2014 1985 1996 10.1158/1055-9965.Epi-14-0275 25053711
4 Asthagiri A.R. Pouratian N. Sherman J. Ahmed G. Shaffrey M.E. Advances in brain tumor surgery Neurol. Clin. 25 2007 975 1003 10.1016/j.ncl.2007.07.006 viii-ix 17964023
5 Xiao G. Zhang X. Zhang X. Chen Y. Xia Z. Cao H. Huang J. Cheng Q. Aging-related genes are potential prognostic biomarkers for patients with gliomas Aging (Albany NY) 13 2021 13239 13263 10.18632/aging.203008 33946049
6 Chen P. Zhao D. Li J. Liang X. Li J. Chang A. Henry V.K. Lan Z. Spring D.J. Rao G. Wang Y.A. DePinho R.A. Symbiotic macrophage-glioma cell interactions reveal synthetic lethality in PTEN-null glioma Cancer Cell 35 2019 868 884.e866 10.1016/j.ccell.2019.05.003 31185211
7 Belyaeva E. Kharwar R.K. Ulasov I.V. Karlina I. Timashev P. Mohammadinejad R. Acharya A. Isoforms of autophagy-related proteins: role in glioma progression and therapy resistance Mol. Cell. Biochem. 477 2022 593 604 10.1007/s11010-021-04308-w 34854022
8 Pan T. Wu F. Li L. Wu S. Zhou F. Zhang P. Sun C. Xia L. The role m(6)A RNA methylation is CNS development and glioma pathogenesis Mol. Brain 14 2021 119 10.1186/s13041-021-00831-5 34281602
9 Mu J. Gong J. Shi M. Zhang Y. Analysis and validation of aging-related genes in prognosis and immune function of glioblastoma BMC Med. Genom. 16 2023 109 10.1186/s12920-023-01538-3
10 Wang C. Qiu J. Chen S. Li Y. Hu H. Cai Y. Hou L. Prognostic model and nomogram construction based on autophagy signatures in lower grade glioma J. Cell. Physiol. 236 2021 235 248 10.1002/jcp.29837 32519365
11 Li J. Guo Q. Xing R. Construction and validation of an immune infiltration-related risk model for predicting prognosis and immunotherapy response in low grade glioma BMC Cancer 23 2023 727 10.1186/s12885-023-11222-5 37543576
12 Hou W. Cai J. Shen P. Zhang S. Xiao S. You P. Tong Y. Li K. Qi Z. Luo H. Identification of FXYD6 as the novel biomarker for glioma based on differential expression and DNA methylation Cancer Med. 12 2023 22170 22184 10.1002/cam4.6752 38093622
13 Zhou J. Guo Y. Fu J. Chen Q. Construction and validation of a glioma prognostic model based on immune microenvironment Neuroimmunomodulation 29 2022 402 413 10.1159/000522529 35354148
14 Bao J.H. Lu W.C. Duan H. Ye Y.Q. Li J.B. Liao W.T. Li Y.C. Sun Y.P. Identification of a novel cuproptosis-related gene signature and integrative analyses in patients with lower-grade gliomas Front. Immunol. 13 2022 933973 10.3389/fimmu.2022.933973
15 Xu H. Zhang A. Han X. Li Y. Zhang Z. Song L. Wang W. Lou M. ITGB2 as a prognostic indicator and a predictive marker for immunotherapy in gliomas Cancer Immunol. Immunother. 71 2022 645 660 10.1007/s00262-021-03022-2 34313821
16 Zuo S. Zhang X. Wang L. A RNA sequencing-based six-gene signature for survival prediction in patients with glioblastoma Sci. Rep. 9 2019 2615 10.1038/s41598-019-39273-4 30796273
17 Long J. Wang A. Bai Y. Lin J. Yang X. Wang D. Yang X. Jiang Y. Zhao H. Development and validation of a TP53-associated immune prognostic model for hepatocellular carcinoma EBioMedicine 42 2019 363 374 10.1016/j.ebiom.2019.03.022 30885723
18 Lambrechts D. Wauters E. Boeckx B. Aibar S. Nittner D. Burton O. Bassez A. Decaluwé H. Pircher A. Van den Eynde K. Weynand B. Verbeken E. De Leyn P. Liston A. Vansteenkiste J. Carmeliet P. Aerts S. Thienpont B. Phenotype molding of stromal cells in the lung tumor microenvironment Nat. Med. 24 2018 1277 1289 10.1038/s41591-018-0096-5 29988129
19 Jiang W. He Y. He W. Wu G. Zhou X. Sheng Q. Zhong W. Lu Y. Ding Y. Lu Q. Ye F. Hua H. Exhausted CD8+T cells in the tumor immune microenvironment: new pathways to therapy Front. Immunol. 11 2020 622509 10.3389/fimmu.2020.622509
20 Faubert B. Solmonson A. DeBerardinis R.J. Metabolic reprogramming and cancer progression Science 368 2020 10.1126/science.aaw5473
21 Robinson M.D. McCarthy D.J. Smyth G.K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data Bioinformatics 26 2010 139 140 10.1093/bioinformatics/btp616 19910308
22 Sabatti C. Service S. Freimer N. False discovery rate in linkage and association genome screens for complex disorders Genetics 164 2003 829 833 10.1093/genetics/164.2.829 12807801
23 Yu G. Wang L.G. Han Y. He Q.Y. clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters vol. 16 2012 Omics 284 287 10.1089/omi.2011.0118
24 Subramanian A. Tamayo P. Mootha V.K. Mukherjee S. Ebert B.L. Gillette M.A. Paulovich A. Pomeroy S.L. Golub T.R. Lander E.S. Mesirov J.P. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles Proc. Natl. Acad. Sci. U. S. A. 102 2005 15545 15550 10.1073/pnas.0506580102 16199517
25 Szklarczyk D. Gable A.L. Lyon D. Junge A. Wyder S. Huerta-Cepas J. Simonovic M. Doncheva N.T. Morris J.H. Bork P. Jensen L.J. Mering C.V. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets Nucleic Acids Res. 47 2019 D607 d613 10.1093/nar/gky1131 30476243
26 Shannon P. Markiel A. Ozier O. Baliga N.S. Wang J.T. Ramage D. Amin N. Schwikowski B. Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks Genome Res. 13 2003 2498 2504 10.1101/gr.1239303 14597658
27 Friedman J. Hastie T. Tibshirani R. Regularization paths for generalized linear models via coordinate descent J. Stat. Software 33 2010 1 22
28 Alba A.C. Agoritsas T. Walsh M. Hanna S. Iorio A. Devereaux P.J. McGinn T. Guyatt G. Discrimination and calibration of clinical prediction models: users' guides to the medical literature JAMA 318 2017 1377 1384 10.1001/jama.2017.12126 29049590
29 Newman A.M. Liu C.L. Green M.R. Gentles A.J. Feng W. Xu Y. Hoang C.D. Diehn M. Alizadeh A.A. Robust enumeration of cell subsets from tissue expression profiles Nat. Methods 12 2015 453 457 10.1038/nmeth.3337 25822800
30 Fu J. Li K. Zhang W. Wan C. Zhang J. Jiang P. Liu X.S. Large-scale public data reuse to model immunotherapy response and resistance Genome Med. 12 2020 21 10.1186/s13073-020-0721-z 32102694
31 Qiang H. Li J. Chang Q. Shen Y. Qian J. Chu T. Mining GEO and TCGA database for immune microenvironment of lung squamous cell carcinoma patients with or without chemotherapy Front. Oncol. 12 2022 835225 10.3389/fonc.2022.835225
32 Liberzon A. Birger C. Thorvaldsdóttir H. Ghandi M. Mesirov J.P. Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection Cell Syst 1 2015 417 425 10.1016/j.cels.2015.12.004 26771021
33 Khosla D. Concurrent therapy to enhance radiotherapeutic outcomes in glioblastoma Ann. Transl. Med. 4 2016 54 10.3978/j.issn.2305-5839.2016.01.25 26904576
34 Li L. Liu Y. Aging-related gene signature regulated by Nlrp3 predicts glioma progression Am. J. Cancer Res. 5 2015 442 449 25628952
35 Liao P. Ostrom Q.T. Stetson L. Barnholtz-Sloan J.S. Models of epigenetic age capture patterns of DNA methylation in glioma associated with molecular subtype, survival, and recurrence Neuro Oncol. 20 2018 942 953 10.1093/neuonc/noy003 29432558
36 Xu J. Zhang J. Zhang Z. Gao Z. Qi Y. Qiu W. Pan Z. Guo Q. Li B. Zhao S. Guo X. Qian M. Chen Z. Wang S. Gao X. Zhang S. Wang H. Guo X. Zhang P. Zhao R. Xue H. Li G. Hypoxic glioma-derived exosomes promote M2-like macrophage polarization by enhancing autophagy induction Cell Death Dis. 12 2021 373 10.1038/s41419-021-03664-1 33828078
37 Qiao Q. Wang Y. Zhang R. Pang Q. Autophagy related DNA methylation signature predict clinical prognosis and immune microenvironment in low-grade glioma Transl. Cancer Res. 11 2022 2157 2174 10.21037/tcr-22-310 35966301
38 Bai H. Inoue J. Kawano T. Inazawa J. A transcriptional variant of the LC3A gene is involved in autophagy and frequently inactivated in human cancers Oncogene 31 2012 4397 4408 10.1038/onc.2011.613 22249245
39 Giatromanolaki A. Sivridis E. Mitrakas A. Kalamida D. Zois C.E. Haider S. Piperidou C. Pappa A. Gatter K.C. Harris A.L. Koukourakis M.I. Autophagy and lysosomal related protein expression patterns in human glioblastoma Cancer Biol. Ther. 15 2014 1468 1478 10.4161/15384047.2014.955719 25482944
40 Wang Z. Gao L. Guo X. Feng C. Lian W. Deng K. Xing B. Development and validation of a nomogram with an autophagy-related gene signature for predicting survival in patients with glioblastoma Aging (Albany NY) 11 2019 12246 12269 10.18632/aging.102566 31844032
41 Karimi E. Yu M.W. Maritan S.M. Perus L.J.M. Rezanejad M. Sorin M. Dankner M. Fallah P. Doré S. Zuo D. Fiset B. Kloosterman D.J. Ramsay L. Wei Y. Lam S. Alsajjan R. Watson I.R. Roldan Urgoiti G. Park M. Brandsma D. Senger D.L. Chan J.A. Akkari L. Petrecca K. Guiot M.C. Siegel P.M. Quail D.F. Walsh L.A. Single-cell spatial immune landscapes of primary and metastatic brain tumours Nature 614 2023 555 563 10.1038/s41586-022-05680-3 36725935
42 Vaitkienė P. Skiriutė D. Skauminas K. Tamašauskas A. GATA4 and DcR1 methylation in glioblastomas Diagn. Pathol. 8 2013 7 10.1186/1746-1596-8-7 23320456
43 Zhang Y. Yang X. Zhu X.L. Bai H. Wang Z.Z. Zhang J.J. Hao C.Y. Duan H.B. S100A gene family: immune-related prognostic biomarkers and therapeutic targets for low-grade glioma Aging (Albany NY) 13 2021 15459 15478 10.18632/aging.203103 34148033
44 Abdelfattah N. Kumar P. Wang C. Leu J.S. Flynn W.F. Gao R. Baskin D.S. Pichumani K. Ijare O.B. Wood S.L. Powell S.Z. Haviland D.L. Parker Kerrigan B.C. Lang F.F. Prabhu S.S. Huntoon K.M. Jiang W. Kim B.Y.S. George J. Yun K. Single-cell analysis of human glioma and immune cells identifies S100A4 as an immunotherapy target Nat. Commun. 13 2022 767 10.1038/s41467-022-28372-y 35140215
45 Liu Y. Ye G. Dong B. Huang L. Zhang C. Sheng Y. Wu B. Han L. Wu C. Qi Y. A pan-cancer analysis of the oncogenic role of secreted phosphoprotein 1 (SPP1) in human cancers Ann. Transl. Med. 10 2022 279 10.21037/atm-22-829 35433956
46 Molenaar R.J. Radivoyevitch T. Maciejewski J.P. van Noorden C.J. Bleeker F.E. The driver and passenger effects of isocitrate dehydrogenase 1 and 2 mutations in oncogenesis and survival prolongation Biochim. Biophys. Acta 1846 2014 326 341 10.1016/j.bbcan.2014.05.004 24880135
47 Wang Z. Bao Z. Yan W. You G. Wang Y. Li X. Zhang W. Isocitrate dehydrogenase 1 (IDH1) mutation-specific microRNA signature predicts favorable prognosis in glioblastoma patients with IDH1 wild type J. Exp. Clin. Cancer Res. 32 2013 59 10.1186/1756-9966-32-59 23988086
48 Galbraith K. Snuderl M. Molecular pathology of gliomas Surg Pathol Clin 14 2021 379 386 10.1016/j.path.2021.05.003 34373090
49 Gilbert M.R. Wang M. Aldape K.D. Stupp R. Hegi M.E. Jaeckle K.A. Armstrong T.S. Wefel J.S. Won M. Blumenthal D.T. Mahajan A. Schultz C.J. Erridge S. Baumert B. Hopkins K.I. Tzuk-Shina T. Brown P.D. Chakravarti A. Curran W.J. Jr. Mehta M.P. Dose-dense temozolomide for newly diagnosed glioblastoma: a randomized phase III clinical trial J. Clin. Oncol. 31 2013 4085 4091 10.1200/jco.2013.49.6968 24101040
50 Ohba S. Kuwahara K. Yamada S. Abe M. Hirose Y. Correlation between IDH, ATRX, and TERT promoter mutations in glioma Brain Tumor Pathol. 37 2020 33 40 10.1007/s10014-020-00360-4 32227259
51 Chen B. Khodadoust M.S. Liu C.L. Newman A.M. Alizadeh A.A. Profiling tumor infiltrating immune cells with CIBERSORT Methods Mol. Biol. 1711 2018 243 259 10.1007/978-1-4939-7493-1_12 29344893
52 Fridman W.H. Pagès F. Sautès-Fridman C. Galon J. The immune contexture in human tumours: impact on clinical outcome Nat. Rev. Cancer 12 2012 298 306 10.1038/nrc3245 22419253
53 Gentles A.J. Newman A.M. Liu C.L. Bratman S.V. Feng W. Kim D. Nair V.S. Xu Y. Khuong A. Hoang C.D. Diehn M. West R.B. Plevritis S.K. Alizadeh A.A. The prognostic landscape of genes and infiltrating immune cells across human cancers Nat. Med. 21 2015 938 945 10.1038/nm.3909 26193342
54 Wang T. Luo Y. Zhang Q. Shen Y. Peng M. Huang P. Zhou Z. Wu X. Chen K. COX-2-related tumor immune microenvironment in non-small cell lung cancer: a novel signature to predict hot and cold tumor J. Thorac. Dis. 14 2022 729 740 10.21037/jtd-22-257 35399232
55 Yang H. Zhao L. Zhang Y. Li F.F. A comprehensive analysis of immune infiltration in the tumor microenvironment of osteosarcoma Cancer Med. 10 2021 5696 5711 10.1002/cam4.4117 34258887
56 He J. Chen Z. Xue Q. Sun P. Wang Y. Zhu C. Shi W. Identification of molecular subtypes and a novel prognostic model of diffuse large B-cell lymphoma based on a metabolism-associated gene signature J. Transl. Med. 20 2022 186 10.1186/s12967-022-03393-9 35468826
57 Liberti M.V. Locasale J.W. The Warburg effect: how does it benefit cancer cells? Trends Biochem. Sci. 41 2016 211 218 10.1016/j.tibs.2015.12.001 26778478
