Development of An Inammatory Response-Related Gene Signature To Predict Prognosis and Immunotherapy Response For Patients With Glioma

Inammatory response plays a crucial role in the development and progression of gliomas. However, the prognostic value of inammatory response-related genes has never been comprehensively investigated for glioma. In this study, we identied 39 differentially expressed genes (DEGs) between glioma and normal brain tissue samples, of which 31 inammatory response-related genes are related to the prognosis of glioma., The 8 optimal inammatory response-related genes were selected to construct prognostic inammatory response-related gene signature (IRGS) through the least absolute shrinkage and selection operator (LASSO) penalized Cox regression analysis. The effectiveness of the IRGS was veried in the training (TCGA) and validation (CGGA-693 CGGA-325 and Rembrandt) cohorts. The Kaplan-Meier curve revealed a signicant difference in the OS between the high- and low-risk groups. The receiver operating characteristic curve (ROC) shows the powerful predictive ability of IRGS. Meanwhile, a nomogram with better accuracy was established to predict overall survival (OS) based on the independent prognostic factors (IRGS, age, WHO grade, and 1p19q codeletion). In addition, patients in the high-risk group had higher immune, stroma, and ESTIMATE scores, lower tumor purity, higher inltration of immunosuppressive cells, higher expression of immune checkpoints, higher expression of TIDE and Exclusion, and lower expression of MSI Expe Sig. Thus, the patients in the low-risk group had signicantly higher respond rate of immune checkpoint inhibitors (ICIs). A novel prognostic signature incorporated 8 inammatory response-related genes was associated with the prognosis, immune landscape and the immunotherapy response in patients with gliomas. Thus, the signature can be suitable for future clinical application to predict the prognosis of patients with glioma.


Introduction
Glioma is the most common type of malignant tumor in the central nervous system, accounting for 40-50% of brain tumors, with an estimated annual incidence of 3 to 8 per 100,000 individuals [1,2]. Despite progress in comprehensive treatments, such as maximal tumor resection, adjuvant chemotherapy, adjuvant radiotherapy, and immunotherapy, the median survival time of glioma patients is still very short due to the aggressiveness of the glioma, high recurrence rate, and resistance to treatment [3,4]. Glioma is classi ed into 4 grades by the World Health Organization (WHO) according to their prognosis and morphological characteristics [5]. Especially glioblastoma (GBM), WHO grade IV, is the most malignant glioma, with the median overall survival (OS) of only 8 months, and the 5-year survival rate of 7.2% [1,6]. Traditionally, the WHO classi cation is considered an important criterion for the prognosis of gliomas.
However, the greatly differences for the clinical prognosis and treatment effectiveness of patients with the same WHO grade have been found due to the heterogeneity of gliomas [7]. Therefore, to solve this limitation, the construction of new biomarkers for accurately predicting treatment response and prognosis carries great clinical importance for patients with gliomas.
The appearance of white blood cells in tumors was initially discovered by Rudolf Virchow in the 19th century, which provides a connection between in ammatory response and cancers. With the rapid development of molecular biology, the important molecular mechanisms of in ammatory response in tumor formation have gradually been elucidated by some studies in recent year [8,9]. Wu et al. indicated that PTPN2 could be induced by in ammatory response and oxidative stress and its de ciency depressed glioma cell growth [10]. Meanwhile, in ammatory markers, derived from the routinely available parameters in the blood and their derivatives, appear to have the relationship of diagnostic and prognostic value in many neoplastic diseases including intracranial tumors [11][12][13]. For instance, white blood cells are widely existed in various tumors, especially malignant tumors, which are closely related to the important biological characteristics of tumors such as proliferation, migration, immune escape and prognosis [14].Likewise, several peripheral blood-derived biomarkers such as the neutrophil/lymphocyte ratio (NLR), platelet/lymphocyte ratio (PLR), monocyte-to-lymphocyte ratio (MLR) have been validated as prognostic in ammatory markers in various types of cancers [15][16][17]. Except for the in ammatory response factors in the blood, some in ammatory response-related genes are also used as biomarkers to predict the metastatic potential and prognosis of some cancers [18,19]. However, interestingly, the relationship between in ammatory response-related genes and the prognosis of patients with gliomas remains unclear.
Immunotherapy, emerged in recent years, is considered to be one of the most promising treatment methods. For example, a variety of therapeutic antibodies that block immune checkpoints, such as cytotoxic T lymphocyte-associated antigen 4 (CTLA4) and programmed cell death protein 1 (PD1), showed good effects in the treatment of many tumors such as liver cancer, non-small cell lung cancer, glioma and so on [20,21]. Unfortunately, immunotherapy based on immune checkpoints is not always useful for patients due to the heterogeneity of tumors and the complex tumor microenvironment. Therefore, it is greatly crucial and urgent to nd more immunotherapy targets and clarify the molecular mechanism of more detailed immunotherapy response.
With the tremendous progress of modern sequencing technology, many RNA-seq transcriptome data of cancers can be found in some public databases, such as The Cancer Genome Atlas (TCGA), Chinese Glioma Genome Atlas (CGGA), and Rembrandt databases. In this study, we systematically analyzed RNAseq transcriptome data and clinical data of glioma patients from these databases. Based on 8 in ammatory response-related genes, a prognostic in ammatory response-related gene signature (IRGS) of glioma patients was constructed in the TCGA database. Meanwhile, a reliable predictive nomogram model incorporating IRGS and clinicopathological characteristics was established and validated in the independent database. In addition, the correlations of IRGS with immune landscape, Tumor mutation burden, and the e cacy of immunotherapy were investigated. We consider that the IRGS has potential in patient management and can serve as potential therapeutic biomarkers for glioma in the clinical practice.

Patient data collection
The RNA-seq transcriptome data and corresponding clinical information of patients with glioma were extracted from TCGA, CGGA (CGGA-693 and CGGA-325), and Rembrandt databases. The RNA-seq Page 4/27 transcriptome data of TCGA and CGGA were quanti ed with Fragments Per Kilobase of transcript per Million mapped reads (FPKM) standardized by using the RSEM algorithm, whereas the Rembrandt data set was normalized microarray format [6]. Patients without survival data or OS less than 30 days were excluded from analysis for which they might die of the acute complications, rather than of the glioma itself. Therefore, a total of 1896 glioma patients were included in this study. The TCGA data set (n=631) was used as the training cohort, the CGGA-693 data set (n=656), CGGA-325 data set (n=309) and Rembrandt data set (n=300) were used as the veri cation cohorts. The clinicopathological characteristics of the included glioma patients are shown in Table 1. Then, the RNA-seq transcriptome data of normal brain tissue samples were searched from the Genotype-Tissue Expression (GTEx) database. In addition, the 200 in ammatory response-related genes were collected from the Molecular Signatures Database V7.4 (https://www.gsea-msigdb.org/gsea/msigdb). Further details on these genes are summarized in the Supplementary Table S1. If the gene expression data is detected in less than 50% of the samples, the gene is excluded in order to facilitate the subsequent data analysis.

Identi cation of differentially expressed genes and cluster analysis
The differentially expressed genes (DEGs) between glioma samples in TCGA cohort and healthy brain tissue samples in GTEx cohort were identi ed through the differential-expression analysis with R package "limma" (log2 fold-change |logFC|>1 and an adjusted false-discovery rate FDR <0.05). Furthermore, Search Tool for the Retrieval of Interacting Genes (STRING) software (http://www.string-db.org/) was used to analyze the interactions and evaluate the level of interactions among the above DEGs.
In order to investigate the best gene subgroups, the glioma patients of TCGA was divided into different clusters through the application of the "ConensusClusterPlus" package in R based on the expression of DEGs. The resampling procedure was used to sample 80% of the sample 50 times, the similarity distance between the samples was estimated by Euclidean distance, and the "km dist" as a clustering algorithm was used to select reliable and stable subgroup classi cation [22]. The "proportion of ambiguous clustering" was used to select an optimal value of clusters with the criteria is that the clusters have high consistency, and the area under the cumulative distribution function curve does not increase signi cantly [23].

Development and validation of a prognostic IRGS
The univariate Cox regression analysis was conducted by using the survival coxph function of the R package to explore the relationship between DEGs and prognosis value (p<0.05 was used as the signi cance threshold). Then the least absolute shrinkage and selection operator (LASSO) penalized Cox regression analysis was performed to construct a prognostic model in order to minimize the risk of over tting [24,25], which was performed within the TCGA cohort by using the R package "glmnet" and "survival". The ten-fold cross-validation was used to selecting the optimal penalty parameter λ of the prognostic model and was followed the minimum criteria (i.e. the value of λ corresponding to the lowest partial likelihood deviance). Subsequently, the risk score of each patient was calculated according to the regression coe cient of the model and the expression of the corresponding gene. The calculation formula is shown below: where represents the number of all the selected gene; represents the serial number of each gene; and refer to the expression level of each selected gene and corresponding coe cient, respectively.
The median risk score was considered as the cut-off value to divide the patients into high-or low-risk group. The Kaplan-Meier survival curve analysis with log-rank test were used to evaluate the ability of predicting OS between the high-and low-risk groups by using the R packages " survival " and "Survminer". The receiver operating characteristic (ROC) curve and the area under the ROC curve (AUC) were calculated to evaluate the accuracy of the in ammatory response-related gene through the R software package "timeROC". The above analyses were performed simultaneously in the training and validation cohorts.

Establishment and pvaluation of a nomogram
The univariate and multivariate Cox regression analyses were performed for identifying the independent prognostic factors in terms of the risk scores and clinicopathological characteristics. Next, a nomogram was established on the independent prognostic factors in the TCGA cohort by using the R packages "rms", "regplot", and "Hmisc". The calibration curves and C-index were performed in the training and validation cohorts to evaluate the availability of this nomogram [26]. In addition, The ROC curves were also ploted to assess the accuracy of the nomogram for OS prediction in the training and validation cohorts.

Functional enrichment analysis
The DEGs between the high-risk group and the low-risk group were identi ed with the criteria of |log2FC|>2 and adjusted p<0.05 by BH method. Then Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted to predict the function of these DEGs through the R packages "clusterPro ler", "org.Hs.eg.db", and "rich lot".

Evaluation of immune landscape of glioma microenvironment
The in ltration levels of immune cells and stromal cells in different tumor tissues were analyzed by immune score and stromal score. The immune scores and stromal scores of glioma patients were calculated using the ESTIMATE algorithm via the R package "estimate" [27]. CIBERSORT, a deconvolution algorithm with 1000 permutations, was applied to estimate the abundance of 22

Evaluation of tumor mutation burden and immune response
Tumor mutation burden (TMB) has been de ned as the total number of somatic, coding, base substitution, and indel mutations per megabase of genome examined. In this study, we selected a novel TME-scoring algorithm to explore the differences in TME between the high-and low-risk groups. Furthermore, the tumor immune dysfunction and exclusion (TIDE) algorithm was applied to evaluate the predictive e ciency of risk scores for the Immune checkpoint inhibitions (ICIs) response in glioma [29]. In addition, the distributions of risk scores between non-respond and respond groups, and the differences of response rates between high-and low-risk groups are also analyzed and compared.

Statistical analysis
All statistical analyses and graph visualization were performed using the R software (version 4.0.1, http://www.R-project.org). The PERL programming language (version, 5.30.2, http://www.perl.org) was used to process RNA-seq transcriptome data. The Chi-square test or two-sided Fisher's exact test was performed for the comparison of categorical variables between the high-and low-risk groups, similarly, correlations between the risk score and clinicopathological parameters were the same. The Student's ttest or one-way ANOVA test was utilized to compare the continuous variables with normal distribution between two groups or more than two groups. While Mann-Whitney U test was used to compare continuous variables with non-normal distribution between two groups, and Kruskal Wallis test was used to compare continuous variables with non-normal distribution in more than two groups. The survival advantage was evaluated using the Kaplan-Meier curve and compared by the log-rank test for the highand low-risk groups. Univariate and multivariate cox proportional hazard models were applied to estimate the hazard ratios of variables and determine independent prognostic factors. The statistical signi cance level was set at two-tailed p<0.05.

Identi cation of the DEGs between glioma and healthy brain tissue samples
These data were pre-processed before identifying the DEGs, the normalization and batch effects removal from GTEx and TCGA cohorts was conducted by R package "sva". 39 DEGs were identi ed between in glioma samples (n = 697) and healthy brain tissue samples (n = 1157), with the absolute value of the log2-transformed fold change (FC) > 1 and the adjusted P-value (FDR) < 0.05 as the threshold levels of signi cance. Among of the identi ed in ammatory response-related DEGs, 17 up-regulated genes and 22 down-regulated genes were identi ed, which indicated that in ammatory response-related genes exerted important biological functions in the tumorigenesis of glioma. The heatmap and volcano plot of the DEGs are shown in Figure 1A and B. In addition, the network and correlation analysis of 39 in ammatory response-related genes were performed to identify their interactions ( Figure 1C and D).

The clinicopathological characteristics of the two clusters of glioma patients.
According to the expression similarity of in ammatory response-related genes and the criteria for selecting the number of clusters, k=2 was considered as the most appropriate value, with clustering stability increasing from k = 2 to 10 in the TCGA datasets (Figures 2A-C). Hence the patients from TCGA cohort were divided into two clusters named cluster 1 and 2. Subsequently, the principal component analysis (PCA) was used to compare the transcriptional pro le between cluster 1 and cluster 2 subgroups. The results displayed that there is a clear boundary between the two subgroups ( Figure 2D).
In order to compare the survival distributions of the two clusters, the Kaplan-Meier survival curves of the them performed with log-test showed the signi cantly shorter OS in the cluster 1 subgroup than the cluster 2 subgroup ( Figure 2E, p<0.001). Moreover, the clinicopathological characteristics, survival status and immune scores between different clusters were identi ed, and 39 in ammatory response-related gene expression values were screened out as a heat map ( Figure 2F). The cluster 1 subgroup is signi cantly correlated with dead (p < 0.001), older patients (p < 0.001), higher grade (p < 0.001), glioblastoma and astrocytoma (p < 0.0001), IDH-Wild type (p < 0.001), Non-codel 1p19q codeletion (p < 0.001), Unmethylated MGMT status promoter (p < 0.001), higher Immune score (p < 0.001), higher Stromal score (p < 0.001), higher ESTIMATES score (p < 0.001), lower Tumor purity (p < 0.001). While cluster 2 subgroup is signi cantly correlated with the opposite results.

Construction and validation of the prognostic IRGS
Univariate Cox regression analyses were rst performed for in ammatory response-related DEGs, and the results showed that 31 genes were signi cantly related to the OS of glioma in the Supplementary Table S2 (p≤0.05). Among these 31 genes, 22 were protective factors and 9 were risk factors for prognosis. Above-mentioned 31 in ammatory response-related genes were incorporated into the LASSO penalized Cox regression analysis in the TCGA cohort. Ultimately, 8 optimal prognostic in ammatory responserelated genes were acquired and included in the risk model, including GABBR1, CALCRL, EBI3, BTG2, SEMA4D, SELL, SERPINE1, and MMP14 ( Figures 3A-B). The survival analyses of these 8 in ammatory response-related genes were indicated that high expression of these genes was all signi cantly correlated with poor OS in Supplementary Figure S1.  Figure 3C). Subsequently, the median risk score, as the cut-off value, was used to strati ed the glioma patients into the high-and low-risk groups. In the TCGA cohort, the Kaplan-Meier curve demonstrated that the survival in the low-risk group was signi cantly better than in the high-risk group (log-rank test p < 0.001; Figure 3D). The distribution plot of the risk score and survival status showed that the higher the risk score, the more deaths of glioma patients ( Figure 3H). A satisfactory prediction performance of the prognostic model was con rmed by the AUC for 1-, 3-, and 5year OS (AUC = 0.876, 0.903, and 0.842, respectively; Figure 3L).
To test the stability of the prognostic model constructed on the TCGA cohort, the same analyses were performed in the CGGA-693 cohort, CGGA-325 cohort and the Rembrandt cohort. Similarly, patients in the low-risk group had better survival outcomes than the patients in the high-risk group (Figures 3E-G, I-K).

Correlation analysis between the prognostic model and clinicopathological characteristics
The association of the 8 selected genes, and clinicopathologic parameters between high-and low-risk groups was analyzed in the TCGA cohort, and the results showed that the genes were up-regulated in the high-risk group and the protective genes were up-regulated in the low-risk group ( Figure 4A). Meanwhile, the differences of age, gender, histology, WHO grade, IDH status, 1p19q codeletion and MGMT promoter status between high-and low-risk groups were compared in all training cohort and validation cohorts and showed in the Supplementary Tables S3-6. In addition, the values of risk score between glioma patients strati ed by clusters and various clinicopathological characteristics were also compared. In the TCGA cohort, glioma patients with the clinicopathological characteristics of the cluster 1 subgroup, age >45 years, more malignant type of histology, higher WHO grade, IDH wild type, 1p19q codeletion non-codel and MGMT promoter unmethylated showed signi cantly higher levels of risk score, while no differences were observed between patients satis ed by gender ( Figures 4B-I). Likewise, the risk score of glioma patients in the CGGA-693, CGGA-325 and Rembrandt cohorts was identi ed the similar results to the TCGA cohort ( Supplementary Figures S2A-C).
We performed subgroup survival analyses to determine whether the clinicopathological characteristics would decrease the prediction accuracy of the prognostic IRGS. The results showed that patients with high-risk score had worse survival outcomes than those with low-risk score in all subgroups, except for the grade IV subgroup in the TCGA cohort ( Supplementary Figures S3A-M). The other cohorts showed the same results but IDH-Wild subgroup in the CGGA-693 cohort and the 1p19q-Noncodel subgroup in the Rembrandt cohort ( Supplementary Figures S4A-C).

Establishment and evaluation of a nomogram
The univariate Cox regression and multivariate Cox regression were analysed subsequently for identifying the OS-related factors in the TCGA, CGGA-693, CGGA-325, and Rembrandt cohorts (Figures 5A-D). And the risk score based on 8 selected genes was con rmed as an independent prognostic factor without being affected by other clinicopathological characteristics (Figures 5A-D). Then, the independent prognostic factors in the TCGA cohort (Age, WHO grade, 1p19q and Risk score) were used to establish the nomogram for predicting the 1-, 3-, and 5-year prognoses of the glioma patients ( Figure 6A). The internal evaluation was initially performed. The C-index was 0.857 (95%CI: 0.812-0.896) and the calibration curves indicated an excellent match between the actual and nomogram-predicted probability of 1-, 3-, and 5-year OS ( Figure 6B). The ROC curve analysis showed that the nomogram had highest predictive accuracy of 5-year OS than risk score or clinicopathological features ( Figure 6F).  Figures 6C-E). Furthermore, the ROC curve analysis showed that the nomogram had highest predictive accuracy of 5-year OS than risk score or clinicopathological features in three validation cohorts ( Figures 6G-I). Meanwhile, The ROC curve analysis also exhibited that the nomograms had the highest accuracy in predicting 1-, and 3--year OS in comparison to other independent prognostic factors (Supplementary Figure S5). Thus, the nomogram has the potential as a quantitative method to predict the prognosis of patients with glioma.

Functional enrichment analysis
The GO function and KEGG pathway enrichment analyses were conducted to characterize the biological functions of DEGs between the high-and low-risk groups. As a result, the GO function analyses revealed signi cant enrichment of in ammatory response-related functions, including response to molecule of bacterial origin, response to lipopolysaccharide, and positive regulation of cytokine production ( Figure  7A). Meanwhile, the DEGs were also signi cantly enriched in several immune-related biological processes, for instances T cell activation and regulation of immune effector process ( Figure 7A). Besides, the KEGG pathway analyses identi ed that some KEGG pathways were correlative with immune-related pathways including the Cytokine-cytokine receptor interaction, Chemokine signaling pathway and Toll-like receptor signaling pathway ( Figure 7B).

Correlation of the prognostic IGRS with the immune landscape of glioma microenvironment and immunotherapy
The correlation of the prognostic IGRS with the immune landscape of glioma microenvironment was further investigated for the immune-related functions enriched by DEGs between high-and low-risk groups. Firstly, the immune system scores were compared between high-and low-risk groups, and the result indicated that the high-risk group showed signi cantly higher immune, stroma and ESTIMATE scores and lower tumor purity than the low-risk group (Figures 7C-F). Then, we analyzed the different extent of immune cell in ltrations between the low-and high-risk groups for all the samples. The differences were observed in the high-risk group with lower abundance of naive B cells, activated NK cells, activated mast cells, and eosinophils, but higher abundance of memory B cells, plasma cells, CD8+ T cells, CD4+ memory activated T cells, regulatory T cells, gamma delta T cells, M0-type macrophages, M1type macrophages, resting mast cells and neutrophils ( Figure 7G).
The expressions of immune checkpoint proteins were all up-regulated in the high-risk group ( Figure 8A). We next determined whether a correlation existed between the immune checkpoint protein and prognostic in ammatory response-related genes. The correlation matrix of and was demonstrated that the immune checkpoint protein had a signi cant negative correlation with the in ammatory response-related genes such as GABBR1, CALCRL, EBI3, BTG2, SEMA4D, and SELL, except for SERPINE1 and MMP14 ( Figure  8B). Moreover, a comparative analysis of tumor mutation burden in the two subgroups showed the highrisk group had signi cantly higher tumor mutation burden ( Figure 8C).
Recently, more and more study has reported the ICIs targeting immune-checkpoints such as PD-1 and PD-L1 could improve the e ciency in treatment of tumors. We found that the expression of TIDE and Exclusion was signi cantly higher in the high-risk group, while the expression of MSI Expe Sig was signi cantly lower, with the comparison to the low-risk group ( Figure 8D). However, there was no statistically difference about the expression of Dysfunction between the two groups ( Figure 8D). Then, the distributions of risk scores between non-respond and respond groups was compared and indicated the non-respond group with signi cantly higher risk scores, which inferred that patients with no respond to hemotherapy had the poor prognosis ( Figure 8E). The TIDE algorithm, which was conducted to predict the immunotherapy responders through transcriptomic data, was used to explore whether prognostic genes could predict immunotherapeutic bene t in glioma. The ROC curve showed the superior predictive performance with AUC value of 0.833 in predict immunotherapeutic bene t ( Figure 8F). Meanwhile, the result revealed that the response rates to ICI treatment were signi cantly higher in low-risk patients (92%) compared with high-risk patients (44%) (p<0.001, Figure 8G). Thus, the more survival bene t was likely to obtained from immunotherapy for patients in the low-risk group.

Discussion
Glioma, the most common malignant brain tumor, is notorious for its characteristics of low cure rate and high recurrence rate [30]. Currently, the problem faced by most neurosurgeons is that the prognosis of most patients with glioma has not been signi cantly improved by through the existing universal treatment plan. The high heterogeneity of glioma is the main limitation, which leads to inconsistent treatment response and prognosis. Therefore, it is very important and urgent to develop personalized treatment plans for glioma patients. The rapid development of bioinformatics and sequencing technology provides a new perspective to solve this problem. Previous studies indicated that hypoxiarelated gene signature [31], m6A-related gene signature [32], ferroptosis-related gene signature [33], immune-related gene signature [34], and energy metabolism-related gene signature [35] are used as prognostic biomarkers to predict the prognosis of glioma. No matter prognosis of 1-year OS, 3-year OS or 5-year OS, they all exhibited superior predictive performance, which makes researchers con dent in individually predicting the prognosis of patient by genomics. In ammatory response, just like the above biological processes, is also an important factor in the occurrence, development and treatment of various tumor cells[36]. Lin et al. [19] constructed prognostic signature based on 8 in ammatory response genes to predict the prognosis of liver cancer, and the signature showed the OS of 1,3,5-year with AUC is 0.685, 0.626, 0605, respectively. Nevertheless, the in ammatory response-related gene signature as prognostic marker for glioma is still unclear.
In this study, 39 DEGs between glioma and normal brain tissues were differentiated from the 200 in ammatory response-related genes. Then 31 prognostic DEGs were identi ed after univariate Cox analysis, 8 of them were selected to construct prognostic IRGS. Whether in the training cohort (TCGA) or the validation cohorts (CGGA-693, CGGA-325 and Rembrandt), IRGS has been con rmed the strong capability in predicting survival outcomes of glioma patients. Subsequently, the prognostic IRGS was incorporated with other independent prognostic factors (age, WHO grade, and 1p19q codeletion) to construct a nomogram with superior OS prediction ability. The functional enrichment analysis revealed the potential difference in in ammatory response sensitivity between the high-risk group and the low-risk group. Meanwhile, immune-related biological processes and pathways were also observed through functional enrichment analysis. Thus, we further revealed the differential immune landscape between the two risk subgroups by comparing the immune, stromal and ESTIMATES scores, tumor purity, abundance of immune cells, and expression levels of immunoregulatory molecules. In addition, the assessment of tumor mutation burden, immune checkpoints, and immunotherapy response are also analyzed in detail based on IRGS risk strati cation. To our knowledge, this is the rst study to predict the prognosis of glioma with IRGS.
The prognostic IRGS constructed in our study incorporated 8 in ammatory response-related gene such as GABBR1, CALCRL, EBI3, BTG2, SEMA4D, SELL, SERPINE1, and MMP14. Except for GABBR1 and SEMA4D, other genes are up-regulated in glioma tissues. Among these genes, SERPINE1 and MMP14 are associated with poor prognosis, whereas the remaining 6 genes with good prognosis. GABBR1 is a metabotropic G protein-coupled receptor that mediates the inhibitory effect of γ-aminobutyric acid, Yang et al. demonstrated that 5 miRNAs promote the proliferation and invasion of colorectal cancer cells by inhibiting GABBR1 [37]. As a G protein-coupled seven transmembrane domain receptor, CALCRL has been identi ed as a potential tumor suppressor for lung adenocarcinoma by Lu et al.[38]. A β-subunit of IL-12 family member (EBI3) was con rmed to promote T-and B-cell differentiation through the gp130-STAT3 signal pathway, and thereby playing a crucial role in inhibiting tumor cell proliferation and invasion [39]. In the non-small-cell lung cancer models, overexpression of BTG2 prevented the lung cancer cell metastasis by inhibiting cell invasion [40]. SEMA4D is an important member of the semaphorin subfamily and plays an important role in the nervous and immune systems. Chen et al. [41] rstly reported that highly expressed SEMA4D plays an important role in osteolytic bone metastasis of lung cancer by inhibiting osteoblast differentiation, thereby providing a potential strategy for targeting osteoblasts to treat bone metastasis. However, some studies declared to promote tumor cell invasion and metastasis by tumor angiogenesis in some cancers such as prostate cancer and colon cancer [42]. SELL, encoding the selectin L protein (CD62L), mediates the adhesion of cells to the vascular endothelium [43]. While the role of SELL in tumorigenesis and development has not been reported yet. Conversely, SERPINE1, as a member of the superfamily encoding serine protease inhibitors, works by inhibiting proteolytic activity and promoting angiogenesis. Previous studies demonstrated that the overexpression of SERPINE1 may lead to the spread and metastasis of colon cancer, and is a poor prognostic indicator for malignant tumors such as breast cancer and gastric adenocarcinoma [44,45]. MMP14 is also named as membrane type 1 matrix metalloproteinase, which participates in the degradation of basement membrane and extracellular matrix, and thus promotes tumor invasion and metastasis [46]. It was of great signi cance that these genes closely correlated to the origin and metastasis of some cancers, while the detailed action mechanism remains to further investigate.
These in ammatory response genes are involved in the biological processes of many cancers and have predominant predictive performance as a prognostic signature. However, whether these genes affect the prognosis of glioma patients through in ammatory response remains to be clari ed. The KEGG analysis revealed that DEGs between different risk subgroups were signi cantly enriched in pathways related to in ammatory response and cancer, which further con rmed that the in ammatory response has a closely connection with tumor progression. Furthermore, the DEGs between different risk subgroups were also enriched in many immune-related biological processes such as T cell activation and regulation of immune effector process and signal pathways, which remind us that in ammatory response-related genes to may be related to immunity. Thus, we subsequently analyzed the immune scores and immune cell in ltration between the two risk subgroups. Further analyses found that the high-risk group exhibited higher immune scores, higher abundance of immunosuppressive cells (Treg, CD8+ cell) [47], and higher expression levels of immune checkpoints and macrophage associated molecules. In contrast, tumor killer cells (activated NK cells) showed a higher abundance in the low-risk group. To some extent, these results illustrated that in ammatory response-related genes is related to the immune landscape of the glioma microenvironment. From the above results, it can be concluded that the antitumor immunity of high-risk group is signi cantly weakened, and thus we speculated that this may be an important reason for their poor prognosis. However, the potential molecular mechanism between in ammatory response-related genes and glioma immunity has not been elucidated, and further research is needed to investigate.
With the continuous improvement of gene sequencing technology in the last few decades, the knowledge of researchers for molecular structure has continued to strengthen. Cancer treatment has shifted from chemotherapy and radiotherapy targeting tumors broadly to antibody-based immunotherapies that modulate immune responses against tumors more precisely [48]. ICIs, as the rst generation of immunotherapy, play a pivotal role by blocking receptor and/or ligand interactions of molecules, such as PD-1/ PD-L1 pathway and CTLA-4 [49]. Some prior studies have been described the therapeutic antibodies targeted therapy for immune checkpoints, which indicated robust and durable responses in patients with many cancers [50], including glioma [51]. In the current study, the high-risk group had a signi cantly higher immune checkpoints than the low-risk group, and the immune checkpoints was positively correlated with the expression of genes related to poor prognosis, and negatively correlated with genes related to good prognosis. In addition, the response rate to ICIs in the high-risk group was signi cantly lower than the lowrisk group. Therefore, our prognostic model showed superior capacity in predicting the expression level of immune checkpoints, and had the potential to guide immunotherapy.
Certainly, there are some limitations warranted mention in this study. Firstly, it is a retrospective study, thus the results should be further validated its clinical utility in prospective studies. Secondly, these four cohorts have various degrees of de ciencies in clinical information, especially, the Rembrandt cohort is more serious, leading to insu cient veri cation of some research results. In addition, we were not able to validate in clinical practice owing to the lack of glioma samples. Thirdly, GO and KEGG enrichment analysis and subsequent analysis of immune microenvironment and immune checkpoints are only evaluated in the TCGA cohort, not veri ed in other cohorts. Lastly, the potential molecular mechanism between in ammatory response-related gene and immunotherapy has not yet been elucidated, and further experiments are needed to explore it.

Conclusion
In conclusion, we performed a comprehensive analysis of the in ammatory response-related gene in the four independent cohorts and established a new prognostic signature incorporating 8 in ammationrelated genes for glioma patients. The signature exhibited robust capacity in predicting survival outcomes of glioma patients, and was correlated with functional analysis, immune landscape of glioma microenvironment, and immunotherapy. Accordingly, our study provides insight in in ammatory response-related gene for predicting the prognosis of glioma, which is crucial for the development of personalized glioma treatment therapies.

Declarations
Ethics approval and consent to participate: The study received full approval from the ethics committee of