The transcriptomic data and clinical data of 970 glioma patients were downloaded from the Chinese Glioma Genome Atlas (CGGA) database (http://www.cgga.org.cn/) as the training cohort. For the validation cohort, transcriptomic data and the clinical data of 666 glioma patients were obtained from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). Samples with gene expression 0 were excluded. A total of 1811 immune-related genes were downloaded from the Immunology Database and Analysis Portal (IMMPORT; https://www.immport.org/).

Based on gene expression data, the ESTIMATE algorithm was used to calculate the stromal and immune scores for each tumor sample, predicting the levels of stromal and immune cells. The stromal score and immune score calculated by this algorithm helped quantify the stromal and immune cells in the tumor. The glioma patients were then divided into high- and low-stromal score subgroups and high- and low-immune score subgroups based on the median stromal score and immune score.

The R package limma was used to identify DEGs between high- and low-stromal score subgroups, as well as high and low immune scores. The cut-off criteria were set as |fold change| greater than 2 and p value less than 0.05. To reduce the false positive rate, the eBayes test was used to adjust the p value. The heatmaps and volcano plots of the DEGs were displayed using the ComplexHeatmap and “ggplot2” packages in R, respectively. Shared DEGs were analyzed using immune gene sets, DEGs of the stromal score group and the immune score group. A Venn diagram to display the crossover genes was constructed using the R package “venn”.

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were used to identify the characteristic biological attributes and functional attributes of DEGs, respectively. GO and KEGG analyses were performed using the R package “clusterProfiler” and “org.Hs.eg.db”. Enrichment results were visualized using the R package “ggplot2”, and a p-value of less than 0.05 and a false discovery rate (FDR) of less than 0.05 were considered as the criteria for significant enrichment.

The shared DEGs with p-values less than 0.05 were screened by univariate Cox analysis. The screened immune-related DEGs were further selected by Lasso regression with the R package “SIS” to identify the final variables used to construct the Cox model. A multivariate Cox regression risk model was then constructed using the final 21 DEGs. The results of the multivariate Cox analysis were presented with the R package “forestplot”. The risk score of the Cox model was calculated based on a linear combination of coefficients and gene expression levels. Its formula is as follows: Risk score = Coef 1 * Exp 1 + Coef 2 * Exp 2 + … + Coef n * Exp n . (1) Here, n represents the gene number. Exp represents the gene expression level. coef represents the coefficient value. The median risk score was set as the cut-off value, and all glioma patients were divided into high-risk and low-risk score groups. Receiver operating characteristic (ROC) curves of the multivariate Cox analysis were used to evaluate the performance of the model for the CGGA data and TCGA data. The model with an area under the curve (AUC) value of greater than 0.7 was considered a reliable model. Risk curves and survival status scatter plots were drawn using the risk scores and survival status of each patient, with the median risk value taken as the cut-off value, which were used to evaluate the predictive effect of the model on the survival prognosis of patients. The nomogram and calibration curves were generated using the “rms” package in R. The decision curve analysis (DCA) was conducted using the “ggDCA” package, specifically designed for creating decision curves.

Correlations were analyzed between several key variables and overall survival: the stromal score, the immune score, the expression levels of DEGs, the risk score, and the infiltration level of immune cells. The impact of these variables on the overall survival was calculated using Kaplan–Meier survival curves. For these analyses, two R packages, namely, “survival” and “survminer,” were employed. A p value less than 0.05 was considered statistically significant.

The difference in the risk score between the different genders, WHO grades, and IDH1 mutation states was compared. The statistical significance of the differences between the two groups was estimated using the Wilcoxon test, and the differences between the three groups were assessed by the Kruskal–Wallis test. Univariate and multivariate Cox regression analyses of age, gender, WHO grade, IDH1 mutation status, and risk score were performed and displayed using the R package forestplot. Factors with a p value less than 0.05 were considered to be independent prognostic factors.

According to gene expression data, the CIBERSORT algorithm was used to assess the immune infractions in each tumor sample and to provide an estimation of the abundances of each immune cell type in a mixed-cell population. The CIBERSORT algorithm was based on a known reference set that provided gene expression signatures for 22 leukocyte subtypes (LM22). The CIBERSORT algorithm was used to quantify the immune cells between high- and low-risk score groups. The results of inferred scores for immune cell populations were considered accurate at a threshold of p value 0.05. The distribution of immune cells in the two groups was shown using the ggpubr package. The univariate Cox analysis was performed using the plyr package and visualized using the forestplot package.

In this study, all statistical computations and figures were performed using R software (version 3.6.3). The Wilcoxon rank-sum test was used as a non-parametric test for comparison between the two groups, and the Kruskal–Wallis test was used to estimate the differences between three groups. The Wald test was used in the univariate and multivariate Cox analyses. The log-rank test was used to generate p-values in Kaplan–Meier survival analysis. A p value less than 0.05 was considered to be statistically significant.

In this study, we obtained gene expression data and clinical information of glioma patients from the CGGA database for the training cohort and the TCGA database for the validation cohort. After excluding patients with incomplete clinical information and normal controls, 970 and 666 glioma patients with gliomas were included in this study (Table 1). These glioma patients were diagnosed pathologically at ages 8–89 years, with a median age of 44 years. Men accounted for 58.37%, and these patients were at different WHO grades (31.48% WHO II, 35.57% WHO III, and 32.64% WHO IV). To evaluate the abundance of stromal cells and immune cells in the tumor microenvironment, we calculated the stromal score and the immune score of all patients in the training cohort using ESTIMATE. The results showed that the stromal score ranged from −2567.35 to 1358.71 and the immune score ranged from −2557.44 to 2611.02 (Figures 1A, B). We then examined the relationship between stromal/immune scores and the different WHO grades in the CGGA database. The results displayed a positive correlation between stromal/immune scores and WHO stages, indicating a higher prevalence of advanced-grade patients within the high-stromal score and high-immune score groups (Supplementary Tables S1, S2). To investigate the association between the stromal score, immune score, and overall survival, we divided 970 glioma patients into a high-score group and a low-score group and performed a Kaplan–Meier survival analysis between the two groups. The results displayed that the median overall survival of glioma patients with a low stromal score (−1400.96) was longer than those with a high stromal score (−677.93), and the median overall survival of low-immune score (−589.07) patients was longer than that of high-immune score (403.25) patients (p <0.0001) (Figures 1C, D). These results indicated that stromal score and immune score could be used as prominent predictors of overall survival in patients with gliomas.

To determine key factors in the immune microenvironment that affect prognosis, we analyzed immune-related genes that were differentially expressed in gliomas. First, we compared the DEGs between the high-score group and the low-score group with RNA-seq data of the training cohort, and we identified 9,589 DEGs from the stromal score groups and 9,777 DEGs from the immune score groups. The heatmaps and volcano plots of DEGs are shown in Figures 2A–D. The shared genes between the DEGs of the stromal score groups, the DEGs of immune score groups, and immune genes were analyzed. Among the 1,811 immune-related genes obtained from the ImmPort database, 667 genes were shared with the DEGs of the stromal score groups and 670 were shared with the DEGs of the immune score groups. A total of 610 immune-related DEGs were identified between the stromal and immune score groups (Figure 2E). To explore the functions and enrichment pathways of these immune-related DEGs, GO and KEGG analyses were performed on 610 DEGs. As for GO analysis, we analyzed three subontologies of DEGs: cellular components, molecular functions, and biological processes. We exhibited the top 20 significantly enriched GO terms (Figure 2F). The results indicated that the external side of the plasma membrane and the MHC protein complex were the main enriched cellular components of immune DEGs. The molecular functions of the immune DEGs were primarily focused on receptor–ligand activity and signaling receptor activator activity. The biological process of the immune DEGs mainly included the pathways of leukocyte migration, cell chemotaxis, and positive regulation of leukocyte activation. Subsequently, the results of KEGG pathway enrichment analysis indicated that the DEGs were mainly involved in the pathway of the cytokine–cytokine receptor interaction and the Epstein–Barr virus infection (Figure 2G).

To clarify the potential value of immune-related DEGs in the prognosis of glioma patients, we used univariate Cox regression and Lasso regression analysis to screen for the most relevant DEGs. Finally, 21 immune-related DEGs were screened out from the 610 DEGs. A multivariate Cox regression risk model was constructed with the 21 genes (Supplementary Figure S1). We further investigated whether these 21 genes independently affected the overall survival of glioma patients. The Kaplan–Meier survival curves exhibited that the high expression level of CHGB, BMP2, SSTR2, IL17D, ADCYAP1R1, UCN, GDF10, and ARRB1 was significantly correlated with long overall survival (Figure 3A), indicating that these genes are good prognostic genes. Meanwhile, a significant negative association between the expression levels of TMSB4X, IFNGR2, SERPINA3, VIM, JUN, TXLNA, CDC42, HDAC1, PSMC2, IL-13RA2, PAK4, PDK1, and TMSB15A and overall survival is observed in Figure 3B. We then divided the patients into LGG and GBM sets and explored the relationship between the 21 genes and survival rates. The results indicated that the 21-gene signature exhibited significant correlations with the OS in LGG patients, revealing discernible patterns in gene expression, potentially indicative of disease progression (Supplementary Figure S2). However, many genes within our signature lacked a marked impact on OS in GBM subsets (Supplementary Figure S3). This suggested the inherent complexity and heterogeneity of GBM, which might attenuate the influence of individual genes or challenge the signature’s ability to encapsulate the nuanced biology of this subtype.

We developed a prognostic risk model using 21 immune-related DEGs to evaluate the prognosis of glioma patients. The risk scores were calculated based on gene expression levels and their corresponding regression coefficients. The receiver operating characteristic plot for the prognostic model over a 5-year period is shown in Figure 4A, and the area under the curve was 0.848. We then used TCGA data as a validation cohort, and the AUC value of the ROC plot was 0.846 (Figure 4B). These results suggest that the prognostic model based on 21 immune-related DEGs is effective in predicting the prognosis of glioma patients. To further evaluate the model’s efficacy, we divided glioma patients into high-risk and low-risk groups based on their median risk score. The risk score distribution and survival status of glioma patients in the CGGA database showed that the patients in the high-risk score group had more deaths (Figure 4C). Additionally, we observed a significant difference in overall survival between the high-risk and low-risk groups through Kaplan–Meier survival analysis (p <0.0001). The low-risk group had a good prognosis, while the high-risk group had a poor prognosis (Figure 4D). Furthermore, we analyzed the expression levels of 21 immune-related genes between the high-risk and low-risk groups. We found that genes associated with a good prognosis (CHGB, BMP2, SSTR2, IL17D, ADCYAP1R1, UCN, GDF10, and ARRB1) were highly expressed in the low-risk group. Conversely, the genes associated with a poor prognosis (TMSB4X, IFNGR2, SERPINA3, VIM, JUN, TXLNA, CDC42, HDAC1, PSMC2, IL-13RA2, PAK4, PDK1, and TMSB15A) were highly expressed in the high-risk group (Figure 4E). Subsequently, we developed a nomogram utilizing the clinicopathological parameters and risk scores to estimate the 3-year and 5-year survival probabilities for glioma patients in the CGGA database (Supplementary Figure S4A) and TCGA database (Supplementary Figure S4B). Furthermore, we conducted a decision curve analysis and analyzed calibration curves to assess the performance of our predictive nomogram. The DCA plots revealed favorable net benefits across the majority of threshold probabilities, specifically indicating optimal threshold probabilities between 0 and 0.92 for the CGGA database (Supplementary Figure S5A) and between 0 and 0.75 for the TCGA database (Supplementary Figure S5C). Additionally, the calibration curves corroborated the concordance between the predicted survival probabilities and the observed outcomes, further substantiating the model’s predictive accuracy (Supplementary Figures S5B, D). The results indicated that the risk score derived from the prognostic model emerged as a molecular indicator of poor prognosis in glioma patients.

We used the constructed risk model to predict the risk scores for different subgroups of glioma patients. When comparing the risk scores between female and male patients, there was no statistically significant difference (p = 0.37; Figure 5A). The WHO stage is positively associated with poor prognosis in glioma patients. When comparing the risk scores of patients with different WHO stages, including WHO II, WHO III, and WHO IV, the results showed that the risk score of WHO IV had a significant difference compared with WHO III (p <0.05) or WHO II (p <0.05). Moreover, the risk score of WHO III also had a significant difference with the WHO II risk score (p <0.05) determined by the Kruskal–Wallis test. A significant positive association between WHO grades and the risk score is displayed in Figure 5B. We also analyzed the risk scores of the IDH1 mutant and wildtype glioma patients. The results indicated that the IDH1 mutant group had significantly decreased risk scores (p <0.05), indicating a better prognosis for patients with IDH1 mutations (Figure 5C). To verify whether the risk scores derived from the prognostic model can be used as an independent prognostic factor, univariate and multivariate Cox regression analyses were performed with the CGGA dataset. From the forest plots, the hazard ratios (HRs) of the risk score were 2.718 (95% confidence interval (CI), 2.467–2.995, p <0.05) and 2.179 (95% CI, 1.943–2.444, p <0.05), respectively. This result suggests that risk score could be considered an independent prognostic factor for the survival of glioma patients. Meanwhile, WHO grade was also an independent prognostic factor, the HRs of WHO grade III were 2.761 (95% CI, 2.143–3.357, p <0.05) and 2.206 (95% CI, 1.705–2.856, p <0.05), respectively, and the HRs of WHO grade IV were 7.605 (95% CI, 5.949–9.723, p <0.05) and 3.155 (95% CI, 2.386–4.172, p <0.05), respectively (Figure 5D). These results confirm that the risk score and WHO staging grades III and IV can be used as independent clinical prognostic factors for patients with gliomas.

To clarify the impact of immune cells on glioma prognosis, we analyzed immune cell infiltration in the TME using the CIBERSORT algorithm and then compared the differences of immune cell presence in gliomas between the high- and low-risk score subgroups. The results showed that infiltration levels of memory B cells, plasma cells, naïve CD4⁺ T cells, activated NK cells, monocytes, resting dendritic cells, activated mast cells, and eosinophils were increased in the low-risk score group. In contrast, CD8⁺ T cells, regulatory T cells (Tregs), γδ T cells, M0 macrophages, M1 macrophages, M2 macrophages, activated dendritic cells, and neutrophils were increased in the high-risk score group (Figure 6A). These results suggested that while the former set of immune cells might be linked to a favorable prognosis, the latter may be associated with a poorer outcome. To further investigate the potential risk of infiltrating immune cells on survival, we performed univariate Cox regression analysis and calculated the HRs with 95% CI. M0 macrophages, M2 macrophages, and neutrophils emerged as independent risk factors for unfavorable glioma prognosis (HR >1; p <0.05). Conversely, CD4 naïve T cells, monocytes, memory B cells, resting dendritic cells, and activated NK cells appeared as potential protective factors (HR <1; p <0.05) (Figure 6B). Kaplan–Meier survival analysis was further performed in these immune cells. We found pronounced survival probability differences associated with nine immune cell types. Notably, increased levels of memory B cells, resting dendritic cells, monocytes, plasma cells, and naïve CD4 T cells correlated with prolonged survival, while high levels of activated dendritic cells, M0 macrophages, M2 macrophages, and Tregs were linked to shorter survival (Figure 6C). In conclusion, our study underscored the variability in immune cell infiltration across different risk score groups. Crucially, specific immune cell infiltration levels emerged as influential determinants in glioma prognosis.