Developing an Immune-Related Signature for Predicting Survival Rate and the Response to Immune Checkpoint Inhibitors in Patients With Glioma

Background: Glioma is one of the most aggressive cancer types affecting the central nerve system, with poor overall survival (OS) rates. The present study aimed to construct a novel immune-related signature to predict prognosis and the efficiency of immunotherapy in patients with glioma. Methods: The mRNA expression data and other clinical information of patients with glioblastoma multiforme (GBM) and low grade glioma (LGG) were obtained from The Cancer Genome Atlas and Chinese Glioma Genome Atlas databases. The immune-related genes were obtained from the Immunology Database and Analysis Portal database. Subsequently, an immune-related signature was created following the results obtained from the Least Absolute Shrinkage and Selection Operator regression model. To validate the predictability of the signature, Kaplan-Meier survival curves and time-dependent receiver operating characteristic curves were created. Moreover, both univariate and multivariate analyses were carried out using the OS between this signature and other clinicopathologic factors, and a nomogram was constructed. In addition, the association between signature, immune cell infiltration, tumor mutation burden and immunophenoscore were determined. Results: Results of the present study using 118 GBM and LGG samples uncovered 15 immune-related genes that were also differently expressed in glioma samples. These were subsequently used to construct the immune-related signature. This signature exhibits the ability to predict prognosis, the infiltration of immune cells in the tumor microenvironment and the response of patients with glioma to immunotherapy. Conclusion: Results of the present study demonstrated that the aforementioned novel immune-related signature may accurately predict prognosis and the response of patients with glioma to immunotherapy.


INTRODUCTION
Glioma is one of the most aggressive and common malignant cancer types affecting the central nerve system (CNS) (Chen et al., 2021a), , which poses a serious threat to human health worldwide. In elderly patients (age, >65 years), the 5-year survival rate is~5% (Hendriks et al., 2019). According to the 2021 World Health Organization classification standard, gliomas are divided into grades, namely I-IV (Louis et al., 2021). Notably, grade I glioma exhibits slow proliferation and an improved prognosis, and grade IV glioma (also known as glioblastoma multiforme; GBM) is highly proliferative and invasive, exhibiting rapid recurrence and poor rates of overall survival (OS) and progression-free survival (PFS) (Soeda et al., 2015;Karim et al., 2016). Despite the development of GBM therapy over numerous years, the prognosis of patients with GBM remains poor (Janjua et al., 2021). In patients with GBM, the OS and PFS are~14 and 7 months, respectively (Ostrom et al., 2014;Kelly et al., 2017). The standard therapy used for the treatment of GBM is surgery combined with radiotherapy and chemotherapy (Di Cintio et al., 2020). However, the current therapy model cannot prevent the recurrence of GBM effectively. As a result of these factors, a novel model to predict prognosis and novel therapy options for patients with glioma are required.
Following technological advances and progress made in immunotherapies, advances have been made in the treatment of multiple cancer types, including brain metastasis (Tawbi et al., 2018;Hendriks et al., 2019). Notably, there have been two main obstacles preventing the development of effective immunotherapies for GBM: 1) Development of a highly immunosuppressive microenvironment and 2) high tumor heterogeneity (Lim et al., 2018). However, a number of recent studies have suggested that anti-tumor responses to immunotherapy may also occur in the micro-environment in the brain, which may be promising for the development of novel treatment options for malignant gliomas Xun et al., 2021). Three types of immune checkpoint inhibitors (ICIs) have been tested to treat GBM; namely, anti-PD-1 antibody, anti-CTLA-4 antibody and anti-LAG-3 antibody (Preusser et al., 2015). Although the results of previous studies have suggested that the response rates of ICIs for glioma were not optimal (Berghoff et al., 2015;Kurz et al., 2018), these developments are useful to determine the reasons for therapy failure and to distinguish which patients with glioma may benefit from the current immunotherapy options. The tumor immune microenvironment (TIME) and immune-related genes (IRGs) play important roles in immunotherapy (Remon et al., 2021), and can be used as biomarkers for predicting the survival rates of patients (Pan et al., 2020) {Wan, 2021 #6}. However, few studies have previously explored whether IRGs and TIME may act as biomarkers for survival rate prediction and determining the response to immunotherapy in patients with glioma.
The present study aimed to create a novel signature for predicting both the prognosis and response to immunotherapy in patients with glioma, according to IRGs and TIME data obtained from public databases. Following the development of this signature, the association between clinicopathological factors and prognosis in glioma was also detected. Furthermore, the association between immune cell infiltration, tumor mutation burden (TMB), immunophenoscore (IPS) and this signature was determined in patients with glioma. Results of the present study demonstrated that the signature may aid in improved prognosis prediction, and may provide a novel theoretical basis for further elucidating immunotherapy options for glioma.

Expression and Clinical Information of Patients
Clinical characteristics and mRNA expression levels of GBM and LGG samples were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) and the Chinese Glioma Genome Atlas (CGGA) database (http://www. cgga.org.cn/). After matching the expression information with the clinical information, a total of 1,008 samples were used in the present study, which included 1,003 tumor and 5 healthy samples.

Differentially Expressed Immune-Related Genes (DEIRGs)
A total of 1,793 IRGs were downloaded from the Immunology Database and Analysis Portal (ImmPort) database (https://www. immport.org/home). Using the R Package "limma", a total of 2,321 differently expressed genes (DEGs) were determined using the transcription information of 690 samples, which followed the thresholds (|log2(Fold Change)|>1 and false discovery rate (FDR) < 0.05. A total of 158 genes were extracted as DEIRGs, both existing in IRGs (1,793 genes) and DEGs (2321 genes); this process was visualized using a Venn diagram ( Figure 1B). In the DEIRGs group, the expression levels of 80 genes were upregulated and 78 genes were downregulated. The volcano plot of DEIRGs was created using the R package "ggplot2".

Functional Enrichment Analyses
Two types of functional enrichment analyses, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG), were carried out to determine the biological functions of DEIRGs, using the online tool Database for Annotation, Visualization, and Integrated Discovery (DAVID) (version, 6.8; https://david. ncifcrf.gov/). The enrichment of GO terms and KEGG signaling pathways followed the threshold (FDR<0.05). Whole GO terms and the top 20 most significant KEGG pathways were visualized using the R package "ggplot2" (Figures 1C,D).

Development and Accuracy Evaluation of the Immune-Related Signature
The 690 tumor samples obtained from the TCGA database were allocated into two sets randomly, named the TCGA Training Set (n = 345) and the TCGA Validation Set (n = 345). Moreover, the 313 tumor samples obtained from the CGGA database were named the CGGA Set (n = 313). All sample details are displayed in Table 1. The TCGA Training Set was used to construct the risk signature. TCGA Testing Set, Entire TCGA set and CGGA Set were applied to verify the prognostic accuracy. A univariate Cox proportional hazard regression analysis was carried out to select a total of 118 DEIRGs from the whole DEIRGs cohort, which associated with the OS of patients with glioma (p < 0.05; Supplementary Table S1). A least absolute shrinkage and selection operator (LASSO) test was performed on the 118 DEIRGs to minimize the error rate using the "glmnet" R package, and 15 hub genes were used to construct the prognostic signature of patients with glioma (Supplementary Figure  S1). The formula of risk score was as followed: Risk score = [Expression level of Gene 1×coefcient]+[Expression level of Gene 2×coefcient]+. . .+[Expression level of Gene n × coefficient]. The risk scores from patients were calculated using the aforementioned formula. Patients were subsequently allocated into high-risk and low-risk groups by the median of risk scores. To verify the accuracy of this 15-gene risk signature, the K-M survival curves and receiver operating characteristic (ROC) curves were used to compare the survival rates of patients in the low-risk and high-risk groups using the "survival" and "survivalROC" R package.

Evaluation of Immune Cell Infiltration
Details of immune cell infiltration were obtained from the TCGA RNA-sequencing database and were calculated using the online tool CIBERSORT (https://cibersort.stanford.edu/).

Analyses of IPS and TMB
IPS data were obtained from The Cancer Immunome Atlas (TCIA) (https://tcia.at/home). TMB data were obtained from the TCGA database (https://portal.gdc.cancer.gov/) and stored in Mutation Annotation Format (MAF). TMB analyses were carried out using the R package "maftools" and "TCGAmutations".

Statistical Analysis
Univariate and multivariate Cox regression calculations were applied using the R package "survival" to validate the predictability of the signature and clinicopathological factors. Differences between clinical features were determined using the independent t-test. p < 0.05 was considered to indicate a statistically significant difference.

Construction of the Immune-Related Signature
According to the criteria (FDR<0.05 and |log2 (fold change)|>1), a total of 2,163 DEGs (811 upregulated and 1,352 downregulated genes) were identified and used in subsequent analyses. A Venn diagram was created to extract 158 DEIRGs (80 upregulated and 78 downregulated genes) from the DEGs cohort ( Figures 1A,B). Results of the KEGG functional enrichment analysis demonstrated that DEIRGs exhibited a high significant correlation with the "Axon guidance" signaling pathway ( Figure 1C). Moreover, results of the GO enrichment analysis suggested that "cellular process", "binding" and "cell" were the most enriched terms of "BP", "MF" and "CC", respectively ( Figure 1D).

Immune-Related Signature and the Survival of Patients With Glioma
Patients were allocated into high-risk and low-risk groups according to the median of risk scores (Figure 2A). Results of the present study demonstrated that patients in the low-risk group exhibited a prolonged OS, compared with those in the high-risk group (p < 0.05; Figure 2C). Moreover, the accuracy of the immune-related signature was tested using time-dependent ROC curves ( Figure 2B).
In the TCGA training set, the under areas of 1-year, 3-year and 5-year survival were 0.88, 0.94 and 0.93, respectively, which demonstrated the ability of the immune-related signature to predict the prognosis of patients with glioma. The same ROC curves were also drawn using the data in the TCGA Validation set, TCGA Entire set and CGGA set. All results demonstrated a high predictive ability of the immunerelated signature. The analysis of the correlations between risk score and survival rate of samples in GBM cohorts or LGG cohorts, which showed the same results (Supplementary Figure S2).

Immune-Related Signature as an Independent Prognostic Factor
To verify whether the immune-related signature may act as an independent factor to predict the survival of patients with glioma, both univariate and multivariate analyses between risk score, OS and other clinical indicators were carried out ( Table 3). Results of the multivariate analysis using the Entire TCGA set and CGGA set demonstrated that risk score, age, histology, tumor grade and 1p19q-co-deletion-status may act as independent prognostic indicators for patients with glioma (p < 0.05). As there is no TNM stage in patients with glioma, the category of TNM stage was omitted.

Immune-Related Signature and Clinicopathological Factors
To determine the association between IRDEGs prognostic risk score and the other clinicopathological factors, independent t-tests were performed between risk score, gender, histology, tumor grade, 1p19q-co-deletion, MGMTp-methylation and IDH mutation. Risk scores were notably increased in advanced grade tumors (such as  glioblastoma), deleted-1p19q, methylated-MGMTp and IDHwildtype ( Figures 3B-F). However, there was no significant difference between males and females for risk scores ( Figure 3A). Moreover, a nomogram was created according to risk score and the aforementioned clinicopathological indicators, for constructing an accurate tool to predict the survival of patients with glioma ( Figure 3G). Notably, the calibration curves of the prognostic nomogram suggested a positive consistency between predicted and exact 1-year, 3-year and 5-year survival ( Figure 3H).

Immune-Related Signature and Immune Cell Infiltration
To determine functions of the immune-related signature in the TIME, CIBERSORT was used to estimate immune cell infiltration, and to estimate the relative proportion of 22 types of immune cells in glioma. Immune cell type abundance in the high-risk group and low-risk group are displayed in Table 4. As shown in Figure 4 and Table 4, 5 types of immune cells (T cells CD8, T cells CD4 memory resting, Macrophages M0, Macrophages M2 and Dendritic cells resting) exhibited a positive correlation with risk score (p < 0.05; Figure 4; Table 4). In total, 7 types of immune cells (Plasma cells, T cells CD4 memory activated, NK cells resting, NK cells activated, Macrophages M1, Eosinophils and Neutrophils) were negatively correlated with IRDEGs prognostic risk score (p < 0.05; Figure 4; Table 4). Moreover, the samples were divided into high and low groups according to the median of different types of immune cell infiltrations, and the K-M survival curves were used to determine the association between the OS and immune cell infiltration. As the risk score demonstrated a negative correlation with glioma prognosis, it was hypothesized that the association between immune cell infiltration and risk score, and between the immune cell infiltration and survival would demonstrate the opposite trends. In total, 8 types of immune cells (T cells CD8, NK cells resting, NK cells activated, Macrophages M0, Macrophages M1, Macrophages M2, Dendritic cells resting and Neutrophils) proved this hypothesis. The K-M survival curves of these are displayed in Figure 5.

Immune-Related Signature and TMB
Results of previous studies have demonstrated that TMB acts as a biomarker of prognosis and immunotherapy (Halbert and FIGURE 4 | Association between immune cell infiltration and the immune-related signature in glioma. The blue columns represent the low-risk groups. The red columns represent the high-risk groups. Einstein, 2021). Thus, in the present study, the association between immune-related signature and TMB were determined, in order to explain the predictive capability of the immunerelated signature. The top 20 most frequently mutated genes in the low-risk and high-risk groups are displayed in Figure 6A and Figure 6B, respectively. Moreover, TMB in the high-risk group was compared with that in the low-risk group. Results of the present study demonstrated that TMB exhibited a positive correlation with risk score (p < 0.001; Figure 6C). The samples were further allocated into high-TMB and low-TMB groups, and the OS was compared between these two groups (p < 0.001; Figure 6D). Notably, OS in the high-TMB group was decreased compared with that in the low-TMB group.

Immune-Related Signature and Immune Checkpoint Inhibitors
The present study utilized four types of score (IPS, IPS-CTLA4, IPS-PD1-PDL1-PDL2 and IPS-CTLA4-PD1-PDL1-PDL2) to evaluate the potential of patients with glioma to be placed on ICIs. Firstly, we evaluated the association between these four types of score and OS in glioma. The results showed these four scores significantly had the negative correlation with the OS (p < 0.001, Figure 7A). Then these four scores were significantly elevated in the high-risk group compared with those in low-risk group (p < 0.001; Figure 7B), which explained the prognostic ability of our signature. Moreover, the same analyses were performed on the expression levels of CTLA4, PD1, PDL1 and PDL2, which were significantly related with the OS of glioma patients and increased in the high-risk group (p < 0.001; Figures 7C,D).

DISCUSSION
Due to the high level of malignancy and poor prognosis of glioma (Xun et al., 2021), the development of accurate methods is required to predict the associated prognosis and to create effective therapy for patients with glioma. At present, immunotherapy demonstrates high levels of effectiveness in the treatment of malignant tumors (Chen et al., 2021b;Sprooten et al., 2021). Results of previous studies have demonstrated that TIME exhibits a high association with the occurrence and development of several tumor types, and this may act as an important prognostic indicator (Berraondo et al., 2016;Elola et al., 2018). With the development of bioinformatics, many signatures have been established, which demonstrate a high capability to predict the prognosis of patients with tumors (Hermansen et al., 2017;Kang et al., 2020). Therefore, a novel and more accurate immune-related signature was created in the present study to predict the prognosis and response to immunotherapy in patients with glioma. However, the efficiency of immunotherapy in glioma is not significant at present, indicating that a low number of patients can benefit from it (Touat et al., 2020). The reasons for immunotherapy failure are likely to be complex, but appear to be associated with a lack of biomarkers guiding ICB in different patients (Topalian et al., 2015;Lim et al., 2018). Therefore, a novel signature was constructed to predict the survival and response to ICIs in glioma.
In the present study, the immune-related signature was constructed using GBM and LGG samples from the TCGA database ( Figure 8). In total, 118 genes named DEIRGs were identified to be associated with immune processes and were differently expressed between glioma and healthy tissues. Following the GO and KEGG enrichment analyses, the association between DEIRGs and immune processes were further determined, such as "Antigen processing and presentation" and "T cell receptor signaling pathway". After performing the Cox analysis and LASSO analysis on 118 DEIRGs, a total of 15 hub genes (MET, SSTR2, CDK4, S100A16, TNFRSF10B, BIRC5, BMP2, ADCYAP1R1, BMP1, APOBEC3C, ANGPTL2, TNFRSF19, VGF, NRG3, and JAG2) were highlighted. Among them, the expression levels of 5 genes (MET, SSTR2, VGF, NRG3 and JAG2) were downregulated, and 10 genes (CDK4, S100A16, TNFRSF10B, BIRC5, BMP2, ADCYAP1R1, BMP1, APOBEC3C, ANGPTL2 and TNFRSF19) were upregulated. Using these 15 hub genes, the immune-related signature was constructed and the risk score of each sample was determined. Samples were allocated into highrisk and low-risk groups according to the median of the risk scores. Moreover, univariate and multivariate analyses were performed using risk score and other clinical indicators.
Results of the present study demonstrated that risk score may act as an independent indicator to predict the prognosis of patients with glioma. To verify the predictive accuracy of the constructed signature, the association between risk score and clinicopathological factors were determined, and a nomogram was constructed by integrating risk score and these clinicopathological factors.
To determine the mechanism underlying the predictability of the immune-related signature in the present study, the association between 22 immune cells, risk scores and the prognosis of patients with glioma was determined using CIBERSORT. Results of the present study demonstrated that T cells CD8, Macrophages M0, Macrophages M2 and Dendritic cells resting were negatively associated with the OS of patients with glioma, and NK cells resting, NK cells activated, Macrophages M1 and Neutrophils demonstrated the opposite of these results. Notably, CD8 + T cells have previously demonstrated an improved prognosis in the majority of tumor types, which was not demonstrated in the present study. However, this may be due to the immune environment in the CNS. In addition, the association between TMB, risk score and prognosis was determined in patients with glioma. Results of the present study demonstrated that risk score exhibited a positive correlation with TMB, and TMB was negatively associated with the OS of patients with glioma. This may also indicate why the immune-related signature constructed in the present study can predict the prognosis of patients with glioma.
To date, immunotherapy has been used to treat a number of different tumor types, such as lung cancer (Tang et al., 2022), melanoma (Baruch et al., 2021), cervical cancer (Ferrall et al., 2021), liver cancer (Feng et al., 2021). Some immune checkpoints, such as PD-1 and CTLA-4 may lead to anti-tumor immunity (Lee and Seong, 2021). Results of a previous study demonstrated that nivolumab and pembrolizumab, two therapeutics against PD-L1/PD1, have been approved for subsequent line therapy (Leone et al., 2021). However, according to the Phase III clinical trial, the efficiency of immunotherapy in glioma is not significant (Wu and Lim, 2021). These results may still indicate which patients with glioma are suitable for undertaking immunotherapy. Moreover, results of previous studies have demonstrated that TIME and TMB are associated with prognosis in many tumor types, and this is associated with the efficiency of immunotherapy (Cheng et al., 2021;Wei et al., 2021). These results led to the development of the signature in the present study. IPS is an accurate method to predict the response to ICIs, and this has been widely used to guide immunotherapy in numerous tumor types (Meng et al., 2021;Xu et al., 2021). Thus, the association between IPS and risk score was determined in the present study. The results suggested that patients in the highrisk group exhibited increased levels of IPS, which meant that they may exhibit a poor prognosis, but their response to immunotherapy may be improved. Finally, the association between risk score and the other 6 ICBs (LAG3, TIGIT, TIM3, CD27, ICOS, and IDO1), which were predicted by CIBERSORT, was determined. Results of the present study demonstrated an improved immune response compared with traditional PDL1/CTLA4 immune checkpoints in patients with glioma.
In the present study, the accuracy signature was used in a large number of samples, which highlighted the reliability of the signature. Moreover, the predictive mechanism of the signature was determined from several perspectives, including the association between risk score with immune cell infiltration, TMB and IPS. However, the present study exhibits a number of limitations. An increased number of patients with glioma are required to participate in clinical trials to further validate the predictive capability of prognosis and the efficiency of immunotherapy of this signature.
In conclusion, an immune-related signature was constructed in the present study that can accurately predict prognosis and determine the response to immunotherapy in patients with glioma. The constructed signature may provide a deeper theoretical basis, and increase the accuracy of immunotherapy for patients with glioma.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.