The flowchart shows the RNA sequencing and clinical data sources used in this study (Figure 1): The Cancer Genome Atlas (TCGA; cancergenome.nih.gov/), Gene Expression Omnibus (GEO; ncbi.nlm.nih.gov/geo/) and China Glioma Genome Atlas (CGGA; cgga.org.cn/) databases. While the CGGA database shared the cohort of 1,018 glioma samples, the TCGA database provided the cohort of 703 samples (698 glioma and 5 adjacent normal tissues). The GEO database covered 4 cohorts, including the GSE108474 cohort (148 astrocytoma, 228 glioblastoma multiforme, 67 oligodendroglioma, and 28 normal brain tissues), GSE4290 cohort (26 astrocytoma, 81 glioblastoma, 50 oligodendroglioma, and 23 normal brain tissues), GSE4412 cohort (85 glioma samples) and GSE43378 cohort (50 glioma samples). In addition, the ImmPort (immport.org/home) and InnateDB (innatedb.ca/) databases provided us with a gene list of 6196 immune-related genes (IRGs).

Based on this list, we extracted the RNA sequencing data of 2365, 1879, 1879, 1638, 1879, and 1700 IRGs from the TCGA, GSE108474, GSE4290, GSE4412, GSE43378, and CGGA cohorts, respectively. To identify differentially expressed immune-related genes (DE-IRGs) between tumor and normal tissues from TCGA, we set |log2 fold change| (|log2FC|) >1 and false discovery rate (FDR) < 0.05 as the filter condition. Similarly, FDR <0.05 was used as a new filter condition to identify DE-IRGs from the GSE108474 and GSE4290 cohorts, respectively. Finally, the R package “venn” was used to visualize the overlapping process of DE-IRGs lists and IRGs lists.

The RNA sequencing values of 163 common DE-IRGs were extracted from the TCGA cohort to construct the co-expression network among them. After clustering the TCGA samples for eliminating free samples, we used the function pickSoftThreshold to select the best soft power β = 4 to construct the best scale-free network. Based on the formula: a i j = | S i j / β

(a_ij: adjacency matrix between gene i and gene j, S_ij: similarity matrix which is done by Pearson correlation of all gene pairs, β: softpower value), we created an adjacency matrix, and converted it to a topological overlap matrix (TOM) and a corresponding dissimilarity (1-TOM) (Yip and Horvath, 2007). Then, we took 1-TOM as the distance to cluster the genes, and built a dynamic pruning tree to identify the gene modules (Yip and Horvath, 2007). Finally, after merging similar modules with 75% similarity, we identified 3 modules. Likewise, we chose the best soft power β = 4 and β = 5 to identify two modules from the GSE108474 and GSE4290 cohorts, respectively. The turquoise modules of the GSE4290, GSE108474, and TCGA cohorts all exhibited the strongest positive correlation with tumor status. The genes of these three modules were extracted separately to obtain the common DE-IRGs highly associated with glioma.

Samples with complete overall survival (OS) and RNA sequencing data from 6 cohorts were extracted for subsequent analysis, respectively. To obtain corresponding PR-DE-IRGs in the TCGA, GSE4412, GSE43378 and CGGA cohorts, we performed univariate Cox analysis with a cutoff value of p < 0.05. At the same time, we also performed Kaplan-Meier survival analysis on the common PR-DE-IRGs in these four cohorts to explore the relationship between their expression and OS.

668 samples from TCGA cohort were randomly matched to the training set (n = 334) and test set (n = 334) on average. Lasso regression analysis can screen out highly relevant PR-DE-IRGs from the training set, thereby minimizing the risk of overfitting the screening features, and achieving the purpose of accurately predicting the patients’ prognosis (Lu et al., 2022). The optimal penalty parameter (λ) obtained by the minimum 10-fold cross validation was used to determine five PR-DE-IRGs and corresponding coefficients for constructing the prognostic model. We calculated each sample’s risk score from four cohorts by the formula: R i s k   s c o r e = ∑ ( P R − D E − I R G s   e x p r e s s i o n   v a l u e s × c o r r e s p o n d i n g   c o e f f i c i e n t )

CCGA set, GEO’s GSE4412 and GSE43378 sets, and TGGA’s training, test and whole sets were used for the verification of model’s predictive ability. We used the median risk score of samples from each cohort as a cutoff point to classify samples into high-risk and low-risk groups. Kaplan-Meier survival analysis was used to validate the model’s ability in differentiating the glioma patients’ prognosis. We also plotted the receiver operating characteristic (ROC) curves to evaluate the accuracy of the model in predicting prognosis. Univariate and multivariate Cox regression analyses further identified the role of risk score as an independent predictor of prognosis. In addition, we also performed principal component analysis (PCA) and t-distributed random neighborhood embedding (t-SNE) to evaluate the model’s ability to discretize samples through the R package “Rtsne".

To enrich for potential biological functions and pathways involved in different risk groups, we ran Gene Set Enrichment Analysis (GSEA) based on the R package “clusterProfiler”. In this process, we used “c5.go.v7.4.symbols.gmt” and “c2.cp.kegg.v7.4.symbols.gmt” as reference gene sets and set nominal p value < 0.05 as the filter condition.

First, the overall stromal and immune cell scores for each TCGA sample were calculated using the ESTIMATE algorithm (sourceforge.net/projects/estimateproject/). We compared their differences between different risk groups, and analyzed their correlations with risk score through Spearman. Next, based on the single-sample gene set enrichment analysis (ssGSEA) of the R packages “GSEABase” and “gsva”, we quantified the scores of 16 immune cells and 13 immune functions in each sample. After visualizing their distribution across all samples using a heatmap, similar differential expression and correlation analyses were again applied to explore their relationship with the model.

Six types of immune infiltration were identified in human tumors, which corresponded from tumor promoting to tumor suppressing respectively, namely C1 (wound healing), C2 (INF-g dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (TGF-b dominant) (Lin et al., 2021). The two-way ANOVA was used to explore association between risk score and different types of immune infiltration. To compare the distribution of C3, C4, and C5 subtypes between different risk groups, we ran a chi-square test.

The one-class logistic regression machine learning algorithm (OCLR) provided targeted training for stem cells (embryonic stem cells; induced pluripotent stem cells) and their differentiated ectoderm, mesoderm and endoderm progenitor cells (Malta et al., 2018). Malta et al. calculated stemness index (mDNAsi) for each TCGA sample by OCLR, ranging from low (zero) to high (one) (Malta et al., 2018). The mDNAsi data from 564 glioma samples from this study were used for our analysis. Next, we explored the mDNAsi difference between different risk groups/clinical subgroups as well as the correlation between mDNAsi and risk score.

After counting the somatic gene mutation data of gliomas from TCGA, we obtained the top 30 genes with the highest mutation frequency. The R package “GenVisR” was used to visualize the statistical results of these 30 genes’ mutation in different risk groups. A similar method was applied to the mutated DE-IRGs highly associated with glioma. While analyzing the association between tumor mutation burden (TMB) and model, we also explored the effect of TMB on prognosis. Based on the mutation status of IDH1, the samples were divided into wild and mutant groups. In addition to the differences in the expression of 5 PR-DE-IRGs and CD274 between the two groups, the differences in the risk score, survival probability, 16 immune cells and 13 immune functions were compared.

To further explore the possible mechanisms involved in glioma progression from the perspective of ceRNA, we selected BMP2 as the target gene to construct a comprehensive regulatory network composed of interconnected ceRNAs. After annotating the miRNA expression data from TCGA using mature miRNA file from mirbase database (mirbase.org/), we obtained a miRNA expression matrix for 535 samples (530 glioma and 5 adjacent normal tissues). The starBase database (starbase.sysu.edu.cn/) contains multiple target gene prediction programs (PITA, RNA22, miRmap, microT, miRanda, PicTar, and TargetScan). We selected miRNAs that appeared more than 2 times in all prediction programs as candidate miRNAs for BMP2. During this period, Cytoscape (v3.7.2) was used to map the co-expression network of all predicted miRNAs and BMP2. Based on the correlation analyses between these miRNAs and BMP2 expression (filter condition: correlation coefficient < −0.4 and p < 0.001), we further screened out the differentially expressed miRNAs between glioma and normal tissues (filter condition: |log2FC| >1 and p < 0.05). Ultimately, only 3 miRNAs (hsa-miR-365a-3p, hsa-let-7e-5p and hsa-miR-98-5p) showed significant survival differences in the further Kaplan-Meier survival analyses (filter condition: p < 0.05).

Next, starBase was used to predict candidate lncRNAs that may bind to hsa-miR-365a-3p, hsa-let-7e-5p and hsa-miR-98-5p, respectively. Similarly, the correlation analyses between these lncRNAs and hsa-miR-365a-3p/hsa-let-7e-5p/hsa-miR-98-5p expression (filter condition: correlation coefficient < −0.2 and p < 0.001), as well as their difference analyses between glioma and normal tissues (filter condition: p < 0.05) were run separately. Next, 13/13/11 lncRNAs that might bind to hsa-miR-365a-3p/hsa-let-7e-5p/hsa-miR-98-5p were used for further screening. Cytoscape was run again to map the comprehensive regulatory network composed of interconnected ceRNAs. Finally, 4/2/4 lncRNAs upstream of hsa-miR-365a-3p/hsa-let-7e-5p/hsa-miR-98-5p were identified from the following three analysis steps: 1) The correlation analysis between 13/13/11 lncRNAs and hsa-miR-365a-3p/hsa-let-7e-5p/hsa-miR-98-5p expression (filter condition: coefficient < −0.28 and p < 0.001); 2) The correlation analysis between 13/13/11 lncRNAs and BMP2 expression (filter condition: correlation coefficient < −0.4 and p < 0.001); 3) The differential analysis of 13/13/11 lncRNAs expression between glioma and normal tissues (filter condition: |log2FC| >1 and p < 0.05). Kaplan-Meier survival analysis was again used to explore the influence of these lncRNAs on prognosis. We highlighted these lncRNAs with yellow in the comprehensive regulatory network.

To further predict the regulatory network constructed by circRNAs, miRNAs, and BMP2, we also downloaded the circRNAs expression data of the GSE165926 set (12 gliomas and 4 normal brain tissues) and the miRNAs expression data of the GSE138764 set (33 astrocytoma and 9 normal brain tissues). Based on the TCGA and GSE138764 sets, respectively, the differentially expressed miRNAs (DE-miRNAs) between glioma and normal brain tissues were obtained (filter condition: |log2FC| >2 and p < 0.05). Prognosis related differentially expressed miRNAs (PR-DE-miRNAs) were obtained by Kaplan-Meier survival analysis based on TCGA data (filter condition: p < 0.001). We re-identified the common miRNAs (hsa-miR-129-5p and hsa-miR-381-3p) from the PR-DE-miRNAs of the TCGA set, the DE-miRNAs of the GSE138764 set, and the miRNAs upstream of BMP2 predicted by starBase as candidate miRNAs that may regulate BMP2 expression. StarBase was again used to predict circRNAs that might bind to hsa-miR-129-5p and hsa-miR-381-3p, respectively. A similar approach was used to obtain differentially expressed circRNAs (DE-circRNAs) from the GSE165926 set. We got the common circRNAs from DE-circRNAs of the GSE165926 set and the circRNAs predicted by starBase as candidate circRNAs that might bind to hsa-miR-129-5p/hsa-miR-381-3p, respectively. Finally, only the up-regulated hsa_circ_0004662 and hsa_circ_0007548 in gliomas upstream of hsa-miR-129-5p satisfy the regulation of ceRNA. Cytoscape plots the corresponding regulatory network.

We obtained the immunophenotypic score (IPS) of glioma patient from the Cancer Immunology Atlas (TCIA) database (tcia.at/home). Previous studies showed that immunogenicity increases with increasing IPS score (Liu et al., 2020b). By analyzing the gene expression of the four cell types that determine immunogenicity (effector cells, immunosuppressive cells, major histocompatibility complex molecules and immunomodulators), the IPS of the sample was obtained (Liu et al., 2020b). Spearman correlation analysis was used to analyze the correlation between 4 types of IPSs and 5 PR-DE-IRGs expression/risk score. We used the violin plot to show the difference of Tumor Immune Dysfunction and Exclusion (TIDE), microsatellite instability (MSI), T cell dysfunction and exclusion score between high and low risk groups. TIDE scores were calculated online at the TIDE website (tide.dfci.harvard.edu/).

We accessed the NCI-60 database containing 60 different cancer cell lines from 9 different types of tumors through the CellMiner interface (discover.nci.nih.gov/cellminer). The Pearson correlation analysis was run to explore the relationship between 5 PR-DE-IRGs expression and 263 drugs approved by the food and drug administration or clinical trials.

We used the R package “regplot” to combine the risk score with two clinical prognostic factors (grade, age) to establish a combined nomogram, which can more accurately predict the survival probability of glioma patients at 1-, 2-, and 3- years. To better prove the accuracy and effectiveness of nomogram, we constructed internal calibration curves and ROC curves based on the training, test, and whole sets.

As the most common covalent modification of RNA at the posttranscriptional level, N6-methyladenine (m6A) mRNA modification plays a key role in gliomas through various mechanisms and associates with clinicopathological features and prognosis of gliomas patients, showing a great clinical significance (Zhang et al., 2020; Zhang et al., 2021). To explore the correlation between risk score and Immune Checkpoint Inhibitors (ICIs)/m6A related gene expression, we again used Spearman correlation analysis. In addition, we also compared the expression differences of ICIs/m6A related genes between high and low risk groups. Studies have shown that the multidrug resistance-associated protein (MRP or ABCC) expression may be related to the inherent multidrug resistance in human gliomas (Mohri et al., 2000), so we respectively analyzed the correlation between the ABCC1/ABCC3 expression and the risk score, as well as their expression differences between different risk groups.

We used a heatmap to show the distribution of clinical characteristics across all TCGA samples from different risk groups. To visualize the differences in risk score between different subgroups for each clinical characteristic, we drew boxplots. Kaplan-Meier survival analysis was used to evaluate the model’s predictive ability in each subgroup with different clinical characteristics. We obtained the immunohistochemical (IHC) staining results from the human protein atlas database (proteinatlas.org/) to verify the differences in protein levels of 5 modelled genes between normal tissues and glioma tissues. Finally, only 4 PR-DE-IRGs (HMGB2, CCNA2, CASP3, and APOBEC3C) IHC staining images were obtained.

We used Student’s t-test or Wilcoxon signed-rank test to compare the differences between continuous variables, while the differences between categorical variables were compared by Chi-square test or Fisher’s exact test (Fan et al., 2021). Univariate cox regression analysis was used to select PR-DE-IRGs. Lasso and multiple Cox regression were used to construct prognostic model and nomogram. Univariate and multivariate Cox regression analyses were also used to identify independent prognostic factors. Spearman or Pearson correlation analyses were used to analyze the correlation among variables. The above analyses were run in the R software (version 4.0.3), Perl and SPSS Statistics 22.

We extracted 455 DE-IRGs (189 genes were down-regulated and 266 genes were up-regulated in glioma) from the TCGA cohort (Supplementary Figure S1) and 968 DE-IRGs (492 genes were up-regulated and 476 genes were down-regulated in glioma) from the GSE108474 cohort (Supplementary Figure S2). In the same way, we also identified 986 DE-IRGs from the GSE4290 cohort (Supplementary Figure S3). Finally, by Venn diagram based on the 6 cohorts, we obtained 163 common DE-IRGs (Figure 2A).

The merged dynamic pruning tree in Figures 3A,C,E showed the results of genes and corresponding modules matched in the TCGA, GSE4290 and GSE108474 cohorts, respectively. Then we identified 3 modules from the TCGA cohort (Figure 3B), 2 modules from the GSE4290 cohort (Figure 3D), and 2 modules from the GSE108474 cohort (Figure 3F). The heatmaps (Figures 3B,D,F), reflecting the correlation between modules and tumor status, showed that the turquoise module from the GSE4290 cohort had the strongest correlation with tumor status (|correlation coefficient| = 0.64, p < 0.001, Figure 3D). Finally, we got 41 common DE-IRGs highly associated with glioma using the Venn diagram (Figure 2B).

668 samples from the TCGA cohort, 85 samples from the GSE4412 cohort, 50 samples from the GSE43378 cohort, and 983 samples from the CGGA cohort were used in the subsequent analyses. Supplementary Tables S1,2 showed the statistical results of these samples’ clinical characteristics. During the overlapping process of TCGA cohort’s 41 PR-DE-IRGs, GSE4412 cohort’s 17 PR-DE-IRGs, GSE43378 cohort’s 25 PR-DE-IRGs, and CGGA cohort’s 40 PR-DE-IRGs, we obtained 15 common PR-DE-IRGs (Supplementary Figure S4A). The forest plots showed the univariate Cox analyses’ results of 15 common PR-DE-IRGs from 4 cohorts (TCGA cohort: Figure 2D; CGGA cohort: Supplementary Figure S4B; GSE4412 cohort: Supplementary Figure S4C; GSE43378 Cohort: Supplementary Figure S4D). Except for BMP2 and PRKX (HR < 1), the remaining 13 genes (HR > 1) were identified as prognostic risk factors in all cohorts. This conclusion was supported by further Kaplan-Meier survival analysis results (TCGA cohort: Figure 4; CGGA cohort: Supplementary Figure S5; GSE4412 cohort: Supplementary Figure S6; GSE43378 Cohort: Supplementary Figure S7). In addition, we also used a heatmap to visualize the expression levels of these 15 PR-DE-IRGs in the TCGA cohort (Figure 2C).

Finally, lasso regression analysis determined 5 PR-DE-IRGs and corresponding coefficients (Supplementary Table S3) to build the model according to the optimal penalty parameter (λ) (Supplementary Figures S8A,B).

With survival point plots, we observed more deceased samples in the high-risk group (Figures 5A–C and Figures 6A–C). Kaplan-Meier survival curves further showed the lower survival probability for high-risk group samples (Figures 5D–F and Figures 6D–F). To better verify the model’s predictive performance, we drew the ROC curves. In the ROC curves of three TCGA’s sets, all area under curve (AUC) values were greater than 0.85 (Figures 5G–I), showing our model’s excellent performance. This conclusion was also supported by the ROC curves of the three external cohorts (all AUC >0.7) (Figures 6G–I). The risk score was identified as a prognostic risk factor by the univariate Cox regression analysis (Figures 5J–L and Figures 6J–L). After adjusting for various clinical factors, the multivariate Cox regression analysis further determined the risk score as an independent prognostic risk factor (Figures 5M–O and Figures 6M–O). Finally, both PCA and t-SNE dot plots showed that the risk grouping separated the samples well (Supplementary Figure S9).

In the low-risk group, we found GO analysis results (neurotransmitter secretion, neurotransmitter transport, regulation of neurotransmitter levels, regulation of postsynaptic membrane potential, and regulation of synaptic plasticity) (Figure 7B) and KEGG analysis results (calcium signaling pathway, cardiac muscle contraction, long-term potentiation, and neuroactive ligand-receptor interaction) (Figure 7D) were all closely related to the nervous system. In the high risk group, both GO analysis results (CD8-positive alpha-beta T cell activation, interferon gamma mediated signaling pathway, interferon gamma production, positive regulation of B cell activation, regulation of T cell activation, T cell activation involved in immune response, T cell mediated immunity and toll-like receptor 2 signaling pathway) (Figure 7A) and KEGG analysis results (cytokine-cytokine receptor interaction, JAK-STAT signaling pathway, natural killer cell mediated cytotoxicity, NOD-like receptor signaling pathway, p53 signaling pathway and toll-like receptor signaling pathway) (Figure 7C) contained numerous immune-related biological processes. Because GSEA exhibited a strong link between high-risk group and immunity, we further analyzed the relationship between the model and the immune microenvironment. Figures 7E,F showed a positive correlation between risk score and overall immune cell score, and higher overall immune cell score in the high-risk group. Similar results were observed in the overall stromal cell score (Figures 7G,H). After visualizing 16 immune cells and 13 immune functions scores for each sample (Figure 7I), we further analyzed their correlation with risk score, and their differences between high-risk and low-risk groups. The results showed that the remaining immune cells and functions scores were positively correlated with risk score except for natural killer (NK) and mast cells (Figure 7J). The above results were consistent with the differential analysis results presented in Figure 7K.

Figure 8A showed the distribution of C3, C4, and C5 between the high-risk and low-risk groups. We observed the most C5 subtypes in the low-risk group and the most C4 subtypes in the high-risk group. The results in Figure 8B further verified the above situation. Cell stemness analysis was of great significance to the epigenetic characteristics of glioma patients. We observed a significant correlation between the mDNAsi and risk score (Figure 8C), and higher mDNAsi in the high-risk group (Figure 8D). Finally, we observed higher mDNAsi in the age >60 group (Figure 8E), deceased group (Figure 8F) and G3 group (Figure 8H). But, there was no significant difference in mDNAsi between different gender subgroups (Figure 8G).

From Figure 9A, we found that 32 of the 41 DE-IRGs highly associated with glioma were mutated. The mutation frequency of IDH1 was the highest in different risk groups (Figures 9B,C). We observed that TMB was positively correlated with risk score, which was supported by further differential analysis (Figures 9D,E). The survival analysis in Figure 9F showed that the high TMB group samples have significantly lower survival probabilities (p < 0.001). Lower risk score and higher survival probability were observed in the IDH1 mutant group (Figures 9G,H). Except for BMP2, the expression levels of other genes (CD274, APOBEC3C, CASP3, CCNA2 and HMGB) in the wild-type group were higher (Figures 9I,K–O). Figure 9J showed that NK cells infiltrated more abundantly in the mutant group, while the others infiltrated more abundantly in the wild-type group.

After reasonable prediction, we obtained 52 upstream miRNAs that may bind to BMP2. Figure 10J showed the regulatory network between these 52 miRNAs and BMP2. Supplementary Table S4 showed all predicted miRNAs’ correlation and differential analysis results. Figures 10A–I visualized the results of the correlation, difference and Kaplan-Meier survival analyses for 3 miRNAs (hsa-miR-365a-3p, hsa-let-7e-5p, and hsa-miR-98-5p). Therefore, these three miRNAs may become the most potential miRNAs regulating BMP2 expression in glioma.

Based on the starBase database, we successfully obtained 112/101/112 possible upstream lncRNAs of miR-365a-3p/let-7e-5p/miR-98-5p, respectively. Figures 10K–M showed the comprehensive regulatory network composed of interconnected ceRNAs (13 lncRNAs, hsa-miR-365a-3p and BMP2; 13 lncRNAs, hsa-let-7e-5p and BMP2; 11 lncRNAs, hsa-miR-98- 5p and BMP2). Then, 4 lncRNAs (TUG1, LINC00689, NNT-AS1 and ZEB1-AS1), 2 lncRNAs (SNHG16 and AL590666.2), and 4 lncRNAs (PXN-AS1, ZEB1-AS1, TUG1 And NNT-AS1) were screened out separately through the analysis of next three steps. Supplementary Figures S10A–J showed the results of the correlation and difference analyses for these 10 lncRNAs. Supplementary Figure S10K showed the Kaplan-Meier survival curves of these seven lncRNAs.

We obtained 2 common miRNAs from 59 PR-DE-miRNAs of TCGA set, 92 DE-miRNAs of GSE138764 set, and 52 miRNAs of starBase (Supplementary Figure S11A). Supplementary Figures S11B–E showed the results of difference and survival analyses corresponding to these two miRNAs. 11/1 common circRNAs were obtained from 269 DE-circRNAs of the GSE165926 set and 6618/3647 circRNAs upstream of hsa-miR-129-5p/hsa-miR-381-3p predicted by starBase, respectively (Supplementary Figures 11F,G). Supplementary Figure S11H visualizes the differential expression results of these 12 circRNAs. The mechanisms by which miRNAs regulate target gene expression suggest that there should be negative correlations between miRNAs and BMP2, and between miRNAs and circRNAs. Only hsa_circ_0004662 and hsa_circ_0007548, which were up-regulated in gliomas, satisfied this regulation. Finally, the regulatory network composed of hsa_circ_0004662/hsa_circ_0007548, hsa-miR-129-5p and BMP2 was shown in Supplementary Figure S11I.

IPS is used to assess which patients are more inclined to respond to ICIs therapy. A higher IPS evaluation value implies a better therapeutic effect of ICIs (Liu et al., 2020b). IPS showed broad correlation with risk score, CASP3, APOBEC3C, HMGB2, and BMP2 expression (Figures 11A,B). Such results suggested that BMP2, CASP3, APOBEC3C, HMGB2 and risk score can be used to predict the efficacy of ICIs. TIDE score was negatively correlated with the therapeutic effect of ICIs (Liu et al., 2020b). Figures 11C,F showed that there was no significant difference in TIDE score and T cell dysfunction score between the high-risk and low-risk groups. But, the T cell exclusion and MSI scores in the low-risk group were higher (Figures 11D,E).

The correlation analysis results of the 5 drugs most closely related to 5 PR-DE-IRGs were visualized, respectively (Figures 12A–E). These results implied the guiding value of the five PR-DE-IRGs in these drugs’ clinical application. Unexpectedly, the expression of CASP3 was positively correlated with the sensitivity of carmustine, a drug recommended by the latest National Comprehensive Cancer Network (NCCN) for the treatment of glioma (Figure 12A).

Finally, we built a nomogram based on age, grade and risk score to predict the glioma patients’ survival probability (Figure 13A). The internal calibration curves showed that the predicted results were basically consistent with the actual results (Figures 13B–D). The ROC curves further confirmed the excellent predictive performance of the nomogram (AUC >0.8) (Figures 13E–G).

Except for VTCN1, CD200, TNFSF9, HHLA2, and ADORA2A, which were negatively correlated with the risk score, the expressions of the other ICIs-related genes were all positively correlated with the risk score (Figure 14A). We also observed that half of the m6A-related genes’ expression were positively correlated with risk score, while the other half were negatively associated with the risk score (Figure 14A). The results of the above ICIs-related and m6A-related genes were verified in the further difference analysis (Figures 14B,C). We observed a positive correlation between ABCC1/ABCC3 expression and risk score as well as their higher expression in the high-risk group, further supporting the potential value of the model in drug resistance prediction (Figures 14D–G).

The heat map showed the distribution of clinical characteristics across all TCGA samples from different risk groups (Figure 15A). Subsequently, we further explored the relationship between survival status (Figure 15B)/age (Figure 15C)/gender (Figure 15D)/grade (Figure 15E) and risk score. The risk score of patients in the deceased group, age >60 group and G3 group were higher (Figures 15B,C,E). Figure 15F showed that the model still maintained the excellent performance of discriminating prognosis in all clinical characteristics subgroups. Figures 15G–J showed the IHC staining images reflecting the protein expression levels of the four PR-DE-IRGs (HMGB2, CCNA2, CASP3, and APOBEC3C) in the model, respectively. We only observed higher expression of three genes (HMGB2, CCNA2, and CASP3) in the glioma tissues (Figures 15G–I), consistent with the conclusions from the previous analysis based on the TCGA cohort. The above results further verified the model’s stability.

Figures 15G–J showed the IHC staining images reflecting the protein expression levels of the four PR-DE-IRGs (HMGB2, CCNA2, CASP3 and APOBEC3C) in the model, respectively. We only observed higher expression of three genes (HMGB2, CCNA2 and CASP3) in the glioma tissues (Figures 15G–I), consistent with the conclusions from the previous analysis based on the TCGA cohort. The above results further verified the model’s stability.