Prognostic value and immune relevancy of a combined autophagy-, apoptosis- and necrosis-related gene signature in glioblastoma

Glioblastoma (GBM) is considered the most malignant and devastating intracranial tumor without effective treatment. Autophagy, apoptosis, and necrosis, three classically known cell death pathways, can provide novel clinical and immunological insights, which may assist in designing personalized therapeutics. In this study, we developed and validated an effective signature based on autophagy-, apoptosis- and necrosis-related genes for prognostic implications in GBM patients. Variations in the expression of genes involved in autophagy, apoptosis and necrosis were explored in 518 GBM patients from The Cancer Genome Atlas (TCGA) database. Univariate Cox analysis, least absolute shrinkage and selection operator (LASSO) analysis, and multivariate Cox analysis were performed to construct a combined prognostic signature. Kaplan–Meier survival, receiver-operating characteristic (ROC) curves and Cox regression analyses based on overall survival (OS) and progression-free survival (PFS) were conducted to estimate the independent prognostic performance of the gene signature. The Chinese Glioma Genome Atlas (CGGA) dataset was used for external validation. Finally, we investigated the differences in the immune microenvironment between different prognostic groups and predicted potential compounds targeting each group. A 16-gene cell death index (CDI) was established. Patients were clustered into either the high risk or the low risk groups according to the CDI score, and those in the low risk group presented significantly longer OS and PFS than the high CDI group. ROC curves demonstrated outstanding performance of the gene signature in both the training and validation groups. Furthermore, immune cell analysis identified higher infiltration of neutrophils, macrophages, Treg, T helper cells, and aDCs, and lower infiltration of B cells in the high CDI group. Interestingly, this group also showed lower expression levels of immune checkpoint molecules PDCD1 and CD200, and higher expression levels of PDCD1LG2, CD86, CD48 and IDO1. Our study proposes that the CDI signature can be utilized as a prognostic predictor and may guide patients’ selection for preferential use of immunotherapy in GBM.

Nervous System, the diffuse gliomas include WHO grade II and grade III astrocytic tumors, grade II and III oligodendrogliomas, grade IV glioblastomas, and related diffuse gliomas of childhood [1]. Various grades of gliomas differ considerably in tumor pathology, tumor development, and patient prognosis. Glioblastoma (GBM) is considered the most malignant and invasive primary intracranial tumor, with a high risk of recurrence [2][3][4]. Patients with GBM have a very poor prognosis, with an average overall survival of merely 12-15 months [5]. In spite of recent advances in standard treatment, including surgery, chemotherapy, radiotherapy, and the achievement in targeted therapies and immunotherapies over the past several years, GBM still carries a dismal prognosis with poor survival [6][7][8][9][10]. Therefore, novel prognostic approaches to pick out patients with high risks are warranted to further help therapeutic options for GBM patients.
Cell death is a critical process that maintains physiological homeostasis in multicellular organisms. Recently, numerous studies have revealed that the tumor microenvironment (TME) could be affected by dying and dead cancer cells for their potent immunomodulatory effects [11,12]. Dying/death cell leads to redundant bioactive factors release, which can either improve or weaken anticancer immunity. Cell death can also result from severe conditions existing in the TME and may significantly alter tumor progression. Researches established that multiple cell death pathways tended to play a part in the treatment response of tumors [13]. The three classically known cell death pathways are autophagy, apoptosis, and necrosis [14,15]. Autophagy, the process of self-degradation of cellular components, is upregulated when stimulated by extra-or intracellular stress and signals, such as starvation and growth factor deprivation [16]. Consequently, the chronic stress induction can cause irreversible damage, leading to apoptosis or necrosis [17]. Apoptosis is a programmed cell death process with distinct morphological characteristics and energy-dependent biochemical mechanisms [18]. It represents a critical pathway for eliminating cells that are not vital and protects against cells that have received significant genotoxic damage, and is instrumental in immune function [19]. Necrosis, the aftereffect of irreversible cellular damage, is recognized by the rapid destruction of plasma membrane followed by cytoplasmic leakage and the spilling of inflammatory cellular contents into the TME [20]. In short, the cell death processes dysregulation can significantly affect tumorigenesis.
GBM is a highly heterogeneous tumors with multiple subtypes, functionally different for their specific immunological landscapes, such as differences in T cell infiltration and macrophage composition, which require different treatment regimens [21,22]. Immune checkpoints, widely studied in recent years, are immunosuppressive molecules that avoid normal tissue damage and destruction primarily by modulating the immune response of T cells. Therefore, activating immune checkpoints may cause immune tolerance during tumor progression. Immune checkpoint inhibitors (ICI) can evade anti-tumor immune response, act on the tumor, and restrict its growth. The most effective ICI, anti-PD-1/ PD-L1, has been approved in non-small cell lung cancer, colon cancer, and melanoma [23]. However, recent clinical trials indicated that anti-PD-1/PD-L1 treatment might not benefit the clinical outcome of GBM without patients' selection [24]. Besides, contrary to other cancers, there is still no immunotherapy approved by Food and Drug Administration (FDA) for GBM. One of the arguments challenging GBM immunotherapy is its highly immunosuppressive TME. Thus, identifying regulators of the brain TME could help discover promising new targets for therapeutic intervention. Studies analyzed current clinical trial failures and demonstrated that biomarkers for appropriate patient selection for immunotherapy appeared hopeful in GBM treatment [25,26]. Recently, several novel prognostic markers for GBM patients have been identified through multiomic analysis and differential expression profiles. However, most of these studies are mathematical analyses based on whole-scale genetic or transcriptomic data, and there is still a lack of specific research on multiple biological pathways [27][28][29]. Therefore, comprehensive recognition of the characteristics of TME cell infiltration mediated by multiple cell death pathways is needed to deepen our understanding of TME immune regulation and help design enhanced treatment for GBM patients.
In this study, GBM patients were stratified based on a combination of autophagy-, apoptosis-and necrosisrelated gene signatures along with the characteristics of their immune response to facilitate the prediction of individualized survival and a superior treatment scheme.

Patient population and multiomic data acquisition
The genomic expression and clinical data of GBM patients in the TCGA database were retrieved from GlioVis (http:// gliov is. bioin fo.cnio.es/) [30]. The RNA sequencing data of the Illumina HiSeq 2000 platform and the clinical data were accessed from the Chinese Glioma Genome Atlas (CGGA) database (http:// www. cgga. org. cn) [31]. We included 518 GBM patients from TCGA and 137 patients from the CGGA database after excluding those without survival information. The data in the TCGA database were analyzed as the training cohort, and data from the CGGA dataset were used for validation. The Trans Per Million values of RNA-Seq and robust multichip analysis-processed values of microarray data were log2 transformed and then normalized by the scale method in R to make the data comparable between platforms [32]. Furthermore, 505 GBM patients' copy number alteration (CNA) data were obtained from the TCGA database. Using the RCircos package in R, Circos plots visualized the CNA summary results and determined chromosomal alterations [33]. Additionally, the somatic mutation data of 390 GBM patients were acquired on the basis of the whole-exome sequencing platform from the TCGA database. The data were analyzed and uncovered by utilizing the maftools package in R [34].

Gene Expression Analysis to Determine Cell Death Index (CDI)
To clarify the prognostic association of cell death-related genes in GBM, autophagy-, apoptosis-and necrosisrelated gene lists were accessed from the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases through the Gene Set Enrichment Analysis (GSEA) website(https:// www. gsea-msigdb. org) [35]. The KEGG dataset of apoptosis-related genes (n = 87) (Table S1) and GO gene list of necrosis (n = 49) (Table S2) were obtained. Autophagy-related genes were downloaded from the GO dataset and the Human Autophagy Database (HADb, http:// www. autop hagy. lu/index.html). The two gene sets were merged into one(n = 495) (Table S3). Univariate Cox regression models were used in each cell-death pathway to screen genes associated with OS in the TCGA datasets. The prognostic gene combination for establishing the index was screened out with LASSO regression. To further determine the optimal genes, a multivariate cox regression model was then performed using the "step" function in R. Subsequently, 140 patients, with the highest or lowest expression level of specific pathway genes, were selected from each cell death group. Significant gene signatures from individual cell death pathways were chosen to create a combined prognostic model to construct a cell death index (CDI). The latter was formed on the basis of a linear combination of the regression coefficient acquired from the multivariate Cox regression model and the genes expression levels. The CDI formula was calculated as follows: Risk score = (exprgene1 × Coef-gene1) + (exprgene2 × Coefgene2) + … + (exprgene n × Coefgene n). GBM patients were assigned to the low risk and high risk groups according to the median value of the risk scores. Kaplan-Meier survival analyses were conducted to compare the overall survival (OS) and progression-free survival (PFS) in the two groups. The Kaplan-Meier (K-M) method and ROC were performed to evaluate the index efficiency. From 518 patients, 40 patients demonstrated the highest expression level of cell death-related genes, and 40 patients with the lowest expression level of cell death-related genes were picked for further analysis.

Differential analysis of the high and low CDI groups
The differentially expressed genes (DEGs) between the high and low CDI groups were determined using the limma package in R with conditions of adjusted P < 0.05 and |fold change (FC)|> 1 [36]. The volcano plot was constructed by using the ggplot2 package in R.

Clinico-Pathological Analysis and Cox-Proportional Hazard
Pearson's chi-square (χ2) test was conducted to compare categorical variables of clinico-pathological characteristics between groups. Univariable and multivariable Cox proportional hazards models were carried out to assess the performance of the CDI in predicting prognosis. The hazard ratios (HR) with 95% confidence intervals (CI) were based on OS.

Evaluation of Cytokines
To assess the immune activity of GBM patients, cytokine gene list was acquired using the keyword: 'KEGG cytokine-cytokine receptor interaction' (n = 265 genes) (https:// www. gsea-msigdb. org) [37]. The differential expression of cytokines between high and low risk groups of individual cell death pathway and functional enrichment analysis was conducted in the web-based application of STRING ver.11.0 (http:// strin gdb. org) [38].

Estimation of TME (Tumor Immune Environment) cell infiltration
Immune infiltration information, including macrophages, neutrophils, B cells, CD4 + T-cells, CD8 + T-cells, and dendritic cells, etc., were accessed based on the tumor immune estimation resource (TIMER2.0) (http:// timer. cistr ome. org/) [39]. Single-sample geneset enrichment analysis (ssGSEA) was utilized to analyze the subgroups of tumor-infiltrating immune cells between the high CDI and low CDI groups of individual cell death pathway and explore their immune function [40]. At the same time, CIBERSORT [41], xCell [42], MCP-counter [43], quan-TIseq [44] and TIMER [45] algorithms were compared between the two groups, and a heatmap was used to display their differences in the immune response. The correlation between tumor immune cell infiltration and the CDI was analyzed to investigate the performance of CDI in the TME of GBM.

Assessment of the Role of CDI in Immune Checkpoint Blockade (ICB) treatment
Recent researches reported that the expression level of ICB key targets might have a close association with the clinical outcome of ICI [46]. Therefore, the potential immune checkpoints were derived from previous studies [47]. To evaluate the role of CDI in ICB therapy of GBM, we performed correlation analysis between the gene signature and expression level of these ICB key targets.

Functional enrichment analysis
Functional enrichment on gene level was completed by using the g:Profiler program (https:// biit. cs. ut. ee/ gprof iler) [48]. It interprets and maps genes to the corresponding enriched pathways based on well-established data sources. The search tool for the Retrieval of Interacting Genes/Proteins (STRING) database was utilized to conduct the protein-protein interaction (PPI) network to uncover the relationships of DEGs [49]. Cytoscape (Ver 3.8.2) and the plugin of Cyto-Hubba were used for visualizing the PPI network and identifying the top 100 highly connected protein nodes (hubs) by degree, betweenness centrality, and closeness centrality of DEGs [50].

Gene set variation analysis (GSVA)
GSVA was carried out, using the "GSVA" package, to assess the difference in the biological process of the CDI risk groups [51]. Differential analysis of the enrichment scores of KEGG pathways between the high and low risk groups was performed using the limma package in R [35,36,52,53]. The gene sets of "c2.cp.kegg.v7.4.symbols" were retrieved from the MSigDB database for analysis. The molecular pathways enriched differentially between the two groups were determined by |log2FC|> 0.1 and adjusted P < 0.05.

Connectivity Map (CMap) analysis
We used the CMap database (https:// clue. io/) to investigate candidate compounds targeting the molecular pathways and genes associated with CDI for GBM patients [54]. The Connectivity Map (CMap) analysis can also display the mechanism of action (MoA) of compounds. The top 148 most upregulated genes (P < 0.01) and top 148 most downregulated genes (P < 0.01) between the high and low CDI groups were utilized to inquire the CMap database. Potential compounds were identified by the most significant highly expressed genes of each group. The compounds enrichment scores were obtained, and candidate therapeutic drugs for the high risk group were screened out by a negative enrichment score.

Statistical analysis
All the statistical analyses of this study were executed by the R 4.0.4 software, GraphPad Prism (version 7 GraphPad Software), and SPSS 23.0. To compare the clinico-pathological parameters between groups, the independent Student's t test was utilized for continuous data, while the Pearson's chi-square (χ2) test was utilized for categorical data. Statistical differences were compared by the Wilcoxon and Kruskal-Wallis H tests. A two-tailed p-value < 0.05 was considered statistically significant.

Construction of the cell death index
To determine the prognostic signature of each cell death pathway, higher expression of genes in the training set (TCGA) were evaluated for prognostic correlation with OS, and various gene combinations were tested. Finally, the 4-gene apoptosis signature, 8-gene autophagy signature, and 4-gene necrosis signature displayed prognostic association in GBM (Table 1) and the combined 16-gene cell death index (CDI) was generated (Fig. 1A). The interaction between CDI signature genes in GBM is displayed in Fig. 1B.

RNA-Seq Analysis of Patients in TCGA
Seven hundred sixty-nine differentially expressed genes were identified and 621 of them were upregulated at > onefold in high CDI group compared to the low CDI group (Table S4). In the low CDI group, 148 genes were upregulated at > onefold compared to the high CDI group (Table S5). The volcano plot of differentially expressed genes between two groups is illustrated in Fig. 2A.

High CDI patients possessed a higher CNA burden and lower TMB
CNA and somatic mutation analyses were carried out to investigate the genomic variations in the two CDI risk groups. The differential analysis of CNA between the two groups showed that, compared with the low CDI group, 225 (35.3%) genes were significantly amplified, and 309 (48.4%) genes had deletion variation in the high CDI group (Fig. 2B). Moreover, the boxplots exhibited more copy number amplifications(p = 3.6e-05) and deletions(p = 0.042) burdens in the high CDI group (Fig. 2C). Somatic mutation analysis revealed that 31 out of the 769 DEGs (4.0%) had a mutation frequency > 1%, and most (83.9%, 26/31) of them were upregulated in the high CDI group. The top 15 most frequently mutated DEGs altered in the 126 GBM samples were visualized by oncoplot and illustrated in Fig. 2D. Somatic mutation analysis of two CDI groups separately uncovered that each risk group had distinct top mutated genes (Fig. 2E).
In the low CDI group, TP53 (41%) was most frequently mutated, higher than that in the high CDI group, in which PTEN (33%) was most frequently mutated ( Fig. 2F). Figure 2G demonstrated that the CDI risk score was significantly higher in IDH-wide-type samples than in IDH-mutant samples (P < 0.001) and patients with MGMT promoter methylation (P = 0.002) showed significantly lower CDI risk scores. All these findings might uncover the underlying differences in response to immunotherapy between the two groups.

Associations between CDI and clinical features
We calculated the CDI of 518 GBM patients and ranked them from low to high to analyze the associations between the CDI and clinical features. The demographics and follow-up data of the patients in the high CDI and low CDI groups were compared and presented in

Fig. 1 A)
Combined cell death index (CDI) was generated, which included the highest expression of genes involved in autophagy, apoptosis, and necrosis. B) The interaction between CDI signature genes in GBM. The circle size represented the effect of each signature gene on the prognosis, and the range of values calculated by Log-rank test was p < 0.001, p < 0.01, p < 0.05 and P < 0.1, respectively. The Autophagy, apoptosis and necrosis signature gene was marked with blue, yellow and red respectively. Green dots in the circle represent protective factors of prognosis and black dots in the circle represent risk factors of prognosis. The lines linking signature genes showed their interactions, and thickness showed the correlation strength between genes. Negative correlation was marked with blue and positive correlation with red 95%CI:1.541-2.858, P < 0.001) ( Table S6), suggesting that the CDI high risk could independently predict the OS and PFS of GBM patients. In addition, radiotherapy was a favorable factor for patients' survival (OS, HR = 0.224, 95% CI:0.125-0.403, P < 0.001; PFS, HR = 0.420, 95% CI:0.299-0.589, P < 0.001).
As illustrated in Fig. 3, patients were assigned into either the high CDI group or the low CDI group according to their median risk score (Fig. 3A). With the risk score increasing, the number of alive patients reduced (Fig. 3B). Kaplan-Meier analyses exhibited that patients with a high CDI had significantly worse OS and PFS than those with a low CDI (Fig. 3C, P < 0.0001; Fig. 3E, P < 0.0001). The predicted efficiencies of CDI were evaluated by ROC curves (Fig. 3D and Fig. 3F). The results illustrated that the AUCs of CDI for predicting the 1.5-, 3-and 4.5-year OS were 0.727, 0.833 and 0.844, respectively.

Validation of CDI
One hundred thirty-seven GBM patients in the CGGA database were enrolled as validation data set for assessing the signature performance. The CDI risk score of the validation cohort patients was then calculated according to the risk score formula and subsequently allocated into either the high CDI or the low CDI group (Fig. 4A). The demographics and clinicopathological data of the patients in the two groups were compared and presented in Table 2. Univariable and multivariable Cox analyses found that the high CDI group was significantly associated with OS (HR = 1.498, 95%CI:1.037-2.166, P < 0.05) ( Table 3). Consistent with the findings in TCGA, the Kaplan-Meier survival curves demonstrated that high CDI patients had a poorer prognosis than those with a low CDI (Fig. 4C, P < 0.0001). The AUCs of ROC curves for predicting the 1.5-, 3.0-and 4.5-year survival of GBM in the dataset were 0.607, 0.600, and 0.721, respectively (Fig. 4D). Concomitantly, the OS and alive patients deceased with an increase in risk scores (Fig. 4B). These results implied a satisfactory performance of CDI for survival prediction in GBM patients.

Cytokine Gene Expression Analysis
In the autophagy group, TNFRSF12A, OSMR, LIF, CCL2 and VEGFA displayed higher expression levels in the high risk group, whereas BMP2 showed a higher expression level in the low risk group (Fig. 5B). Analyses conducted on the apoptosis and necrosis groups were depicted in Fig. 5A and 5C. In the CDI group, the expression level of CCL2, LIF, FAS, IL1B, CXCL10, CCL20, TNFSF13, IL6, IL1R2, IL10RA, and IL13RA1, among others, were higher in high CDI patients whereas the BMP2 expression level was higher in the low CDI patients (Fig. 5D). Furthermore, functional enrichment analysis showed that the high CDI group seemed to have a higher inflammatory cytokines proportion than the low CDI group (Fig. 5E). Figure 6A displays the heatmap of immune responses between low CDI and high CDI groups, based on CIB-ERSORT, xCell, MCP-counter, quanTIseq, ssGSEA and TIMER algorithms. The results established that the high CDI group exhibited higher immune scores, indicating a significantly increased immune cell infiltration. To estimate the correlations between the CDI and immune cells, spearman analysis was performed, and the results are illustrated in Fig. 6B. Correlation analysis based on the ssGSEA of the TCGA dataset showed that higher Treg, T helper cell and macrophages infiltration correlated with the high risk group of all the three cell death pathways (Fig. 6C, 6D and 6E). Notably, neutrophils and aDCs infiltration were higher in patients with high necrosis, apoptosis, and CDI groups. In contrast, B cells were enriched in patients with low necrosis, apoptosis, and CDI. Furthermore, the ssGSEA algorithm demonstrated a significant difference between the low CDI and high CDI groups in T cell functions, including checkpoint (inhibition), cytolysis, CCR, regulation of inflammation, co-stimulation of T cell, co-inhibition of T cell, and type I INF response, as illustrated in Fig. 6F. Altogether, these results imply that CDI may present a novel understanding of the immune response in GBM.

CDI correlated with Key Genes of ICB Therapy in GBM
Correlation analysis showed that the expression levels of 20 (i.e., PDCD1, IDO1, etc.) ICB-related genes were significantly different between high CDI and low CDI groups (Fig. 7D). Similar analyses were also performed Fig. 2 A) Volcano plot showing the differential expression of genes between high CDI and low CDI patients (p < 0.05, |log2 fold-change|> 1). B) The differential analysis of copy number variations between two groups was visualized by Circos plot, which revealed that compared with the low CDI group, 225 (35.3%) genes were significantly amplified, and 309 (48.4%) were significantly deleted in the high CDI group.  in the other three classically known cell death pathways (Fig. 7A, 7B and 7C). We singled out seven key ICI genes (CD86, CD40, PDCD1, CD48, CD200, IDO1, and PDCD1LG2) for further research. To investigate the potential effect of the CDI in the ICB therapy of GBM, the correlation analyses were performed between the expression level of the ICB key genes and the CDI risk scores (Fig. 7E). The result revealed that the CDI risk score had close relationship with CD86 (r = 0.28; p = 7.0e − 11), PDCD1LG2 (r = 0.23; p = 1.4e − 07), CD48 (r = 0.21; p = 1.3e − 06), and PDCD1(r = -0.23; p = 1.4e − 07) (Fig. 7F), indicating CDI might have a nonnegligible role in the outcome prediction of ICB treatment in GBM.

Enrichment analysis
Functional enrichment analysis of differential gene expressions between high CDI and low CDI groups identified 769 genes based on an inclusion filter of > onefold for the DEGs. 621 genes of them were upregulated in high CDI group, while 148 genes of them were upregulated in the low CDI group (Fig. 8A and 8B). Patients in high CDI group had significant immune-related pathways, whereas the patients in low CDI group lacked enrichment in immune-related pathways. The enriched Gene Ontology terms in high CDI group were receptor ligand activity, signaling, cytokine, and chemokine activity ( Fig. 8C and 8D). Conversely, the enriched Gene Ontology terms in low CDI group were predominated by transmembrane transporters and gated channel activity ( Fig. 8E and 8F). Potential molecular pathways and underlying mechanisms related to the CDI of GBM were investigated by performing GSVA. The results identified that 19 pathways were positively correlated with the high CDI group, while 1 pathway was positively correlated with the low CDI group (Fig. 8G). High CDI patients mainly correlated with toll-like receptors signaling, JAK-STAT signaling, chemokine signaling pathway, complement and coagulation cascades, ECM receptor interaction, and others, whereas the low CDI group mainly correlated with glyoxylate and dicarboxylate metabolism.
The PPI network of the 769 overlapping DEGs was selected from the STRING database. The top 100 hub proteins were mined using Cyto-Hubba, and the PPI network was then performed and visualized in the Cytoscape software (Fig. 8H). The network suggested strong interactions among these hub genes.

Potential compounds targeting the two CDI groups
CMap analysis was performed to investigate potential drugs targeting the two CDI groups of GBM (Table S7,  Table S8, and Figure S1). Compounds with a score less than -90.0 or higher than 90.0 were selected, the predicted drugs are listed in Table 4. It was determined that nine compounds targeting 7 molecular pathways had a score of less than -90 and were potential inhibitors for high CDI patients. On the other hand, 5 pathways targeted by 7 compounds with a score higher than 90 were considered potential inhibitors for low CDI patients.

Discussion
GBM is the most common primary intracranial tumor with high malignancy and poor outcomes. Relevant molecular biomarker models for predicting the prognosis of GBM, such as hypoxia, ferroptosis-related gene, and stem cell signatures, have been established to improve patients' survival and develop novel individualized treatments [55][56][57]. However, their accuracy and predictive abilities remain limited, and most of them were developed based on single transcriptomics data with inadequate focus on biological pathways. Autophagy, apoptosis  and necrosis represent three significant biological hallmarks of tumors that are valuable in predicting GBM patients' prognosis. In this study, we constructed a combined gene signature involved in autophagy, apoptosis, and necrosis, and evaluated its prognostic effect and correlation with immune cells and mediators.
We first integrated the genomic information of 518 GBM patients from the TCGA database to comprehensively evaluate the prognostic of genes involved in the three classically known cell death pathways and identified a novel 16-gene CDI signature. Among these significant genes of the signature, BID, CHP1, MAPK3, PRKAB2, PRKAG2 and TRAP1, were favorable factors for GBM patients' survival in this study. BID is a proapoptotic member of the Bcl-2 protein family and could promote Ca2 + -induced neuronal injury [58]. Inhibition of BID during acute endoplasmic reticulum stress may protect against cell death [59]. In addition, studies reported that BID was associated with the prognosis of immunoglobulin A nephropathy [60]. TRAP1 which affects the mitochondrial protein quality control system and mitochondrial metabolism, has multiple functions in mitochondria. Researchers found that suppressing the function of TRAP1 would benefit temozolomide therapy in GBM in vitro, revealed great potential and practical value in GBM treatment. [61]. CHP1, MAPK3, PRKAB2 and PRKAG2, these four genes, were rarely focused on. Other ten genes of the CDI gene signature, including CFLAR, PRKAR1B, EEF1A2, LAMTOR3, MET, SER-PINA1, SREBF1, CASP3, NOL3, TRAF3, were risk factors for GBM patients' survival. EEF1A2 is a translation factor selectively expressed by heart, skeletal muscle, nervous system and some specialized cells. It is a putative oncogene highly expressed in ovarian cancer [62]. Studies have shown its negative prognostic role in breast cancer, nonsmall cell lung cancer and gastric cancer [63][64][65][66]. MET plays a well-defined role as a selectable oncogenic driver of tumor proliferation. Preclinical studies found that inhibition of MET could reduce cell survival, local invasion and metastasis to distant sites [67]. CASP3 plays an important role in the development of the brain. Moderate active CASP3 levels were found in human GBM samples and down regulation of CASP3 may inhibit the migration GBM cells, suggesting that CASP3 inhibition may serve as a novel therapeutic strategy for GBM [68]. NOL3 protects against oxidative stress-induced cell death. Researches showed a tumor suppressor role of NOL3 in the pathogenesis of myeloid malignancies [69]. LAMTOR3 and SERPINA1 both showed prognostic effect in tumors [70,71]. Compared with the previous autophagy-related gene signatures for GBM, the 8 autophagy-related genes included in CDI were all novel [72][73][74]. In addition, their data source and bioinformatic processes are different. Our findings could propose several potential targets for glioblastoma therapy. As the first prognostic model with three cell death pathway-related genes for GBM patients, the CDI signature can simultaneously reflect the tumor's cell death characteristics and is convenient for clinical use.
We classified patients into two groups based on their CDI, and analyzed the correlation between CDI subtypes and clinical features. The results indicated that low CDI patients had longer OS and PFS. Moreover, the correlation between CDI and clinicopathological parameters, as well as survival outcomes, suggests that CDI has an effective prognostic prediction in GBM patients. Afterward, the CDI signature was verified in the CGGA cohort patients. ROC curves and survival analysis revealed an efficient performance of CDI.
Further in this study, we explored the correlation of CDI with the TME cell infiltration and ICI in the prognosis of GBM. GBM is one of the most immunologically "cold" tumors, which presents an unsatisfactory treatment effect with immunotherapy [75]. Although the results so far are mostly disappointing, a large number of clinical trials indicated that immunotherapy remains conceptually promising for GBM treatment [76][77][78]. Previous studies revealed that the poor results of immunotherapy were mainly due to the profound immunosuppressive characteristic of GBM and a TME that is challenging for immune cells [26]. Profiling the tumor immune microenvironment has lately been regarded as one of the future breakthroughs to improve immunotherapy of GBM [79,80]. Herein, neutrophils, macrophages, T helper cells, Tregs, TIL, and aDCs were enriched in patients with a high CDI. The CDI risk score was significantly correlated with immune cell infiltration (i.e., neutrophils, macrophages, T helper cells, Tregs, TIL, and aDCs). Besides, the ssGSEA analysis also pointed that the infiltrating immune cells (i.e., neutrophils, macrophages, T helper cells, Tregs, TIL, and aDCs) were remarkably increased, and immune signatures (i.e., checkpoint, cytolytic, CCR, regulation of inflammation, co-stimulation of T cell, co-inhibition of T cell, and type I INF response) were significantly activated when the risk score was elevated, indicating that the CDI signature possesses an unneglectable role in the TME of GBM. The higher neutrophils and macrophages infiltration in high CDI risk patients reflected high inflammation and local immune dysfunction in TME [81]. In fact, the presence of macrophages has been approved to be a negative predictor for survival in high-grade gliomas murine models and it could impair antitumor immunological functions of TME in glioblastoma [82,83]. Effector T cell infiltration in the tumor has been identified positively correlated with patients' survival whereas the Treg fraction increase in GBM patients implied a deficit of patient effector T-cell responsiveness, suggesting that higher Tregs might be related with poor survival of patients [84,85]. All the above studies suggest correlations between the CDI and TME. Therefore, targeting specific GBM patients based on CDI selection may help immunotherapy achieve a superior therapeutic effect.
With the development of ICB treatment, ICIs have greatly affected the cancer treatment scheme [86,87]. Indeed, ICB treatment has brought a novel field for GBM patients' clinical management [88,89]. Even though current clinical trials of ICB treatment in GBM failed, demonstrating appropriate patients, benefiting from ICB treatment, through various biomarkers, still appears promising in GBM treatment. In this study, the CDI signature was negatively associated with the ICB treatment key target genes(i.e. PDCD1 and CD200), and the expression level of some ICB-related genes ( i.e. PDCD1LG2,CD86, CD48, CD40, and IDO1) elevated significantly with an increase of the CDI risk score. The activation of IDO is involved in tumorigenesis by helping tumor cells evade immune surveillance [90]. In addition, studies have found that high IDO1 expression was associated with the poor prognosis of GBM [91]. Therefore, ICB therapy targeting patients with high IDO1 might achieve superior treatment effect. Besides, our data also revealed that low CDI risk patients were more likely to have a higher tumor mutation burden. Researchers construed that the low TMB in GBM was a fundamental cause for its failure in immunotherapy except for the immunosuppressive microenvironment [92,93]. Thus, low CDI patients with high TMB were speculated to show a satisfactory therapeutic effect with ICB treatment. Overall, our findings indicated that CDI might well predict the clinical outcome of ICB therapy in GBM patients.
CMap analysis was performed to find potential compounds targeting genes related to CDI, enabling more practical results. In the case of the high CDI group, ruxolitinib was uncovered as one of the candidate drugs. Ruxolitinib was recently found to prevent glioblastoma invasion and tumorigenesis by inhibiting the IFNinduced JAK/STAT signaling pathway [94,95]. Besides, it has also been proposed for use in patients with chronic neutrophilic leukemia for its safety and efficacy in inhibiting JAK1/2 [96]. For the low CDI group, 3 out of 7 compounds with a score higher than 90 were PKC activators, consistent with a recent study that confirmed PKC as a suitable druggable target to treat recurrent GBM [97]. The above researches demonstrated the reliability of drug screening in our study and the feasibility to applicate in GBM treatment. To further authenticate the therapeutic value of these compounds, more studies are warranted.
Although the results of our study indicate that CDI can serve as an effective prognostic signature for GBM patients, there are also some limitations. First, our research is mainly based on integrating the genomic data from public datasets to comprehensively assess the role of CDI with the TME cell-infiltrating characteristics in GBM  Oligodendroglioma, IDH-mutant and 1p/19q-codeleted; and Glioblastoma, IDH-wildtype [98]. The classification of glioblastoma, IDH wild type includes the category of glioblastoma, IDH wild type defined in the 2016 World Health Organization Classification of Tumors of the Central Nervous System. In addition, diffuse astrocytoma, if accompanied by microvascular proliferation or necrosis or TERT promoter mutation or EGFR gene amplification or + 7/-10 chromosome copy number changes, will also be included in glioblastoma, IDH wild type category. Our research data of GBM patients were derived from TCGA database and the Chinese Glioma Genome Atlas (CGGA) database in which patients were classified according to the 2007/2016 WHO classification system [99,100]. It is really hard to reclassify patients obtained from the previous data with the new WHO classification, due to the complicated new diagnostic criteria with various mutation information and potential diagnostic bias in the second classification. Therefore, analysis of CDI in GBM patients classified by the latest WHO Classification is warranted subsequently with new constructed cohort. Third, an immunotherapy cohort is needed to validate the relationship between CDI and immunotherapy response for the limited number of patients under immunotherapy.
In conclusion, this research provides a potential approach for screening patients at higher risks of mortality based on the cell-death-based gene signature. Moreover, CDI was significantly associated with immune cell infiltration as well as ICB treatment key genes in GBM. Thus, this study brings a novel signature to promote the individualized prediction of overall survival and deepen the understanding of the TME in GBM, further providing promising clinical applications in GBM therapy.

Funding
None.