Identification of GBM endothelial cells (ECs) Immune-Related Genes
The detailed workflow of the risk model construction and downstream analysis is shown in Figure 1. As shown in the Venn diagram, a total of 11 prognostic ECs immune-related genes (EIRGs), which were indicated to have significant prognostic value (P<0.05), were screened out after intersecting the identified 59 GBM ECs related prognostic genes and the identified 438 immune-related prognostic genes (Figure 2a and 1c). Subsequently, we performed to investigate the correlation among the identified EIRGs with a threshold of correlation coefficient was 0.2. The correlation analysis results in this study showed that both SBDS and RBP1 were negatively correlated with PLXND1, while other EIRGs exhibited different positive correlation(Figure 2b).
Construction and validation of the prognostic-related 6-EIRG signature
To build the EIRGs signature for forecasting the OS of GBM patients, we performed the least absolute shrinkage and selection operator (LASSO) Cox regression analysis on the basis of the 11 prognostic EIRGs in the TCGA cohort. Next, we identified 6 prognostic-related EIRGs (PLXND1, CALCRL, RBP1, SBDS, CFH, and SOCS3) to build the risk model, and the coefficients of these genes were used to calculate the risk score. (Figure 3a). The risk score= PLXND1×0.0521 + CALCRL× (-0.0292) + RBP1×0.0878 + SBDS×0.0912 + CFH×0.0656+ SOCS3×0.0831. According to the median value of the risk scores, GBM patients were divided into the high-risk group and low-risk group, and then performed Kaplan-Meier survival analysis. Risk score, survival status distributions, and the relative expression standards of the 6 prognostic-related EIRGs for each patient are shown in Figure 3b. The Kaplan-Meier survival analysis showed that patients in the high-risk group displayed worse overall survival (OS) compared to those in the low-risk group (P < 0.001), indicating that the risk score could predict the prognosis in the TCGA cohort (Figure 3c). To validate the prognostic ability of the prognostic risk model, we calculated risk scores for GBM patients in the CGGA cohort using the same formula (Figure 3d). The results of the Kaplan-Meier survival curve in the CGGA cohort were consistent with our findings in the TCGA cohort: GBM patients in the high-risk group had significantly worse OS than those in the low-risk group (P < 0.001), which made our results convincing (Figure 3e).
Evaluation the correlation between risk score and the clinicopathological characteristics
The heatmap shows clinical characteristics and the expression of the 6 prognostic-related EIRGs in high- and low-risk patients in the TCGA cohort (Figure 4a). We observed significant differences between the high- and low-risk groups associated with IDH mutation status, stromal score, immune score, estimate score as well as tumor purity. To determine whether the prognostic value of the 6 prognostic-related EIRGs signature-based risk score is independent of other clinical features, univariate and multivariate Cox regression analyses were performed to analyze with risk score and other clinical features, such as including age, gender, IDH mutation status, and radiation therapy. Univariate and multivariate Cox analysis indicated risk score remained to be an independent prognostic factor (P < 0.05, Figure 4b and 4c).
Nomogram establishment and evaluation
As shown in Figure 5a-5d, a nomogram for predicting 0.5-, 1- and 2-year survival rates of GBM patients was established according to age, gender, IDH mutation status, radiation therapy, and risk score. The C-index of the prognostic nomogram was 0.707, and the calibration curve of the 0.5-, 1- and 2-year OS indicated that our nomogram model had a good predictive ability. Importantly, we found similar results in the CGGA cohort (Figure 5e-5g).
Estimation of the tumor immune microenvironment using the prognostic-related EIRGs risk model
As shown in Figure 6a, the immune, stroma, and ESTIMATE scores were significantly higher (p < 0.05) whereas tumor purity was lower in the high-risk group compared with the low-risk group (p < 0.05). To elucidate the difference in the immune microenvironment of GBM, we analyzed immune cell infiltrating density between high- and low-risk groups. The results indicated that high immune infiltration was strongly related to the high-risk group (Figure 6b). Consistent with our previous results, the risk score of the prognostic model had a significant positive correlation with immune cells infiltration, including regulatory T cell (Tregs), Th1 (T helper cell, type 1), activated dendritic cell, macrophages, Mast cell, NK natural killer T cells (NKT), MDSCs, and neutrophils (Figure 6c). Notably, the results from the correlation analyses showed that the risk score was positively correlated with activated CD4+ T cells, activated CD8+ T cells, activated B cells, activated dendritic cells, macrophages, regulatory T cells (Treg), and natural killer cells (P < 0.05) (Supplementary Figure 1). Among the 29 immune gene sets, immune-related functional cells such as activated dendritic cells (aDCs), B cells, DCs, interdigitating dendritic cells (iDCs), macrophages, neutrophils, T helper cell, follicular helper T cells (Tfh), Th1 cells, Th2 cells, and regulatory T cells (Treg) showed higher ssGSEA scores in the high-risk patients (Figure 6d). Similarly, immune pathways such as antigen-presenting cell (APC) co-inhibition, APC co-stimulation, CC chemokine receptor (CCR), Check-point, cytolytic activity, human leukocyte antigen (HLA), inflammation-promoting, para-inflammation, T-cell co-inhibition, T-cell co- stimulation, Type II IFN Response, and tumor-infiltrating lymphocytes (TIL) exhibited higher ssGSEA scores in the high-risk group (Figure 6d). Additionally, similar results were obtained in the CGGA cohort of GBM patients (Supplementary Figures 2 and 3). These data indicated that the 6 prognostic-related EIRGs signature were highly associated with immune infiltration.
GSEA Analysis of 6-EIRGs signature
Next, GSEA analysis was performed to assess the involvement of biological processes, immune-related pathways, and cancer-related hallmarks of each GBM sample. The top 10 KEGG pathways, GO biological processes (GOBP), and hallmark gene sets were identified based on the P-value < 0.05 in the TCGA and CGGA cohorts. The GSEA analysis showed that most of the KEGG pathways, GOBP, and hallmark gene sets enriched in the high-risk group were associated with the immune and inflammatory responses as well as the cancer-related pathway, including Toll-like receptor signaling pathway, JAK/STAT signaling pathway, NOD-like receptor signaling pathway, IL-2-STAT5 signaling pathway, IL-6-JAK-STAT3 signaling pathway, TNFA signaling via NFKB, inflammatory response, activation of the immune response, and others in the TCGA and CGGA cohorts (Figure 7a and 7b). These findings indicate that prognostic-related EIRGs may take part in regulating the tumor immune environment and immunologic biological processes of GBM.
Sensitivity of Different Therapies in the high- and low-risk groups.
Next, we further predicted the likelihood of response to immunotherapy in the high-risk and low-risk groups using the TIDE algorithm. TIDE scores were significantly higher in the high-risk group compared with the low-risk group, which indicated the low-risk group was much more sensitive to immunotherapy (Figure 8a). Subsequently, we found that GBM patients with high-risk group exhibited higher Microsatellite instability (MSI) and Dysfunction, and lower expression of Exclusion compared to the patients with the low-risk group (Figure 8b-d). Additionally, a significant correlation relationship between the expression level of the key EIRGs and chemotherapy drug sensitivity are displayed in Figure 8 and ranked by the p-value, selected by p<0.05. Notably, SOCS3 was positively correlated with the sensitivity of Bleomycin, Sonidegib, Irofulven, and Olaparib and the resistance of Selumetinib, Encorafenib, Cobimetinib, and ARRY−162. Highly expressed PLXND1 were more resistant to Palbociclib (Cor=−0.366, p=0.004) and METHOTREXATE (Cor=−0.355, p=0.005) (Figure 9). Afterwards, we further explore the differences in the treatment response to temozolomide (TMZ) chemotherapy in the high- and low-risk groups. As shown in Figure 10a, the risk score was significantly positively correlated with the TMZ (P<0.05). The estimated IC50 corresponding to the low-risk group was significantly lower than that of the high-risk group, indicating that GBM patients from the low-risk group were more sensitive to TMZ (P < 0.05, Figure 10b). Notably, we observed similar results in the CGGA cohort (Figure 10c and 10d). These results indicated that the prognostic-related EIRGs signature could predict the treatment response to immunotherapy and TMZ chemotherapy, and guide effective therapy.