Development and Validation of an Immune-Related Gene Pairs Signature in Grade II/III Glioma

Introduction

Gliomas are prevalent primary intracerebral malignant tumors.¹ To date, surgical resection combined with postoperative chemoradiotherapy is the standard treatment for Grade II/III gliomas. Although much effort has been put forth to improve the clinical outcome, more than one-half of people diagnosed with Grade II/III glioma experience recurrence or progression to high-grade glioma over time.² Moreover, in addition to conventional surgical therapy, patients with Grade II/III glioma have limited sensitivity to chemotherapy and radiotherapy.^3,4 Hence, it is suggested that other prognostic factors for Grade II/III glioma remain to be discovered. Several biomarkers, including codeletion of chromosome arms 1p and 19q (1p/19q codeletion), O-6-methylguanine-DNA methyltransferase (MGMT) methylation and isocitrate dehydrogenase (IDH) mutation, have been added to the 2016 WHO classification to clarify the histological features and to guide the therapeutic plan for Grade II/III glioma patients.^5–8 However, these widely recognized biomarkers cannot be effectively used to evaluate individual risk stratification for people with Grade II/III glioma.

The tumor microenvironment (TME) has attracted considerable attention as a vital cellular milieu due to its pivotal roles in tumor initiation, progression, and metastasis. The TME consists of mesenchymal cells, inflammatory mediators, endothelial cells, immune cells and stromal cells.⁹ Additionally, immune cells and stromal cells are two main nontumor components that have great importance in the diagnosis and prognosis of tumors.¹⁰ A prevalent algorithm called ESTIMATE designed by Yoshihara et al has been used to evaluate the types of immune cells and stromal cells residing in many malignant tumors.^11–13 In this study, we also employed this algorithm for our research.

The glioma immune microenvironment plays an essential role in glioma biology.¹⁴ Immune checkpoint inhibitors and immunotherapy play roles in improving the glioma prognosis in patients.^15,16 Regarding prognostic biomarkers, researchers have explored the possibility of stratifying patients with malignant tumors (including gliomas) based on gene expression profiles and have constructed a multigene expression model that can be applied to stratify high-risk subgroups.^17–22 One common drawback of these studies is that gene expression profiles require suitable normalization. However, owing to technical biases across sequencing platforms and biological heterogeneity, the profiles obtained inevitably show variations.²³ Thus, researchers have developed the concept of immune-related gene pairs, which can eliminate the disadvantages of data processing and is considered to be a robust prognostic indicator for some tumors.^24–29

However, the prognostic value of these immune-related gene pairs has not been estimated in patients with Grade II/III glioma. Our study constructed and validated a personal prognostic indicator by integrating immune-related gene pairs found in Grade II/III glioma.

Materials and Methods

Acquisition of Grade II/III Glioma Expression Profiles from TCGA and CGGA Datasets

The RNA-seq fragments per kilobase per million (FPKM) expression data and clinical information of 515 Grade II/III glioma patient samples were downloaded and selected from the TCGA database (https://portal.gdc.cancer.gov). Then we deleted 104 non-frozen specimens and used the caret package in R software to randomly divide the TCGA data into two cohorts. The mRNA sequencing and clinical data of 605 Grade II/III glioma samples were downloaded from the CGGA database (http://www.cgga.org.cn).

Acquisition of Immune-Related Genes

We downloaded a comprehensive list of immune-related genes (IRGs) from the Immunology Database and Analysis Portal (ImmPort) database (https://www.immport.org). Then, the limma package in R software was utilized to extract the RNA expression profiles of these IRGs (determined by median absolute deviation, MAD > 0.5).³⁰

Establishment of a Prognostic Signature Based on Immune-Related Gene Pairs

A pairwise comparison of each sample’s immune-related gene expression levels was performed to acquire a score for each IRGP. IRGP scores were determined as follows: (1) if IRG 1 > IRG 2, then IRGP = 1 and (2) if IRG1 < IRG 2, then IRGP = 0. The superiority of analyzing genes in a pairwise method is that the need for normalization steps for individualized analysis is eliminated.²⁹ If the score of an IRGP is 0 or 1 in more than 80% of the samples obtained from the TCGA or CGGA database, then we discarded the IRGP. Then, we applied least absolute shrinkage and selection operator (LASSO) Cox regression (iteration = 1000) to generate an IRGP index (IRGPI) by using the R package glmnet. To stratify patients into high- and low-risk groups, the optimal cutoff value for the IRGPI was determined by receiver operating characteristic (ROC) curve analysis of 3-year overall survival data in the TCGA training cohort.

Validation of the IRGPI

We further evaluated the IRGPI model with the validation cohorts consisting of Grade II/III glioma samples in the TCGA database and CGGA database by Log rank test. Additionally, we assessed the IRGPI and other clinical variables through univariate and multivariate Cox proportional hazard regression models using TCGA and CGGA datasets.

Profiling of Infiltrating Immune Cells

To clarify immune cell infiltration in different risk groups, we used CIBERSORT,³¹ an algorithm frequently used for exploring tumor tissue cell composition on the basis of gene expression profiles. The gene expression profile was uploaded to the online analytical platform CIBERSORT web portal (http://cibersort.stanford.edu) using the default signature matrix at 1000 permutations. Then, p-values were calculated by using the “wilcox.test” function in R.

Gene Set Enrichment Analysis and Filtration of the Results

To understand the biological functions of the immune-related prognostic signature, we performed gene set enrichment analysis (GSEA/MSigDB, http://software.broadinstitute.org/gsea/downloads.jsp). A false discovery rate-adjusted P<0.05 was considered statistically significant. To obtain a comprehensive understanding of some biological functions, we filtered the outcome of GSEA from the perspective of the tumor microenvironment. Immune and stromal scores were calculated using the ESTIMATE algorithm,¹¹ and immune scores and stromal scores were stratified into low-level groups and high-level groups according to their median value. Then, we acquired the differentially expressed genes (DEGs) common to the stromal scores and immune scores, as determined by the VennDiagram package.³² Additionally, we conducted gene ontology (GO) enrichment analysis based on these DEGs. Finally, we evaluated the results between the GO analysis and GSEA by using the Venn Diagram package.

Statistical Analysis

All statistical analyses were performed using R software (version 3.5.1, https://www.r-project.org). For all tests, a p-value < 0.05 indicated a significant difference. Significance is shown as *P<0.05, **P<0.01, and ***P<0.001.

Results

Prognostic IRGP Signature Construction

Data on 1119 patients with Grade II/III glioma were included in the analysis (Supplementary Table 1). A TCGA cohort was used as the training dataset. Altogether, 1991 immune-related genes (IRGs) were downloaded from the ImmPort database (accessed on 10/5/2020). A total of 372 IRGs were identified in both the TCGA and CGGA datasets (MAD > 0.5). From these 372 IRGs, 15,245 IRGPs were constructed. After discarding IRGPs with relatively small variation (MAD = 0), 22 IRGPs remained, and they were selected as the initial IRGPs. Then, we defined the IRGPI by Lasso Cox proportional hazards regression using the training set and selected 8 IRGPs for the final model. The IRGPI consisted of 15 unique IRGs (Table 1). The optimal cutoff value in the IRGPI for the prediction of Grade II/III glioma prognosis was determined by ROC curve data analysis to be 0.395 (Supplementary Figure 1). ROC curves were also analyzed to validate the impact of clinical characteristics on the risk score in the training cohort of the TCGA dataset (Figure 1A) and in the CGGA dataset (Figure 1B). The five-year ROC curve in the training cohort has an AUC of 0.864 (Figure 1C). The AUC of five-year ROC curve in the external validation cohort is 0.701 (Figure 1D). The IRGPI significantly stratified patients into low- and high-risk groups on the basis of overall survival (OS). We found that in the TCGA training cohort, the high-risk group had a profoundly inferior OS than the low-risk group (P<0.001, Figure 2A). To further investigate the prognostic value of the IRGPI, we applied univariate and multivariate Cox regression to the training cohort. After combining with clinical factors such as age, sex, and tumor grade, the risk score of the IRGPI remained an independent prognostic factor (P<0.001, HR, 2.684 [2.230, 4.046]; Figure 2D and G).

Table 1 Model Information About IRGPI

Figure 1 A comparison of 5-year ROC curves with other common clinical characteristics showed the superiority of the risk score in the training cohort (A) and in the external validation cohort (B). The 1-, 3-, and 5-year survival ROC curves of training (C) and external validation cohorts (D).

Figure 2 Analysis of survival and independent prognostic factors. Kaplan-Meier curves of overall survival (OS) among different IRGPI risk groups (low vs high risk). OS among patients in training (A), internal (B) and external validation cohorts (C). The risk score of IRGPI was an independent prognostic factor among univariate (D) and multivariate (G) analyses in the training database. The risk score of IRGPI was an independent prognostic factor among univariate (E) and multivariate (H) analyses in the internal validation cohort. The risk score of IRGPI was also an independent prognostic factor among univariate (F) and multivariate (I) analyses in the external validation cohort.

Validation of the IRGPI for Use in Survival Prediction

To confirm whether the IRGPI shows consistent prognostic value in different populations, we applied the IRGPI signature to the TCGA validation cohort for internal validation and to the CGGA validation cohort for external validation. Unsurprisingly, in both validation cohorts of TCGA data and the CGGA data, the high-risk patients exhibited a significantly shorter OS than the low-risk patients (P<0.001, Figure 2B and C). Additionally, the IRGPI remained an independent prognostic factor after taking clinical factors such as age, sex, tumor grade, IDH and 1p/19q into consideration (P<0.05, HR, 2.301 [1.164, 4.547]; Figure 2E and H), (P<0.05, HR, 1.369 [1.047, 1.791]; Figure 2F and I). The association between risk score, IDH mutation status and 1p/19q codeletion status was further evaluated, and the results showed that the groups with IDH mutations and 1p/19q codeletion had lower risk scores (P<0.001, Figure 3A and B).

Figure 3 The association between risk score with IDH mutation status (A) and 1p/19q co-deletion status (B).

Immune Cell Infiltration Between Different Risk Groups

We carried out an immune cell infiltration analysis to better understand the difference in immune cell infiltration between the low-risk and high-risk groups. Figure 4A presents a summary of the comparison of the CIBERSORT results for the two risk groups. We found that M1 macrophages and CD8 T cells were significantly more highly expressed in the high-risk group (P=1.757e-04, Figure 4B; P=0.002, Figure 4D), while monocytes were significantly more highly expressed in the low-risk group (P=5.749e-05, Figure 4C). Moreover, Macrophage M0 and T cells CD4 memory activated were also highly expressed in the high-risk group in our study (P=0.004, Supplementary Figure 2A; P=0.033, Supplementary Figure 2B).

Figure 4 Immune infiltration situation between IRGPI risk groups. (A) Summary of the outcome estimated by CIBERSORT in different risk groups. (B) Macrophages M1 was significantly highly expressed in the high-risk group (P=1.757e-04). (C) Monocytes were significantly highly expressed in the low-risk group (P=5.749e-05). (D) T cells CD8 was significantly highly expressed in the high-risk group (P=0.002).

Gene Set Enrichment Analysis and Filtration of the Results

We carried out a GSEA between the low- and high-risk groups to explore the pathways that were significantly different. We found that 299 pathways were quite different (Supplementary Table 2). Then, we selected 50 pathways with minimally adjusted P-values for filtration. According to the statistical analysis, 1199 upregulated DEGs and 296 downregulated DEGs in the stromal score group were selected. A total of 965 upregulated DEGs and 692 downregulated DEGs were selected in the immune score group. Then, 887 upregulated intersecting genes and 277 downregulated intersecting genes were identified in Venn plots (Figure 5A and B). The Gene Ontology (GO) analysis of these intersecting genes is shown in Figure 5C. Finally, we examined the intersection of GO and GSEA results and identified 26 functional pathways (Supplementary Figure 3, Table 2).

Table 2 The Intersection Functional Pathways Between GO and GSEA

Figure 5 Differentially expressed genes analyses in TCGA dataset. (A) 887 up-regulated intersected genes and (B) 277 down-regulated intersected genes were revealed in venn plots. The gene ontology (GO) analysis of these intersected genes (C).

Discussion

We constructed an immune-related gene pair signature in the present study and validated this signature in internal and external datasets. Our research showed that the signature can be used to stratify patients into low- and high-risk groups. Univariate and multivariate Cox regression analyses showed that the risk score was an independent prognostic factor. This prognostic signature related to the tumor immune microenvironment offers a new perspective on developing novel predictive biomarkers and improving Grade II/III glioma patient management in the era of immunotherapy.

Similar to other studies, several immune cells were distributed significantly differently in our research. For example, we found that CD8 T cells and M0 macrophages were greatly enriched in the high-risk group. Wu et al²⁷ also showed the same trend, although it was not statistically significant. In contrast, results from one study indicated that CD8 T cells were significantly higher in the low-risk group.²⁸ This disparity may be a result of differences in immune-related gene pairs and tumor types in the aforementioned studies. Thus, further exploration is needed. Additionally, we identified 26 functional pathways from the intersection of the GO analysis and GSEA results, which may indicate that these immune-related gene pairs function through these pathways.

Despite recently improved prognosis for certain Grade II/III glioma patients, the disease in 70% of patients inevitably progresses to high‐grade aggressive glioma and eventually results in death.³³ Over the past 30 years, the overall survival of the majority of patients with Grade II/III glioma has not significantly improved.² Thus, it is vital to develop an individualized treatment protocol for patients with Grade II/III glioma. Reliable prognostic biomarkers may identify patients with poor outcomes who might benefit from intensive therapy. Recently, significant research on immune-related gene (IRG) expression has shown great prognostic value in Grade II/III glioma. For example, Deng et al³⁴ found that 397 IRGs were significantly associated with Grade II/III glioma survival. Our risk model also showed great prognostic value in Grade II/III glioma. However, there is room for improvement in the effect of current gene-targeted therapy. Through our research, we are trying to find new gene targets. The immune-related gene pair signature depends on relative ranking and paired comparisons of gene expression profiles within a sample. This strategy is favorable for use with combined gene expression profiles obtained from multiple databases. Patients with Grade II/III glioma can be classified into subgroups with different survival outcomes through this prognostic immune signature. Therefore, our signature can be used to evaluate Grade II/III glioma prognosis in a single-sample, individualized form.

Our 8-IRGP signature consists of 15 immune-related genes that play vital roles in regulating the immune microenvironment. Bone morphogenetic protein 2 (MP2) is known to facilitate differentiation and growth inhibition in glioma through the downregulation of both O6-methylguanine-DNA methyltransferase (MGMT) and hypoxia-inducible factor-1a (HIF-1a).³⁵ It has been reported that killer cell lectin-like receptor subfamily C member 2 (KLRC2), a transmembrane receptor that activates natural killer cells, is expressed at higher levels in glioma-initiating cells than in neural stem cell lines and normal adult brain tissue.³⁶ Emerging evidence suggests that chemokine (C-X-C motif) ligand 12 (CXCL12), which is secreted from glioma stem cells (GSCs), plays vital roles in tumorigenesis and proliferation.³⁷ Additionally, many of the selected 26 functional pathways are related to the immune microenvironment, such as regulation of the immune effector process, positive regulation of cytokine production, and cell chemotaxis. Studies have shown that numerous cytokines and chemokines promote the infiltration of various cells, including circulating progenitor cells, endothelial cells, and a range of immune cells, such as peripheral macrophages, microglia, leukocytes, and CD4+ T cells, into gliomas.^38–40 These non-neoplastic cells play essential roles in tumor growth, metastasis, and response to treatment.¹⁴ Thus, these functional pathways are important in the glioma immune microenvironment.

There are also some limitations to this study. First, immune-related gene expression profiles were determined by RNA-seq, and they are difficult to generalize for use in daily clinical applications because of the sophisticated detection program needed and high cost. Second, our study was a retrospective analysis, so prospective research and some related basic experiments are needed to validate the results.

Conclusion

In summary, we systematically estimated the prognostic value of a new IRGP prognostic model, which may provide a legitimate approach to Grade II/III glioma management. Additionally, we identified 26 functional pathways from a GO analysis and GSEA, indicating that these immune-related gene pairs function through these pathways.