We searched the available mRNA expression profiles with prognosis information from several tumor resources, such as the Gene Expression Omnibus (GEO) database, Chinese Glioma Genome Atlas (CGGA) database, and Pan-Cancer Analysis of Whole Genomes (PCAWG) database for glioma prognostic microenvironment identification. For the datasets from the GEO database, the cohorts with at least 40 samples were considered and a total of six public glioma cohorts were downloaded. We also collected two glioma cohorts (LGG and GBM) from TCGA database and the multi-omics data including gene expression, methylation level, and copy number variations (CNVs) from the Human Glioma Cell Culture (HGCC) collection as the independent testing set. Finally, 6,920 brain samples from 27 cohorts were included in our study, and the total information was shown in Table 1 .

We collected 175 immune microenvironment-related signatures from diverse literature for glioma prognostic signature identification. In detail, 28 signatures were obtained from the work of Bindea et al. (10), 11 signatures were obtained from the work of Wolf et al. (11), 24 cell signatures were obtained from the work of Miao et al. (12), 22 immune signatures were obtained from CIBERSORT (13), 29 immune microenvironment signatures were obtained from the work of Bagaev et al. (14), 40 signatures were obtained from Cellmarker database (15), 4 microglia signatures were obtained from the work of Sala Frigerio et al. (16), 8 brain immune cell signatures were obtained from scREAD database (17), and 9 microglia subtype signatures were obtained from the work of Olah et al. (18) More detailed information is listed in Table S1 .

The datasets from GEO, CGGA, and PCAWG databases were treated as the training set, and the datasets from TCGA were treated as the testing set. Based on the 175 tumor microenvironment signatures, we firstly utilized the single sample Gene Set Enrichment Analysis (ssGSEA) method implemented in the R package to calculate the normalized enrichment score (NES) for each glioma sample from the training set (19). To remove the potential batch effects from different cohorts, we further used Combat function to form a merged NES matrix for microenvironment signature identification. Finally, the prognostic performance of this signature was tested using testing datasets.

Based on the merged NES matrix, which consisted of 175 immune microenvironment signatures and 4,887 glioma samples, we firstly performed univariable Cox proportional hazards regression analysis using 3,486 samples with survival analysis. A set of 141 signatures was identified with a prognostic P-value <0.05. The detailed univariable Cox results, including the hazard ratio (HR) value, 95% CI, and P-values, of these microenvironment signatures were provided in Table S2 . The LASSO Cox regression model (20) was used to find the most useful prognostic markers among the 141 immune microenvironment signatures in the training cohort by the R package glmnet. Ultimately, 25 immune microenvironment features with non-zero coefficients were selected through LASSO Cox regression model analysis, and the optimal lambda value was determined by 10-fold cross-validation. The multivariable Cox proportional hazards regression analysis was further performed for these 25 immune microenvironment signatures obtained from the LASSO analysis, and the corresponding coefficients are presented in Table S3 . Finally, a novel glioma prognostic signature, GMEFS, was constructed by comprehensively considering the coefficients of these 25 immune signatures. The formula was provided as follows:

where HR_i is the HR and NES_i is the NES for the ith immune microenvironment signature.

The ESTIMATE score, tumor purity, and stromal and immune scores for each glioma sample were calculated by using ESTIMATE package in R with default parameters (21). The differential expression (DE) analyses for TCGA GBM and LGG cohorts were performed by the limma package in R (22). Functional enrichment analyses, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and subpathway analyses, were performed by using the clusterProfiler package in R (23), and the P-value was adjusted by the Benjamini and Hochberg method.

Based on the HGCC resource (24), we obtained the drug IC50 information, as well as gene expression, methylation, and CNV data for each GBM cell line. Firstly, based on the median IC50 value as cutoff, we defined two cell line groups, high IC50 groups and low IC50 groups. Then, based on these two groups, we respectively identified drug-related genes according to gene expression level, methylation condition, and CNV data. For gene expression profiles, the T-test was used. For methylation and CNV data, the Wilcoxon rank sum test was used. The cutoff for DE analysis was set as adjusted P-value <0.05. Finally, we evaluated the associations between DE genes and GMEFS subpathways by using the hypergeometric test method. The result with P-value <0.05 was considered as a significant association. An integrated drug–subpathway network was constructed by considering the significant drug–subpathway associations shared by two omics results.

According to the GMEFS score, all glioma samples from both training and testing sets were classified into two groups based on the consistent cutoff as the training set. Then, the Kaplan–Meier (KM) curve and survival P-value calculated by the log-rank test were performed by using R survminer package. The glioma molecular subtypes (Classical-like, Codel, G-CIMP-high, G-CIMP-low, Mesenchymal-like, PA-like) of the glioma patients were obtained from a previous study (25). For both training and testing sets, a nomogram was formulated to provide a visualized risk prediction after each factor was assigned a score. A calibration plot was generated to assess the calibration ability of the nomogram (26). Nomograms and calibration plots were generated by using the rms package. The decision curve analysis (DCA) was performed by using dca package. All of the P-value results were considered statistically significant with P-values <0.05.

Based on the training set that consisted of a total of 8 glioma cohorts, we constructed a microenvironment-based prognostic model named GMEFS (see Materials and Methods). The GMEFS consisted of 25 glioma immune microenvironment signatures, and the detailed LASSO results of these 25 signatures were displayed in Figures 1A, B . We calculated the GMEFS score for each glioma patient in the training cohort and stratified the patients into high or low GMEFS groups according to the median cutoff. As shown in Figure 1C , the samples with a high GMEFS score had a significantly shorter overall survival than samples with a low GMEFS score in the training cohort (P < 0.0001; log-rank test). The associations between each of the 25 immune microenvironment signatures and clinical survival are also shown in Figure 1D . To further examine the robustness of the GMEFS model, two independent glioma cohorts, GBM and LGG datasets, were also obtained as testing sets for the prognosis analysis (see Materials and Methods). With the same formula, the samples from the testing set were stratified into high and low GMEFS groups by the cutoff value obtained from the entire training set. As shown in Figures 1E, F , the samples with a high GMEFS score also had significantly worse overall survival than those who displayed a low GMEFS score in both GBM and LGG cohorts (GBM P = 0.0091, LGG P < 0.0001; log-rank test). Meanwhile, survival analysis for progression-free survival (PFS) was also performed, and the predictive performance of GMEFS in GBM ( Supplementary Figure 1A , P = 0.015) and LGG ( Sepplementary Figure 1B , P < 0.0001) was confirmed. Similar results were observed in the entire TCGA cohort for both overall survival and PFS ( Supplementary Figures 1C, D ).

To further test the performance of predicting the glioma patient prognosis, a nomogram that integrated both the GMEFS and clinical factors (including gender and age) was constructed by using patients from the training set. Based on the nomogram results, a score can be calculated for a glioma patient for predicting the 3-, 5-, and 10-year overall survival for an individual, suggesting the power of GMEFS score in contributing the risk point ( Figure 2A ). The calibration curves for 3, 5, and 10 years of the training cohort were respectively illustrated in Figures 2B–D , showing the GMEFS’s performance in the glioma patient survival prediction.

To test the associations between GMEFS and molecular subtypes, we obtained the glioma subtype information (see Materials and Methods). As shown in Figure 3A , the samples from Classic-like, G-CIMP-low, and Mesenchymal-like subtypes displayed high GMEFS scores, whereas samples from the Codel subtype had a low GMEFS score, and the overall difference was significant. Moreover, an imbalance in terms of these five molecular subtypes within the two GMEFS group was observed ( Figure 3B ). Only small proportions of Codel subtype (6.5%) distributed in the high GMEFS group in contrast with 100% of Mesenchymal-like type and 100% of Classic-like type. Similar results were observed between molecular subtype and GBM GMEFS groups. Then, the correlations between tumor purity, stromal characterization, and GMEFS score were also explored. It was observed that the GMEFS score was positively correlated with the stromal score, immune score, and ESTIMATE score, whereas it was negatively correlated with tumor purity for TCGA LGG cohorts ( Figure 3C ), showing that the GMEFS score was a presentation of the immune microenvironment characterization.

Some potential factors that determine brain tumor immunogenicity, including macrophage regulation, homologous recombination deficient (HRD), TGF-beta response, and wound healing, were compared between the high GMEFS group and the low GMEFS group from TCGA LGG cohort ( Figure 3D ). The mean activity values for other factors of the high GMEFS group were significantly higher than those of the low GMEFS group. In addition, the samples with a higher GMEFS displayed a higher proliferation score than that of patients with a lower GMEFS, which was consistent with the prognostic performance of GMEFS. Overall, the differences in brain tumor immunogenicity between the GMEFS groups were significant, and the high GMEFS group had a relatively high immunogenicity. Next, we identified 83 upregulated genes in the high GMEFS group vs. the low GMEFS group by the cutoffs of log-fold change (FC) >1.5 and false discovery rate (FDR) <0.05. The GO biological process terms enriched by these genes included many immune-related processes, such as neutrophil activation, antigen processing and presentation, and neutrophil-mediated immunity ( Figure 3E ). Similar GO biological process terms and KEGG pathway results were illustrated in TCGA-GBM cohort ( Supplementary Figure 2 ).

Furthermore, we quantified these 83 LGG gene signatures in data from the Ivy Glioblastoma Atlas Project (IGAP), which performed RNA sequencing (RNA-seq) on microdissections of glioma anatomical structures from hematoxylin and eosin (H&E) staining (Ivy Glioblastoma Atlas Project; http://glioblastoma.alleninstitute.org). The higher expression activity of these genes was enriched in samples from the hyperplastic blood vessels, microvascular proliferation, and perinecrotic zone. These gene signatures displayed a lower expression level in infiltrating tumor and leading edge ( Figure 3F ). To test the glioma association of 83 genes, we further obtained several glioma-related gene sets from a previous study (27). As shown in Supplementary Figure 3A , these genes displayed a significant overlap with GBM disease genes. These GBM genes were also enriched in the samples with a high GMEFS score compared to those with a low GMEFS score ( Supplementary Figure 3B ).

To further explore the different immune cell proportions within the two GMEFS groups, we performed CIBERSORT analysis with 1,000 permutations for glioma cohorts (13). As shown in Figures 4A, B , the proportions of 22 common immune cell types of the two GMEFS groups of LGG were displayed. The samples from the low GMEFS group had significantly higher proportions of plasma cells, monocytes, activated mast cells, and eosinophils than those of samples from the high GMEFS group (P < 0.05). Correspondingly, the proportions of memory resting CD4 T cells, M2 macrophages, and CD8 T cells in the high GMEFS group were significantly higher than those in the low GMEFS group (P < 0.05). We further studied the correlation between the GMEFS score and the proportion of immune cell infiltration. The results revealed that CD8 T cells and M2 macrophages were positively correlated with the GMEFS, and activated mast cells and T helper cells were negatively correlated with GMEFS ( Figure 4C ). By comparing brain cell-type marker gene from the PanglaoDB database (28), we identified 9 cell types associated with the two GMEFS groups, including microglia, astrocytes, and oligodendrocyte progenitor cells ( Figure 4D ). Parenchymal microglia are the major component of myeloid cells involved in heterogeneous Central Neural System (CNS) immune microenvironment (29, 30), and microglia play important roles in the brain-related tumors, including glioma (31).

To test the performance of the GMEFS in glioma formation and recurrence, we further obtained the glioma datasets from the GEO database ( Table 1 ). The GMEFS scores were calculated for each analytical sample, and the difference was evaluated between glioma and normal samples, as well as recurrence glioma and primary glioma samples. As shown in Figure 5A , among 14 independent glioma cohorts, the glioma samples displayed consistently higher GMEFS scores than those in normal samples in a total of 12 cohorts. These significant cohorts contained glioma samples (11/12) and Diffuse Intrinsic Pontine Glioma (DIPG) samples (1/12). In addition, within 2 of 7 glioma cohorts, the recurrence samples also displayed higher GMEFS scores than those in primary samples ( Figure 5B ). All these findings show that the GMEFS was not only a risk characterization for evaluating high-risk samples from low-risk samples but also could distinguish glioma samples from normal samples, considering the diagnosis index for glioma patients. As shown in Figures 5C, D , within two glioma cohorts with available stage information, we observed that the samples with a higher stage displayed a higher GMEFS score than samples with a lower stage, which was consistent with the prognostic performance of the GMEFS groups (GSE16011, P = 1.9e-08; GSE4290, P < 2.2e-16). Using these testing sets from the GEO, we further confirmed the GMEFS difference in histology subtypes and molecular subtypes, which was also consistent with previous results (GSE60898 and GSE72951 from Figures 5E, F , Figure 3A ). As shown in GSE116520 ( Figure 5G ), the GMEFS score was also related to glioma region distribution; GMEFS was higher in core tumor tissue and lower in non-neoplastic brain tissue.

In our previous studies (32, 33), we observed that the regions of the whole pathway, which is also named as subpathway, was closely related to disease formation and progression. The related results further confirmed that the subpathway displayed more advantages over the whole pathway with providing detailed biological information. We firstly obtained the subpathway list from subpathwayMiner R package with the default parameters (33). Based on the training, TCGA-GBM, and TCGA-LGG cohorts, we utilized the GSEA algorithm to identify the significantly enriched subpathways between the high GMEFS and low GMEFS groups. As shown in Figures 6A, B , a total of 32 risk subpathways were shared by three glioma cohorts. These risk subpathways were derived from Focal adhesion, ECM–receptor interaction, and pathways in cancer. Notably, a total of 8 subpathways derived from Focal adhesion were identified. Take one subpathway (path:04510_9) as an example, and the subpathway activity was calculated for the corresponding tumor samples. As shown in Figures 6C–E , the Focal adhesion subpathway could distinguish all samples into two groups with significant prognosis (Training cohort, p < 0.0001; TCGA-GBM, p = 0.017; TCGA-LGG, p < 0.0001). To explore the detailed associations between antineoplastic compounds and GMEFS subpathways shared by two cohorts, we constructed a multi-omic integrated network based on HGCC resource (see Materials and Methods). As shown in Figure 6F , many subpathways were targeted by many candidate drugs, such as path: 04062_2 from chemokine signaling pathway and path: 04510_9 from focal adhesion. The subpathway 04062_2 was targeted by many approved molecules, such as mesalamine, pamelor, angormin, and amlodipine. Some subpathways from the same whole pathway (path: 04510) were commonly targeted by some molecules, such as etoposide, paludrine, and suloctidil. Notably, as an experimental drug, suloctidil was closely related to many risk subpathways. Furthermore, we obtained several drug target genes for glioma collected from a recent study (34). We observed that many subpathways within networks contained known drug targets. Notably, three subpathways (path: 04510_8, _10, and _12) from Focal adhesion contained four target genes, including PDGFRA, EGFR, and SRC. The systematically reconstructed drug–subpathway network provided a novel framework for identifying subpathway regions as candidate therapy targets and screening-approved or experimental molecules as novel drugs for clinical use.