Endothelial Cells Immune-Related Genes Risk Signature Correlated With Tumor Immune Microenvironment Predicts Therapeutic Response and Prognosis in Glioblastoma

Glioblastoma (GBM) is the most common and aggressive primary tumor of the central nervous system with high recurrence and extremely poor prognosis. Multiple recent studies have indicated a pivotal correlation between GBM prognosis and immune-related risk signature. Nevertheless, the potential value of endothelial cells (ECs) immune-related genes (EIRGs) in prognosis, immune inltration, and their correlation with therapeutic response to immunotherapy and TMZ chemotherapy remain obscure, especially in GBM. Here, we screened out 11 EIRGs after intersecting the identied 59 GBM ECs related prognostic genes and the identied 438 immune-related prognostic genes. A prognostic-related 6-EIRGs signature was established through univariate Cox analysis and least absolute shrinkage and selection operator (LASSO) Cox regression analysis, and patients in the high-risk group were signicantly worse overall survival (OS) compared to those in the low-risk group. Additionally, univariate and multivariate Cox regression analysis conrmed that risk score was an independent predictor of OS in patients with GBM. The nomogram which comprised age, gender, IDH mutation status, radiation therapy, and risk score yielded a strong predictive ability of 0.5-, 1-, and 2-year OS for GBM patients. Our results demonstrate that the EIRGs signature, which is associated with immune cell inltration, may play a regulatory role in the immunobiological process of TIME. Prognostic-related 6-EIRGs signature is a promising classication index for predicting the drug sensitivity to immunotherapy and TMZ chemotherapy, suggesting that EIRGs signature may serve as a biomarker to stratify patients who will benet from immunotherapy and chemotherapy. role of prognostic-related EIRGs in the prognosis, tumor immune microenvironment, and drug sensitivity to immunotherapy and TMZ chemotherapy of GBM, which provides a new research strategy for the early diagnosis, treatment, and prognosis of GBM. We constructed a prognostic-related 6-EIRGs signature with good prognostic value and which can serve as an independent prognostic indicator for GBM. Notably, this is the rst study to explore the relationship between prognostic-related EIRGs and TIME, and drug sensitivity to immunotherapy and TMZ chemotherapy in the GBM. Our results indicate that the risk stratication based on risk signature might help distinguish GBM patients that could benet from immunotherapy and chemotherapy drug TMZ. related gene pairs; KEGG: Kyoto Encyclopedia of Genes and Genomes; LASSO: Least Absolute Shrinkage and Selection Operator; MSI: Microsatellite instability; OS: Overall survival; ssGSEA: Single Sample Gene Set Enrichment Analysis; scRNA-seq: Single-cell RNA sequencing; TCGA: The Cancer Genome Atlas; TIME: Tumor immune microenvironment; TMZ: Temozolomide; TIDE: Tumor Immune Dysfunction and Exclusion.


Introduction
Glioblastoma (GBM) is the most common and aggressive primary tumor of the central nervous system, which is characterized by a high degree of intra-tumoral heterogeneity and resistance to treatment [1][2][3].
Despite recent advances in multimodal standard therapy, including surgical resection, radiotherapy, and temozolomide (TMZ) chemotherapy, and immunotherapy, the prognosis of patients with GBM is still poor, with a median survival time of around 15 months [4][5][6]. So far, the alkylating drug temozolomide (TMZ) is the mainstay of rst-line chemotherapy for the treatment of GBM, but drug resistance to TMZ chemotherapy and drug delivery across the blood-brain barrier (BBB) are the major obstacles to the effective treatment of GBM [7,8]. Additionally, growing evidence has indicated that immunotherapy has revolutionized the treatment of many solid-tumor malignancies and generated great hope for GBM, but without success until now [9][10][11]. Therefore, novel biomarkers that can predict the drug sensitivity to immunotherapy and TMZ chemotherapy in GBM should be sought. Meanwhile, the tumor immune microenvironment (TIME) also plays a critical role in tumorigenesis, progression, therapeutic e cacy, and clinical outcomes of GBM patients [12][13][14]. Thus, elucidating the role of the components and signatures of the tumor immune microenvironment and the interaction between cancer cells and the TIME may help to understand the pathogenesis of tumors and guide the development of novel therapeutic interventions.
Cerebral endothelial cells (ECs) are the major components of the BBB. ECs are connected by complex tight junctions and exhibit very low levels of vesicular transcytosis, leading to a very low vascular permeability [15,16]. Additionally, ECs exert the barrier characteristics of the BBB in binding and e ux of various small lipophilic foreign compounds and drugs through the speci c transport network of e ux transporters, thus limiting the development and absorption of neurotherapeutic drugs [17]. Tumor heterogeneity complicates the treatment of patients with cancer [18], therefore, we should further understand how ECs affect the prognosis of cancer and how it fosters the emergence of drug resistance, so as to provide novel breakthroughs and therapeutic strategies for targeted cancer therapy in the future.
Of note, Xie et. al have explored the heterogeneity of the human BBB and its molecular alteration in GBM by scRNA-seq-based molecular atlas of human brain ECs [19]. A previous study has revealed that the 1p/19 codeletion-associated immune prognostic signature could accurately predict the prognosis of glioma [20]. Moreover, our earlier study had demonstrated that immune-related gene pairs (IRGPs) signature, which is associated with immune cell in ltration, can be used for individual OS predictions for glioma patients [21]. Nevertheless, little is known about how ECs-immune-related genes (EIRGs) affect GBM prognosis and how the relationship between EIRGs and in ltrating immune cells in the TIME, and their correlation with therapeutic response to immunotherapy and TMZ chemotherapy.
This study was conducted with an aim to identify EIRGs, whose expression is correlated closely with the prognosis of GBM patients in The Cancer Genome Atlas (TCGA) database. To further assess the prognostic value of these prognostic-related EIRGs, univariate Cox regression analysis and LASSO Cox regression analysis were performed to construct an ECs-immune-related gene prognostic model to facilitate the prediction of GBM prognosis. Subsequently, a validation model from the Chinese Glioma Genome Atlas (CGGA) database was constructed to evaluate the feasibility of the risk signature based on the ECs-immune-related genes. Moreover, we analyzed the abundances of the tumor-in ltrating immune cells in GBM samples between the risk groups to explore the association between prognostic-related EIRGs and the tumor immune microenvironment (TIME). Finally, we explored these EIRGs signature association with therapeutic response to immunotherapy and TMZ chemotherapy.

Construction of EC-immune-related genes (EIRGs) Risk Model
Univariate Cox regression analyses with P-value < 0.05 were performed to identify 59 GBM ECs related prognostic genes and 438 immune-related prognostic genes. To screen the EIRGs for GBM patients, the above-identi ed ECs related prognostic genes and immune-related prognostic genes were intersected in the TCGA and CGGA datasets. From this, we identi ed a total of 11 prognostic EIRGs, with a statistically signi cant correlation (P<0.05), which we developed as a potential risk signature for the selection of candidate risk EIRGs with the LASSO Cox regression analysis. The TCGA set was utilized to construct an EIRGs-related prognostic model, and the CGGA set was used to validate the prognostic accuracy of this established model. Finally, 6 EIRGs with their coe cients (β) were determined according to the minimum criteria. Risk score for each patient in the TCGA and CGGA cohorts were calculated according to the following formula: Risk score = (expression of gene A1*β1) + (expression of gene A2*β2) + (expression of gene A3*β3) +... (expression of gene An*βn). Then, we divided patients into two groups: the high-risk group and the low-risk group based on the median value of the risk score. The Kaplan-Meier survival analysis was used to evaluate the differences in overall survival (OS) between high and low-risk score groups by using the "survival" R package.
Independence of the prognostic-related EIRGs risk model and Nomogram analysis Additionally, univariate and multivariate Cox regression analyses were used to determine whether the EIRGs-related risk scores were independent of other clinicopathological and molecular features (age, gender, IDH mutation status, radiation therapy, and risk score) as independent predictors for patients with GBM. Nomogram analysis for predicting the probability of 0.5-, 1-, and 2-year OS for GBM patients by integrating risk score and clinical features (age, gender, IDH mutation status, and radiation therapy) was constructed and veri ed in the TCGA cohort as well as CGGA cohort using the "rms" R package. Moreover, we plotted the calibration curve to validate the predictive accuracy and reliability of the nomogram.
Afterward, heatmaps were drawn to investigate the correlation between EIRGs-related risk scores and clinicopathological and molecular features using the "ComplexHeatmap" R package [22].

Evaluation of Tumor Microenvironment In ltration Patterns and Immune Cells In ltration
The stromal score, immune score, estimate score, and tumor purity were calculated with the ESTIMATE algorithm using the "estimate" R package [23]. Subsequently, to investigate the immune in ltration landscape of GBM, we utilized the ssGSEA (Single Sample Gene Set Enrichment Analysis) to assess 23 types of tumor-in ltrating immune cells and further revealed a correlation between tumor-in ltrating immune cells and different risk groups. Subsequently, the 'GSVA' R package (v.1.38.2) was used for ssGSEA to evaluate immune-related functions or pathways between high-risk and low-risk groups [24]. And then, Gene Set Enrichment Analysis (GSEA) of prognosis-related EIRGs signature was performed for the assessment of the enrichment of biological processes and pathways among high-risk and low-risk groups [25]. The KEGG gene sets (c2.cp.kegg.v7.4.symbols.gmt), Go-Biological Processes (GO.BP) gene sets (c5.go.bp.v7.4.symbols.gmt), and HALLMARK gene sets (h.all.v7.4.symbols.gmt) were downloaded from the Molecular Signatures Database (MSigDB, http://www.gsea-msigdb.org/gsea/downloads.jsp) to run the GSEA analysis. The "clusterPro ler" R package (V3.18.1) was used to perform the top 10 functional annotations and pathway analysis on KEGG, GO, and hallmark pathway on EIRGs with an adjusted p < 0.05, and FDR < 0.05 [26].

Therapeutic Response Prediction
The Tumor Immune Dysfunction and Exclusion (TIDE) score of patients with GBM from the TCGA dataset was downloaded from the TIDE website (http://tide.dfci.harvard.edu/) for predicting the response to immunotherapy. Higher TIDE predictive scores were associated not only with greater potential of tumor immune evasion but also with worse immunotherapy response. The correlation between the identi ed EIRGs signature and Drug response was predicted by CellMiner (https://discover.nci.nih.gov/cellminer/home.do). To examine the sensitivity of each sample to chemotherapeutic response to TMZ, the "pRRophetic" R package was utilized to predict the half-maximal inhibitory concentration (IC50) of TMZ for the comparison of the different risk groups [27]. P < 0.05 was considered statistically signi cant.

Identi cation of GBM endothelial cells (ECs) Immune-Related Genes
The detailed work ow 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 signi cant prognostic value (P<0.05), were screened out after intersecting the identi ed 59 GBM ECs related prognostic genes and the identi ed 438 immune-related prognostic genes (Figure 2a and 1c). Subsequently, we performed to investigate the correlation among the identi ed EIRGs with a threshold of correlation coe cient 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 identi ed 6 prognostic-related EIRGs (PLXND1, CALCRL, RBP1, SBDS, CFH, and SOCS3) to build the risk model, and the coe cients 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 ndings in the TCGA cohort: GBM patients in the high-risk group had signi cantly 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 highand low-risk patients in the TCGA cohort (Figure 4a). We observed signi cant 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 Cindex 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 signi cantly 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 in ltrating density between high-and low-risk groups. The results indicated that high immune in ltration was strongly related to the high-risk group (Figure 6b). Consistent with our previous results, the risk score of the prognostic model had a signi cant positive correlation with immune cells in ltration, 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), in ammation-promoting, parain ammation, T-cell co-inhibition, T-cell co-stimulation, Type II IFN Response, and tumor-in ltrating 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 in ltration.
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 identi ed 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 in ammatory 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, in ammatory response, activation of the immune response, and others in the TCGA and CGGA cohorts (Figure 7a and 7b). These ndings 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 signi cantly 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 signi cant 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 signi cantly positively correlated with the TMZ (P<0.05). The estimated IC50 corresponding to the low-risk group was signi cantly 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.

Discussion
GBM is the most aggressive and lethal primary brain tumor with heterogeneous clinical features and a generally poor clinical outcome [28]. The treatment of GBM remains a great clinical challenge using temozolomide (TMZ) based on conventional chemotherapy and is mainly due to drug resistance [29]. Previous studies have reported that Xie et. al have performed scRNA-seq on GBM ECs to characterize the heterogeneity of gene expression signatures of different ECs, which will provide key information for designing rational therapeutic regimens and optimizing drug delivery [19]. However, little is known about the prognostic-related EIRGs signature with tumor immune microenvironment and drug sensitivity to immunotherapy and TMZ chemotherapy. In this study, we successfully constructed a prognostic risk signature with six identi ed prognostic-related EIRGs (PLXND1, CALCRL, RBP1, SBDS, CFH, and SOCS3), which divided GBM patients into high-and low-risk groups. Multivariate Cox regression analysis showed that prognostic-related EIRGs signature-based risk score was an independent prognostic factor for OS (HR: 5.207; 95% CI: 1.338-20.260; P < 0.05) after adjusting for other clinical variables. Moreover, combining risk score with age, gender, IDH mutation status, and radiation therapy, we constructed a nomogram, and it had a robust ability to predict 0.5-year, 1-year, and 2-year OS in GBM patients.
Tumor formation involves the co-evolution of neoplastic cells together with the TME surrounded by extracellular matrix, tumor vasculature, and immune cells, which may constitute the basis of tumor microenvironmental heterogeneity due to the diversity of these cells [30]. Numerous studies had con rmed that immune in ltration is related to the prognosis of GBM, and tumor immune microenvironment signi cantly in uences therapeutic response and clinical outcome [31,32]. Intratumor heterogeneity remains a major obstacle for more effective anti-cancer therapy and personalized medicine [33]. Multiple recent studies have shown that the heterogeneity of stromal cells and immune cells in the TIME also plays a critical role in cancer cell proliferation, angiogenesis, immune evasion, metastasis, and responses to treatment [34].In our study, we have demonstrated that patients in the high-risk group had higher immune, stroma, and ESTIMATE scores, lower tumor purity, and a higher abundance of immune cell in ltration. Additionally, the GSEA analysis for the 6 EIRG-signature suggested that patients in the highrisk group were signi cantly enriched in biological processes, immune-related pathways, and cancerrelated hallmarks, including immune responses (Toll-like receptor signaling pathway, TNFA signaling via NFKB, IL2-STAT5 signaling, IL6-JAK-STAT3 signaling, in ammatory response, and interferon-gamma response), cellular stress (hypoxia and apoptosis), and signal transduction (JAK/STAT signaling pathway and cytokine-cytokine receptor interaction). Notably, JAK/STAT signaling pathway plays an important role in immune regulation, apoptosis, drug resistance, tumor proliferation, migration, and invasion in glioma [35,36]. Moreover, previous studies also demonstrated that hypoxia may represent an early sign for GBM recurrence, which might become useful in the development of new combined diagnostictherapeutic approaches of recurrent GBM [37]. These results indicated that the prognostic-related EIRGs risk signature was closely associated with immunologic biological processes and the regulation of tumor immune microenvironment in GBM.
Therapeutic resistance re ects active tumor evolution, and environment-mediated resistance re ects the dynamic interaction between tumor cells and their surrounding microenvironment [38]. A growing body of evidence has demonstrated that temozolomide (TMZ) is the most commonly used chemotherapeutic drug used to treat GBM, and its function is mainly achieved by inducing DNA damage in tumor cells [39,40]. Nevertheless, most GBM patients are resistant to the drug TMZ [41]. In our study, we found that the estimated IC50 corresponding to the high-risk group was signi cantly higher than that of the low-risk group, suggesting that GBM patients from the high-risk group based on prognostic-related EIRGs could be more resistant to TMZ chemotherapy drug, which may partially explain worse survival outcomes of patients in the high-risk group. Interestingly, patients in the low-risk group who exhibited lower TIDE scores showed more likely effectiveness for immunotherapy than the high-risk group, which indicated that patients from the high-risk group might have a higher chance of anti-tumor-immune escape. Our analysis demonstrates that prognostic-related 6-EIRGs signature is a strong discriminative tool for examining the sensitivity of GBM patients to immunotherapy and TMZ chemotherapy. Recent clinical and experimental studies have demonstrated that the neoadjuvant administration of PD-1 blockade might generate enhanced antitumor immune responses and induce functional activation of tumor-in ltrating lymphocytes, which may represent a promising strategy in the treatment of GBM [42,43]. Additionally, Park et al. have reported that PD-1 blockade induces a long-term therapeutic response, and combination with TMZ further enhances antitumor e cacy [44]. Thus, in our study, by exploring the potential relationship between the risk strati cation and drug sensitivity to immunotherapy and TMZ chemotherapy, we can further provide new insights and achieve optimal clinical bene t for the combination of targeted immunotherapy and TMZ chemotherapy in the GBM.
Remarkably, we are also aware of some limitations in this study. First, our prognostic model was constructed based on open accessed databases without our cohort. Second, the biological and regulatory mechanism of prognostic-related EIRGs in the GBM immune microenvironment and their interaction with tumor-in ltrating immune cells needed to be further veri ed using more in vivo and in vitro experiments. Third, there were no recruited cohorts to verify the prognostic value of our prognostic-related EIRGs risk signature and drug sensitivity to immunotherapy and TMZ chemotherapy in high-and low-risk groups.

Conclusion
Collectively, our study provides a comprehensive insight for revealing the important role of prognosticrelated EIRGs in the prognosis, tumor immune microenvironment, and drug sensitivity to immunotherapy and TMZ chemotherapy of GBM, which provides a new research strategy for the early diagnosis, treatment, and prognosis of GBM. We constructed a prognostic-related 6-EIRGs signature with good prognostic value and which can serve as an independent prognostic indicator for GBM. Notably, this is the rst study to explore the relationship between prognostic-related EIRGs and TIME, and drug sensitivity to immunotherapy and TMZ chemotherapy in the GBM. Our results indicate that the risk strati cation based on risk signature might help distinguish GBM patients that could bene t from immunotherapy and chemotherapy drug TMZ.

Declarations Ethics approval
This article is not involved in any studies with human participants or animals performed by any of the authors.

Consent for publication
All authors consent to the publication of this study.

Disclosure of potential con icts of interest
The authors declare that they have no competing nancial interests or personal relationships that could have appeared to in uence the work reported in this paper.          The correlation between chemotherapy drug sensitivity and the key EIRGs (SOCS3, PLXND1, and CALCRL) in TCGA cohort. The scatter plots are ranked by the p-value, selected by p<0.05.

Figure 10
Association between risk score and chemosensitivity of TMZ in GBM. (a) The risk score was signi cantly positively correlated with the TMZ in the TCGA cohort. (b) The estimated IC50 corresponding to the highrisk group was signi cantly higher than that of the low-risk group in the TCGA cohort. (c) The risk score was signi cantly positively correlated with the TMZ in the CGGA set. (d) The estimated IC50 corresponding to the high-risk group was signi cantly higher than that of the low-risk group in the CGGA set. IC50: half-maximal inhibitory concentration; TMZ: Temozolomide.