The bulk mRNA sequencing data of GBM with corresponding demographics were included. Of these, 166 cases retrieved from The Cancer Genome Atlas (TCGA) database (https://www.cancer.gov/) were used as the training dataset, and the remaining cases retrieved from the Chinese Glioma Genome Atlas (CGGA) database (mRNA-seq 325, mRNA-seq 693, and microarray 301, http://www.cgga.org.cn/) were used as the validation dataset. The mRNA-seq data were TPM normalized, and the microarray dataset was log-transformed. A total of 2,498 IRGs were retrieved from the ImmPort Portal (https://www.immport.org/). The somatic mutation profiles and copy number variation (CNV) files were retrieved from the TCGA database.

An unsupervised class discovery technique was employed to find stable clusters based on the top 300 IRGs ranked by their mean absolute deviation (MAD) values using the R package “ConsensusClusterPlus” (Supplementary Material S1) (Wilkerson and Hayes, 2010). The clustering method was set as “km” (k-means), and the maximum evaluated K was 6. The optimal number of clusters was determined using the “proportion of ambiguous clustering (PAC)” method.

The weighted gene coexpression network analysis (WGCNA) was performed to identify genes characterizing an active immune response of GBM based on the R package “WGCNA” (Langfelder and Horvath, 2008). Module Eigenges (ME) was defined as the first principal component of each gene module and served as the representative for all genes in each module. Gene significance (GS) represents the degree of linear correlation between module gene expression and sample features. Genes with weighted correlation coefficients over 0.8 were extracted for further analysis.

The univariate and multivariate Cox regression analyses were performed sequentially to identify genes with prognostic significance as well as to determine their regression coefficients. Genes with p-values <0.05 in the univariate regression analysis were selected for multivariate regression analysis. The coefficients obtained from the regression analysis were used to yield the following risk score equation: risk score = ∑ k = 1 n regression coefficient × gene expression value

Samples were split into the high- and low-risk groups by the median value of the risk score. The Kaplan–Meier (K-M) plots were generated based on the R package “survival”, and the independent prognostic value of the risk score was evaluated by univariate Cox regression analysis. The time-dependent predictive power of the risk model was evaluated using the receiver operative curve (ROC) and the corresponding area under the curve (AUC). The independent prognostic value of the risk group, as well as clinicopathological parameters, was evaluated using the univariate Cox regression analysis.

The differential gene expression profile was conducted between subgroups using the R packages “limma” and “edgeR” (Robinson et al., 2010; Ritchie et al., 2015). The absolute value of logFC >1.0 and adjust p-value <0.05 were set as the cutoff. Functional enrichment analysis was performed based on the R package “clusterProfiler” and the online tool The Database for Annotation, Visualization and Integrated Discovery (DAVID, version 6.8, https://david.ncifcrf.gov/) (Huang et al., 2009a; Huang et al., 2009b). The Gene Ontology (GO) terms including biological process (BP), cell component (CC), molecular function (MF), and Kyoto Encyclopedia of Genes and Genomes (KEGG) and Biocarta pathways were involved. The false discovery rate (FDR) < 0.05 was set as the cutoff. Gene Set Enrichment Analysis (GSEA) was employed to identify differentially activated signaling pathways (Subramanian et al., 2005). Gene sets were permuted 1,000 times. Normalized enrichment scores (NES) were calculated, and FDR < 0.05 was set as the cutoff. To explore immunological characteristics, 103 gene sets involved in the inflammation/innate immune response, antigen presentation, CD8 T cell function, and cytotoxicity were retrieved from the Molecular Signatures Database (MSigDB, v7.4, http://software. broadinstitute.org/gsea/msigdb/index.jsp) (Supplementary Material S1), and the ssGSEA score of the samples was calculated to assess the activation of specific signaling pathways (Liberzon et al., 2011; Hänzelmann et al., 2013). The fraction of immune cells was estimated using the CIBERSORT algorithm (Newman et al., 2015). Samples with p-values over 0.05 were filtered. Tracking tumor immunophenotype (TIP) is a system to quantify the activity of each step of the anti-tumor immune response (Xu et al., 2018). The anti-tumor immune response was conceptualized as seven stepwise events including 1) release of cancer cell antigens, 2) cancer antigen presentation, 3) priming and activation, 4) trafficking of immune cells to tumor, 5) infiltration of immune cells into tumor, 6) recognition of cancer cell by T cells, and 7) killing of cancer cells. The TIP score of each stepwise event was evaluated using the R script of TIP.

Somatic mutation data were analyzed using the R package “maftools” (Mayakonda et al., 2018). The differentially mutated genes were calculated using the functions “mafComapre.” The tumor mutational burden (TMB) was defined as the total number of somatic mutations including common substitutions, insertions, and deletions per megabase, and its calculation has been described before (Chalmers et al., 2017). Significant amplifications and deletions of somatic copy numbers were detected using GISTIC 2.0 (Mermel et al., 2011). FDR-q value ≤ 0.05 was set as the cutoff.

The tumor responsiveness to the ICB therapy was assessed using the machine learning algorithm TIDE, which is a computational method to model the induction of T cell dysfunction and the prevention of T cell infiltration into tumor through analyzing the expression of specific gene signatures in the tumor expression profile (Jiang et al., 2018). Further, the predicted tumor response to ICB was evaluated by an unsupervised subclass mapping method that is designed to determine correspondence or commonality of subtypes found in multiple independent datasets generated on different platforms (Hoshida et al., 2007).

All statistical analysis was conducted using R software (version 4.0.2). The differences in gene expression, ssGSEA score, TIP score, and the fraction of immune infiltration were assessed using a two-tailed Wilcoxon's test. Univariate and multivariate Cox regression analyses were conducted to establish Cox proportional hazard models. Differences in clinicopathological parameters between the high- and low-risk groups were measured using the Wilcoxon test and Fisher's exact test. The K-M analysis and log-rank analyses were employed to evaluate differences in OS. The independent prognostic value of the risk score was assessed using univariate Cox analysis. Pearson correlation analyses were conducted to estimate the correlation between the risk score and the fraction of immune infiltration. The association between epidermal growth factor receptor (EGFR) mutations and immune infiltration was assessed using logistic regression analysis. p < 0.05 was considered statistically significant. We marked p < 0.05 as *, p < 0.01 as **, and p < 0.005 as ***.

To obtain GBM subgroups with potential distinct immunological properties, samples were stratified according to the expression of IRGs. The MAD value of 1,288 IRGs with an intersection with the TCGA profile was calculated; thereafter, 166 samples with survival information were stratified based on the expression of top 300 IRGs in terms of MAD value, and two subgroups [termed group 1 (n = 87) and group 2 (n = 79)] were identified (Figure 1A). Then, we calculated the differentially expressed genes (DEGs, group 2 vs. group 1), performed functional enrichment analysis, and identified inflammatory response, immune response, innate immune response, and neutrophil chemotaxis as the primarily enriched BP, and receptor activity, cytokine activity, and chemokine activity as the main MF (Supplementary Figure S1) implying a difference in immune response between the two groups. Consistently, GSEA analysis revealed that KEGG signaling pathways involved in immune response such as viral protein interaction with cytokine receptor, IL-17 signaling pathway, cytokine–cytokine receptor interaction, chemokine signaling pathway, and tumor necrosis factor (TNF) signaling pathway were significantly enriched in group 2 (Figure 1B), indicating a potentially inflamed microenvironment. Besides, the fraction of immune infiltration was calculated; group 2 was characterized by increased infiltration of B cell memory, T cell CD4 memory activated, Treg, monocyte, and neutrophil, while group 1 had increased infiltration of T cell CD4 memory resting and natural killer (NK) cell (Figure 1C). In particular, quantitative analysis of the activity of the anti-tumor immune response using TIP found that group 2 had increased tumor antigen releasing (step1) and immune cell recruitment activity (step 4) (Figure 1D) but not the infiltration of the immune cell into the tumor (step 5), tumor cell recognition (step 6), and elimination (step 7). In this way, the IRG-based classification yielded two subgroups of GBM, one of which (group 2) contained more inflammatory features in the microenvironment including the activation of immune-related signaling pathways and increased infiltration of neutrophils.

To identify signature genes characterizing such immune-inflamed microenvironment with prognostic significance, the gene expression profile of GBM was divided into 11 modules with functional relevance by WGCNA. Module “magenta” was significantly differentiated by subgroups and highly correlated with group 2 (cor = 0.65) (Supplementary Figures S2A, B). Functional enrichment analysis of the 608 genes included in “magenta” showed that they were mainly involved in leukocyte activation and immune activation (BP) as well as signaling pathways associated with immune cell activation and differentiation (biocarta) (Supplementary Figure S2D). Next, 409 DEGs upregulated in group 2 and the 161 genes that had a weighted correlation coefficient over 0.8 with group 2 were intersected, and 126 genes were obtained as candidate genes for the construction of the risk model (Supplementary Figure S2C). The univariate and multivariate Cox regression analysis identified S100A11, PLAUR, NCF2, FCGR2B, and ARPC1B to be prognostically significant (Table 1).

We first examined the relationship between the risk score and the well-studied clinicopathological and molecular parameters of GBM. We found that the mesenchymal subtype had a significantly increased risk score (p < 0.001) (Figure 2A), which has the poorest prognosis among the four transcriptional subtypes (Ceccarelli et al., 2016). Besides, the risk score was significantly increased in both IDH wild-type and ATRX wild-type tumors (p < 0.01), corroborating the association between the risk score and poor prognosis. Moreover, the risk score was moderately and significantly correlated with age and necrosis percent of the tissue slide (Figures 2B, C), two other factors associated with poor prognosis in GBM patients (Tan et al., 2020).

Thereafter, samples were split into the high- and low-risk groups based on the median of the risk score (Table 2). The distribution and status of OS showed that GBM patients with increased risk scores had unfavorable prognoses (Supplementary Figure S3A). Likewise, the K-M plots showed that high-risk scores predicted a significantly reduced OS in both cohorts (p < 0.05 in the TCGA cohort, p < 0.01 in the CGGA325 cohort) (Figure 3A), as well as in two other independent cohorts (Supplementary Figure S3B). The time-dependent predictive power of the risk model was also assessed. The 1-, 2-, 3-, and 4-year AUC values of the risk stratification were 71.24, 63.16, 73.59, and 73.07 in the TCGA cohort, and 59.15, 68.45, 66.45, and 73.20 in the CGGA cohort, performing better than other biomarkers including molecular subtype, IDH mutation status, gender, and age (Figure 3B). Further, the independent prognostic value of the risk score for different populations was determined. Univariate Cox regression analysis found that a high risk score was indicative of an unfavorable prognosis for the total population (HR = 1.55 and 1.74 in the TCGA and CGGA325 cohort) (Figure 3C). In the TCGA cohort, a high risk score was an independent risk factor for female patients and those aged less than 60 years. In the CGGA325 cohort, a high risk score was an independent prognostic factor for primary and recurrent GBM, male, patients aged less than 60, MGMT promoter methylated patients. Nomogram clearly showed that the high-risk group suggests reduced total score and survival (Figure 3D). In sum, these results suggested that the IRG-based risk model had a reliable prognostic value, with high-risk scores reliably indicating a poor prognosis of GBM patients.

To investigate the immunological differences between the high- and low-risk groups, the activation of immune-related signaling pathways was assessed. One hundred three gene sets involved in the inflammation and innate immune response, antigen presentation, CD 8 T cell activation, and cytotoxicity were retrieved from the MSigDB database. The ssGSEA score of these signaling pathways increased with the risk score, indicating the immune-activated status of the high-risk group (Figure 4A). Given that immune cells are key players in the immune response (Newman et al., 2015), the immune infiltration fraction was assessed. Consequently, neutrophil was positively correlated with the risk score in both cohorts (cor > 0.3) and significantly enriched in the high-risk group (p < 0.001) (Figures 4B, C), suggesting an inflamed microenvironment of the high-risk group. Besides, GSEA analysis found the enrichment of pathways including inflammatory responses, antigen presentation, CD8 T cell activation, and NK cell-mediated cytotoxicity in the high-risk group (Figure 4D), indicating the activation of the adaptive immune response. Further, there are well-established mRNA metrics reflecting the T-cell inflamed immune microenvironment and the cytolytic activity of CD8 T cell and NK cell (Rooney et al., 2015; Danaher et al., 2018). Correlation analysis found a moderate and significant correlation between the risk score and the two metrics (Supplementary Figure S4), corroborating the emergence of the adaptive immune response. Therefore, a high risk score was suggestive of an inflamed tumor microenvironment and the activation of the adaptive immune response.

Given the profound impact of tumor genomic variation on the tumor-immune interaction, we explored the genomic alterations associated with the IRG-based risk stratification. PTEN, TP53, TTN, EGFR, and MUC16 were the top five mutated genes, with a predominance of a missense mutation (Supplementary Figure S5A). We noted some mutations that have been reported to affect immune–tumor interactions; for instance, the mutation frequency of EGFR ranged from 15% in the high-risk group to 39% in the low-risk group, while the frequency of mutations in TP53 and PTEN ranged from 36% to 33% in the high-risk group to 29% and 25% in the low-risk group (Figure 5A). Fisher's exact test found a significantly increased mutation of GRM3, NF1, and RB1 in the high-risk group, and EGFR, TRPM2, and AHNAK2 in the low-risk group (Figure 5C). Among these, EGFR mutation was associated with increased plasma cell infiltration (two-tail Wilcoxon test, p = 0.05) (Supplementary Figure S5B). Besides, investigating co-occurrence and mutually exclusive mutations found EGFR and COL6A3 as a co-occurrence gene pair in the low-risk group (Figure 5B), the latter is involved in the production of type IV collagen and is associated with the strength of the vessel wall (Wu et al., 1996). Given the association of EGFR wild type with disrupted vasculature and inflamed microenvironment (Segura-Collar et al., 2021), and the relatively hypoxic microenvironment and hyper-angiogenesis activity of the high-risk group (Supplementary Figure S5D), the inflamed phenotype of the microenvironment in the high-risk group may be associated with EGFR wild-type-mediated aberrant angiogenesis and hypoxia, which partially explained the correlation between the risk score and necrosis.

In terms of CNV, the high- and low-risk groups shared recurring CNVs including 7p11.2 amplification and 9p21.3 deletion, whereas the low-risk group contained an additional significant focal amplification at 12q14.1, leading to the increased expression of CDK4, TSPAN31, MARCH9, and AGAP2 that located at this region (Figure 5D and Supplementary Figure S5C). Despite the individual gene mutations and chromosome alterations, the statistical indicators that related to immune response were also explored. The high-risk group had significantly increased somatic copy number amplification and deletion burden, while the TMB was comparable between the two groups (Figures 5E, F), consistent with our speculation that necrosis may be the source of the activated immune response in the high-risk group.

The immune phenotype of tumors is related to their responsiveness to ICB therapy (Chen and Mellman, 2017); we, therefore, assessed the ICB responsiveness of the high- and low-risk groups. The high-risk group had significantly increased expression of PD1 and CTLA4 as well as their ligands (PD-L1/2, CD80/CD86) (Figures 6A, C). Besides, the primary immunodeficiency pathway (KEGG), CTLA4 pathway (Biocarta), and cancer immunotherapy by PD1 blockade (WP) were highly enriched in the high-risk group (Supplementary Figure S6), indicating that immune checkpoints play a role in GBM immune evasion. Tumor response to anti-PD1 and anti-CTLA4 therapies was predicted using TIDE and further validated using SubMap. Consequently, a high risk score failed to predict response to anti-PD1 (Bonferroni corrected p = 0.160 and 0.128 in the TCGA and CGGA cohort) or anti-CTLA4 treatment (Bonferroni corrected p = 0.208 and 0.120, respectively) (Figures 6B, D).

In addition, we assessed the relationship between the risk stratification and the effectiveness of radiotherapy and chemotherapy. Samples were firstly split into two groups according to their radiotherapy status. For patients receiving radiotherapy, an increased risk score suggested a significantly poor prognosis (p = 0.013 and 0.0022 in the TCGA and CGGA cohort, respectively), while the prognostic significance of the risk stratification for patients without radiotherapy was not significant (Figure 7A), indicating that the IRG-based risk model had the potential to predict the efficacy of radiotherapy. Similarly, we grouped samples according to the chemotherapy they received. In the TCGA cohort, lower risk scores for patients treated with temozolomide (TMZ) or adjuvant TMZ predicted increased OS (p = 0.038 and 0.048 for patients receiving TMZ and adjuvant TMZ in the TCGA cohort) (Figure 7B). However, the conclusion was not validated in the validation dataset (Figure 7B), suggesting the need to include more samples and to further differentiate chemotherapy modalities to confirm the prediction of chemotherapy efficacy by the risk stratification.