Research Paper Volume 14, Issue 10 pp 4357—4375

Glioblastoma: two immune subtypes under the surface of the cold tumor

Wu Xiong1,2, *, , Cong Li1,2, *, , Guang Kong2, *, , Bowen Wan3, , Siming Wang1,2, , Jin Fan1,2, ,

  • 1 Department of Orthopedics, The First Affiliated Hospital of Nanjing Medical University, Nanjing, Jiangsu, China
  • 2 Nanjing Medical University, Nanjing, Jiangsu, China
  • 3 Department of Orthopedics, Clinical Medical College of Yangzhou University, Subei People's Hospital of Jiangsu, Yangzhou, Jiangsu, China
* Equal contribution

Received: November 8, 2021       Accepted: April 22, 2022       Published: May 23, 2022      

https://doi.org/10.18632/aging.204067
How to Cite

Copyright: © 2022 Xiong et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Glioblastoma is classified as an immunocompromised tumor. The immune pattern beneath the cold tumor surface, however, has yet to be confirmed. Understanding the immune pattern of glioblastoma will aid in the development of effective treatment strategies. We performed weighted gene co-expression network analysis on all immune-related genes in TCGA-GBM transcriptional data and screened 35 prognosis-related immune genes. Unsupervised consistent clustering of these genes was used to analyze the immunological pattern of GBM. A glioblastoma immune prognostic score was developed by using 13 genes discovered by cox regression methods and verified with the GEO dataset to assess the immune profile, prognosis, and immunotherapy effects in individual patients. Glioblastoma has two immune modalities, immune tolerance and immunodeficiency, with distinct immune microenvironments, tumor-associated macrophages being one of the most promising new therapeutic targets. GIPS is a promising biomarker for assessing immune evasion mechanisms, immunotherapy responses, and prognosis in patients.

Introduction

Glioblastoma is a highly malignant tumor of the central nervous system [1], and patients’ prognoses remain dismal despite the availability of numerous treatments such as surgery, radiotherapy, and chemotherapy [2, 3]. Recently, immunotherapy has become increasingly popular, with immune checkpoint inhibitors serving as a new anti-tumor treatment [4, 5]. Immune system function can be suppressed to prevent autoimmune illnesses via immune checkpoints, but if they are suppressed too much, the body's ability to identify and eliminate aberrant cells is compromised [6]. Inhibitors of immunological checkpoints have anti-tumor effects through inhibiting the immune checkpoints' function, but their success in patients with glioblastoma has been uneven. A major cause of immunotherapy's poor efficacy is a lack of understanding of the immunological characteristics of glioblastoma. Individualized immunotherapy for glioblastoma can be achieved by analyzing each patient's immunological profile. However, glioblastoma’s immune milieu is poorly understood. It’s critical to clarify the immune trait of glioblastoma and develop marks that can accurately assess the immunological profile and prognosis of patients.

The goal of this study was to determine the immunological pattern of glioblastoma and create a marker that may be used to assess the prognosis and immune profile of GBM patients. WGCNA analysis of all immune-associated genes in TCGA-GBM transcriptional data was used to filter prognosis-associated immune genes for unsupervised clustering analysis, and we investigated the immunological traits of distinct immune-associated clusters. A Glioblastoma Immunological Prognostic Score (GIPS) was developed to measure individual individuals’ immune characteristics. We then investigated the GIPS’s molecular and immunological properties, examined its ability to predict prognosis, and associated it with immunotherapy response and immunotherapy.

Results

Identification of immune-related prognostic genes

The design of this study is shown in Supplementary Figure 1A. Differential analysis of immune-related genes from ImmPort and InnateDB yielded 802 differentially expressed immune-related genes. Supplementary Figure 1B and 1C show the results of the immune differential gene GO and KEGG pathway enrichment analyses. The candidate genes were subjected to WGCNA analysis in order to identify the immune-related center genes. Based on the scale-free network and the correlation coefficient of greater than 0.9, the soft-thresholding power was optimally set at (Supplementary Figure 1D and 1E). Four modules were identified based on the best soft threshold ability (Figure 1A and 1B). Among them, the blue module had the largest Pearson correlation coefficient and the smallest P value, so the genes in the blue module were selected for further analysis. The Univariate Cox analysis of the genes in the blue module revealed 35 prognostic-related immune genes with p ≤ 0.01 (Supplementary Figure 2). We then looked into the characteristics of the 35 genes in greater depth. These genes had a low mutation frequency but mostly copy number changes, with copy number loss dominating (Figure 1C and 1D). It's worth mentioning that the expression levels of most genes positively correlate with their copy number (Supplementary Figure 3). However, while CNV can explain many observed changes in prognostic-related immune genes expression, CNV is not the only factor involved in the regulation of mRNA expression. Furthermore, we identified BCL3 as the most common transcription factor for these prognostic immune genes and built a protein interaction network for them (Figure 1E and 1F).

Molecular characterizations of immune-related prognostic genes. (A) Weighted gene coexpression network analysis (WGCNA) of immune-related differentially expressed genes with a soft threshold β = 7. (B) Gene modules related to HNSCC obtained by WGCNA. (C) Mutation frequency of 35 immune-related prognostic genes. (D) CNV variation frequencies of 35 immune-related prognostic genes in the TCGA-GBM cohort. (E) Transcription factors for 35 immune-related prognostic genes. (F) The network of the 35 immune-related prognostic genes.

Figure 1. Molecular characterizations of immune-related prognostic genes. (A) Weighted gene coexpression network analysis (WGCNA) of immune-related differentially expressed genes with a soft threshold β = 7. (B) Gene modules related to HNSCC obtained by WGCNA. (C) Mutation frequency of 35 immune-related prognostic genes. (D) CNV variation frequencies of 35 immune-related prognostic genes in the TCGA-GBM cohort. (E) Transcription factors for 35 immune-related prognostic genes. (F) The network of the 35 immune-related prognostic genes.

Two immune patterns of glioblastoma

Unsupervised consensus clustering of the expression of the above prognosis-related immune genes was used to investigate the immune pattern of glioblastoma. The empirical cumulative distribution function (CDF) plots show the consensus distributions for k (1–9) (Supplementary Figure 4A and 4B). Given the consensus matrix for the analysis, k = 2 appeared to be the best option. The consensus matrix demonstrates that an unsupervised algorithm based on these genes is capable of clearly distinguishing samples and that each sample in the cluster has a high degree of correlation (Supplementary Figure 4C4E). As a result, we classified GBM patients into two groups, called immune clusters A and B, based on the expression of prognostically relevant immune genes. The two immune clusters had significantly different transcript expression profiles, according to principal component analysis (Supplementary Figure 4F).

GSVA enrichment analysis revealed 83 differential pathways between the two immune clusters, and the top 20 KEGG pathways are shown in Figure 2A. Immune-related pathways such as the nod-like receptor signaling pathway, the chemokine signaling pathway, and leukocyte transendothelial migration were significantly enriched in immune cluster B, whereas DNA replication and nucleotide excision repair were significantly enriched in immune cluster A.

TME Characterization of the two immune clusters. (A) GSVA enrichment analysis of the two immune clusters. (B) TME scores of the two immune clusters. (C) Differential analysis of immune cell abundance between two immune clusters. (D) Differential of immune function between the two immune clusters. (E) Differential analysis of expression of immune checkpoints in the two m6A clusters. (F) Kaplan-Meier OS analysis in the two immune clusters. P = 0.022. (*, **, and *** indicate p ≤ 0.05,

Figure 2. TME Characterization of the two immune clusters. (A) GSVA enrichment analysis of the two immune clusters. (B) TME scores of the two immune clusters. (C) Differential analysis of immune cell abundance between two immune clusters. (D) Differential of immune function between the two immune clusters. (E) Differential analysis of expression of immune checkpoints in the two m6A clusters. (F) Kaplan-Meier OS analysis in the two immune clusters. P = 0.022. (*, **, and *** indicate p ≤ 0.05, <0.01, and <0.001, respectively).

To reveal the immune landscape of the two immune clusters, we analyzed the immune microenvironmental features of the two immune clusters. Immune cluster B contained a higher number of both immune cells and stromal cells (Figure 2B). Further investigation revealed that almost all adaptive and innate immune cells were more infiltrated in immune cluster B, and immune cluster B's cellular function was more active (Figure 2C and 2D). Furthermore, the expression of immune checkpoints other than LAG3, CD160, and CD200 was more abundant in immune cluster B (Figure 2E). Because immune cluster B's stroma restricted immune cell entry into the tumor parenchyma and overexpressed immune checkpoints inhibited immune cell function, we classified immune cluster B as an immune tolerance phenotype. Immune cluster A had fewer immune cells and lower immune activity, indicating an immunodeficient phenotype. Furthermore, patients in immune cluster A had a better prognosis (Figure 2F).

Construction of GIPS

Thirteen immune genes with independent prognosis and their coefficients were identified based on multivariate analysis of 35 prognosis-associated immune genes (Supplementary Table 1). We then constructed GIPS based on the coefficients for assessing the immune status, immunotherapy response, and prognosis of individual patients. Univariate Cox regression analysis showed that age, IDH1 status and GIPS were significantly associated with the prognosis of GBM (Figure 3A). Multifactorial Cox regression analysis confirmed that GIPS was an independent prognostic factor after adjusting for other factors (Figure 3B).

Prognostic analysis of the GIPS subgroups. (A) Univariate Cox analysis of clinical factors and the GIPS. (B) Multivariate Cox analysis of the factors significant in the univariate Cox analysis. (C) K-M analysis of the GIPS subgroups in TCGA-GBM cohort (P D) K-M analysis of the GIPS subgroups in GEO cohort (P = 0.038).

Figure 3. Prognostic analysis of the GIPS subgroups. (A) Univariate Cox analysis of clinical factors and the GIPS. (B) Multivariate Cox analysis of the factors significant in the univariate Cox analysis. (C) K-M analysis of the GIPS subgroups in TCGA-GBM cohort (P < 0.001). (D) K-M analysis of the GIPS subgroups in GEO cohort (P = 0.038).

Taking the median GIPS as the cut-off value, low-GIPS patients had a better OS than high-GIPS patients (Figure 3C). Then, the role of GIPS was validated by using the GSE13041 dataset. As shown in Figure 3D, the patients in the low-GIPS subgroup had a significantly better prognosis than those in the high-GIPS subgroup, consistent with the results of the TCGA dataset.

Molecular analysis of various GIPS subgroups

GSEA was used to identify the set of genes that were enriched in various GIPS subgroups. The genomes of High-GIPS patients were enriched in cytokine-cytokine receptor interaction and the chemokine signaling pathway, whereas the genomes of Low-GIPS patients were enriched in the cell cycle and DNA replication (Figure 4A and 4B).

Molecular traits of distinct GIPS subgroups. (A) GSEA analysis in High-GIPS subgroup (P B) GSEA analysis in Low-GIPS subgroup (P C) Mutated genes (top 15) in High-GIPS subgroups. (D) Mutated genes (top 15) in Low-GIPS subgroups. (E) The correlation between expression level and mutations of genes (PTEN, TP53, TTN, IDH1 and EGFR).

Figure 4. Molecular traits of distinct GIPS subgroups. (A) GSEA analysis in High-GIPS subgroup (P < 0.05). (B) GSEA analysis in Low-GIPS subgroup (P < 0.05). (C) Mutated genes (top 15) in High-GIPS subgroups. (D) Mutated genes (top 15) in Low-GIPS subgroups. (E) The correlation between expression level and mutations of genes (PTEN, TP53, TTN, IDH1 and EGFR).

We then examined the gene mutations to gain a better understanding of the immunological nature of the GIPS subgroup from a biological standpoint. Figure 4C and 4D depict the top 15 genes in the GIPS subgroups with the highest mutation rates. The most common type of mutation was missense mutation, followed by frameshift deletion and nonsense mutation. PTEN, TP53, TTN, and EGF2 mutation rates were greater than 15% in both groups. PTEN, TP53, and TTN expression levels were not associated with mutations, whereas IDH1 mutations were associated with decreased expression and EGFR mutations with increased expression (Figure 4E). Notably, the mutation rate of IDH1 was 0% in the high GIPS subgroup and 10% in the low GIPS subgroup.

Immunological characteristics of various GIPS subgroups

The High-GIPS subgroup had higher immune and stroma scores (Figure 5A), and correlation analysis between GIPS and TME revealed that GIPS was significantly positively correlated with immune (R = 0.27, p = 2.2e-16) and stroma (R = 0.22, p = 0.0045, Figure 5B and 5C) scores. There were more innate and adaptive immune cells, EMT and Pan-F-TBRS in the High-GIPS subgroup, while there was more DNA damage repair, DNA replication and mismatch repair in Low-GIPS subgroup (Supplementary Figure 5A and 5B). Moreover, almost all common immune checkpoints had higher expression in the High-GIPS subgroup (Figure 5D). According to the alluvial diagram, the majority of low GIPS patients belong to immune cluster A (immune tolerance pattern), while the majority of high GIPS patients belong to immune cluster B (immune deficiency pattern) (Figure 5E).

TME Characterization of distinct GIPS subgroups. (A) TME scores of the distinct GIPS subgroups. (B) Spearman correlation analysis of GIPS scores with immune scores. R = 0.27, P ≤ 0.001. (C) Spearman correlation analysis of GIPS scores with stromal scores. R = 0.22, P ≤ 0.01. (D) Differential analysis of expression of immune checkpoints in the different GIPS subgroups. (E) Alluvial diagram of GBM patient immune cluster and GIPS.

Figure 5. TME Characterization of distinct GIPS subgroups. (A) TME scores of the distinct GIPS subgroups. (B) Spearman correlation analysis of GIPS scores with immune scores. R = 0.27, P ≤ 0.001. (C) Spearman correlation analysis of GIPS scores with stromal scores. R = 0.22, P ≤ 0.01. (D) Differential analysis of expression of immune checkpoints in the different GIPS subgroups. (E) Alluvial diagram of GBM patient immune cluster and GIPS.

Immunotherapy response in various GIPS subgroups

TIDE was used to assess immunotherapy response in various GIPS subgroups. Patients with high TIDE scores have a lower immunotherapy response, indicating that immunotherapy is less likely to benefit them. TIDE scores were higher in the High-GIPS subgroup than in the Low-GIPS subgroup, indicating that ICI therapy may benefit Low-GIPS patients more (Figure 6A). The effects of anti-CTLA4 therapy on different GIPS subgroups were then compared. The Low-GIPS subgroup did fare better in terms of treatment outcomes (Figure 6B).

Immunotherapy response in different GIPS subgroups. (A) Differential analysis of TIDE score between GIPS subgroups. (B) Effectiveness of immunotherapy in GIPS subgroups. (C) Differential analysis of dysfunction score between GIPS subgroups. (D) Differential analysis of MSI score between GIPS subgroups. (E) Differential analysis of T-cell exclusion score between GIPS subgroups. (F) ROC curves of GIPS for predicting 1-, 2-, and 3-year survival in TCGA GBM cohorts. (G) Graphical summary of this study.

Figure 6. Immunotherapy response in different GIPS subgroups. (A) Differential analysis of TIDE score between GIPS subgroups. (B) Effectiveness of immunotherapy in GIPS subgroups. (C) Differential analysis of dysfunction score between GIPS subgroups. (D) Differential analysis of MSI score between GIPS subgroups. (E) Differential analysis of T-cell exclusion score between GIPS subgroups. (F) ROC curves of GIPS for predicting 1-, 2-, and 3-year survival in TCGA GBM cohorts. (G) Graphical summary of this study.

To determine the reasons for the differences in immunotherapy response, we compared microsatellite instability (MSI), T-cell exclusion, and T-cell dysfunction scores in the two subgroups. Low-GIPS subgroup had a higher microsatellite instability (MSI) score, while the High-GIPS subgroup had a higher T-cell dysfunction score, but there was no difference in T-cell exclusion between the 2 subgroups (Figure 6C6E). Furthermore, ROC curve analysis showed that the AUC values of GIPS reached 0.771 at 1 year, 0.833 at 2 years, and 0.916 at 3 years (Figure 6F) and were higher than the AUC values of T-cell inflammatory signature (TIS) and TIDE (Supplementary Figure 5C5E).

Discussion

Glioblastoma is a type of central nervous system tumor that develops from glial stem cells [7]. Despite a variety of treatments, the prognosis for patients is frequently extremely poor [8]. Immunotherapy appears to offer a ray of hope for patients with glioblastoma. However, many patients do not respond well to immunotherapy and not only benefit less but also suffer from a slew of side effects [9]. According to some studies, the tumor microenvironment influences immunotherapy responses [10, 11]. Furthermore, there are no validated biomarkers to predict immunotherapy response and overall survival. As a result, in this study, we attempted to analyze the tumor microenvironment of glioblastoma using big data in order to elucidate the causes affecting immunotherapy response in glioblastoma patients and to propose corresponding treatment plans based on glioblastoma immune microenvironment characteristics.

Based on the GBM immune gene dataset, we identified 35 immune-related prognostic genes using WGCNA and K-M analysis in this study. Although the mutation rate was low, most of these genes had copy number loss. It has been demonstrated that gene copy number loss is linked to tumorigenesis [12, 13]; for example, the loss of ErbB4 receptor tyrosine kinase is linked to the development of glioblastoma [14]. Thus, immune gene copy number loss may play a role in the development of glioblastoma. Based on an unsupervised clustering analysis of these 35 genes, we discovered that glioblastoma has two immune patterns. Immune cluster A is primarily enriched in DNA replication and nucleic acid repair pathways, with less immune cell infiltration and low immune function, defining the immunodeficiency pattern. Immune cluster B had a high infiltration of immune cells and active immune function and was primarily enriched with immune-related pathways. As a result, in immunodeficiency pattern patients, improving immune function may be beneficial in terms of prognosis. Numerous studies have found that a dense infiltration of T cells, particularly cytotoxic [15, 16], predicts a favorable outcome. However, immune cluster B did not show the expected survival advantage. Further research revealed that increased stromal cells in immune cluster B limited the function of immune cells infiltrating the tumor parenchyma, while highly expressed immune checkpoints and increased Treg cells inhibited T cells' tumor killing function. As a result, immune cluster B displayed an immune tolerance pattern. Furthermore, soluble factors secreted by stromal cells can induce microenvironment-mediated drug resistance [17]. As a result, rational drug combinations capable of targeting both tumor cells and the microenvironment may be the key to overcoming therapeutic resistance [18, 19].

GIPS was created using multivariate cox analysis to assess the immune pattern, immunotherapy response, and prognosis of individual patients. When compared to the High-GIPS subgroup, the Low-GIPS subgroup had a higher survival rate and a longer survival time. The validity of the model was subsequently demonstrated by a validation cohort. We then investigated mutations in distinct GIPS subgroups to better understand the immunological landscape of these subgroups. As reported earlier, missense mutation, followed by frameshift deletion and nonsense mutation. Notably, the mutation rate of IDH1 reached 10% in the Low-GIPS subgroup while there were no mutations in the High-GIPS subgroup. According to some studies, the IDH1 mutation is an independent predictor of longer overall survival (OS) and progression-free survival (PFS) in GBM patients [20, 21]. The FAT1-ROS-HIF-1 signaling pathway is activated by the IDH1 mutation, which suppresses malignant glioma [22]. Thus, Low-GIPS patients with high IDH1 mutations have a better prognosis than High-GIPS patients who do not have IDH1 mutations, which is consistent with our survival findings.

Understanding the TME can aid in the discovery of new ways to treat GBM, as well as altering the TME to improve the efficacy of immunotherapy. Tumor microenvironments differ between the two GIPS subgroups. Multiple immune cells, such as natural killer cells, macrophages, and Treg cells, are more abundant in the High-GIPS subgroup. Tumor-associated macrophages contribute to patients' poor prognosis by promoting glioblastoma growth and angiogenesis [23]. Moreover, a large infiltration of Treg cells leads to immunosuppression [24]. As a result, a number of studies have proposed therapeutic approaches that target tumor-associated macrophages [23]. For example, promoting macrophage polarization in combination with immune check inhibitors aids in the treatment of glioblastoma [25]; and blocking macrophage-associated immunosuppression to regulate glioblastoma angiogenesis [26]. These approaches, however, need to be investigated further. The Low-GIPS subgroup has less immune cell infiltration but a greater ability to repair damage.

TIDE was then used to assess the immunotherapy responsiveness of different GIPS subgroups. Patients in the High-GIPS subgroup had higher immune cell infiltration and TIDE and T cell dysfunction scores than those in the Low-GIPS subgroup, suggesting that their lower immunotherapy response could be due to immune evasion caused by T cell dysfunction. The Low-GIPS subgroup, on the other hand, had higher MSI scores and lower TIDE scores than the High-GIPS subgroup, indicating that these patients had less immune evasion and more MSI. The high mutational load caused by MSI has been shown to make the tumor immunogenic and sensitive to immune checkpoint inhibitors. We directly compared the responses of different GIPS subgroups to anti-CTLA4 therapy to further validate GIPS's ability to predict patients’ responses to immunotherapy. Patients in the Low-GIPS subgroup were found to be indeed more sensitive to immunotherapy.

The TIDE score is an algorithm that simulates tumor immune evasion mechanisms [27] and can be used to assess immunotherapy response in patients with a wide range of tumors, including bladder cancer, lung adenocarcinoma, and melanoma [28, 29]. Furthermore, the Tumor Inflammation Signature (TIS) provides quantitative and qualitative information about TME, which has been shown in a pan-cancer cohort to correlate with benefit from anti-PD -1 therapy [30]. Nevertheless, TIDE and TIS are both concerned with assessing the function and condition of T cells and do not truly reflect the impact of the tumor microenvironment on immunotherapeutic responses [2, 31], such as the important role that tumor-associated macrophages play in glioblastoma. Furthermore, both signatures are concerned with the patient's response to immunotherapy rather than the patient's overall survival, and life expectancy is an important consideration when making a clinical decision. GIPS has a higher predictive value than TIDE and TIS, and it may be a more valuable predictor for assessing a patient's immune condition, immunotherapy response, and prognosis.

In conclusion, glioblastoma has two immune modalities, immune tolerance and immunodeficiency, with distinct immune microenvironments, tumor-associated macrophages being one of the most promising new therapeutic targets. GIPS is a promising biomarker for assessing immune evasion mechanisms, immunotherapy responses, and prognosis in patients, but more research is needed to confirm its utility (Figure 6G).

Methods

Data collection and processing

The TCGA-GBM dataset and GSE13041 were used to obtain GBM transcriptome data. GTEx was used to obtain RNAseq transcriptome data for healthy human tissues. The GTEx and TCGA datasets were combined and reconciled using quantile normalization and batch effect removal using svaseq. Immune-related gene lists were obtained from the databases ImmPort and InnateDB.

Identification of immune-related prognostic genes

Differential analysis of immune-related genes from ImmPort and InnateDB yielded 802 differentially expressed immune-related genes. Then, WGCNA was performed to identify these genes. Based on the scale-free network and the correlation coefficient of greater than 0.9, the soft-thresholding power was optimally set at 7. Univariate Cox analysis of the genes in the blue module revealed 35 prognostic-related immune genes (p ≤ 0.01).

Unsupervised consensus clustering based on immune-related prognostic genes

According to expression levels of the 35 prognostic-related immune genes, unsupervised clustering analysis was used for identification of various immune patterns and patient classification for further analyses. Cluster numbers and their stabilities were evaluated using a consistent clustering algorithm.

Gene set variation analysis

Gene set “c2.cp.kegg.v6.2.symbols” was obtained from the MSigDB and GSVA enrichment analysis used to identify different immune patterns using “GSVA” package on R, with adjusted p ≤ 0.05 indicating statistical significance. In non-parametric and unsupervised approaches, GSVA is used for estimation of pathway and biological process changes in gene expression datasets.

Estimation of TME immune trait

The ssGSEA enrichment fraction was used to calculate the relative abundances of each TME-infiltrating cell per sample. Charoentong's study, which annotated human immune cell subtypes, immune checkpoints, and EMT, provided the genomes used to mark each TME-infiltrating immune cell type.

Establishment and subsequent validation of the GIPS

Among 35 immune-related prognostic genes, multifactorial Cox regression analysis was used to identify genes with significant effects on OS. The GIPS for each sample was calculated by multiplying the expression values of specific genes by their Cox model weights and then summing them. Log-rank tests were used to assess GIPS's prognostic ability in the TCGA and GEO cohorts using Kaplan-Meier (K-M) survival curves. GIPS’s independent prognostic value was validated using univariate and multivariate Cox regression analyses.

Statistical analysis

The limma R package was used to analyze differential gene expression. The statistical difference between the two groups was calculated using the Wilcoxon rank sum test. For comparisons of more than two groups, the Kruskal-Wallis test was used. For all statistical analyses, R software was used.

Ethics approval and consent to participate

Patient information is available in public databases that was collected with patients’ informed consent.

Consent for publication

All authors give consent for the publication of this manuscript in Aging -us.

Availability of data and materials

The data used in the study are described in detail in “Data collection and processing”.

Abbreviations

GBM: glioblastoma; TME: tumor microenvironment; WGCNA: Weighted gene coexpression network analysis; GIPS: glioblastoma immune prognostic score.

Author Contributions

This piece was created by JF. WX and CL combined and analyzed the data. This manuscript was written by GK, BW and SW. This manuscript was approved by all of the authors.

Acknowledgments

Thanks to TCGA, GEO, GTEx, ImmPort and InnateDB database for their work.

Conflicts of Interest

The authors declare no conflicts of interest related to this study.

Funding

This study was funded by the Natural Science Foundation of China (81772352), Program for Changjiang Scholars and Innovative Research Team in university (KY216R202103), Postgraduate Research & Practice Innovation Program of Jiangsu Province (KYCX21_1613, KYCX20_1434).

References

  • 1. Omuro A, DeAngelis LM. Glioblastoma and other malignant gliomas: a clinical review. JAMA. 2013; 310:1842–50. https://doi.org/10.1001/jama.2013.280319 [PubMed]
  • 2. Touat M, Idbaih A, Sanson M, Ligon KL. Glioblastoma targeted therapy: updated approaches from recent biological insights. Ann Oncol. 2017; 28:1457–72. https://doi.org/10.1093/annonc/mdx106 [PubMed]
  • 3. Le Rhun E, Preusser M, Roth P, Reardon DA, van den Bent M, Wen P, Reifenberger G, Weller M. Molecular targeted therapy of glioblastoma. Cancer Treat Rev. 2019; 80:101896. https://doi.org/10.1016/j.ctrv.2019.101896 [PubMed]
  • 4. Bagchi S, Yuan R, Engleman EG. Immune Checkpoint Inhibitors for the Treatment of Cancer: Clinical Impact and Mechanisms of Response and Resistance. Annu Rev Pathol. 2021; 16:223–49. https://doi.org/10.1146/annurev-pathol-042020-042741 [PubMed]
  • 5. Kennedy LB, Salama AKS. A review of cancer immunotherapy toxicity. CA Cancer J Clin. 2020; 70:86–104. https://doi.org/10.3322/caac.21596 [PubMed]
  • 6. Chamoto K, Al-Habsi M, Honjo T. Role of PD-1 in Immunity and Diseases. Curr Top Microbiol Immunol. 2017; 410:75–97. https://doi.org/10.1007/82_2017_67 [PubMed]
  • 7. Bhaduri A, Di Lullo E, Jung D, Müller S, Crouch EE, Espinosa CS, Ozawa T, Alvarado B, Spatazza J, Cadwell CR, Wilkins G, Velmeshev D, Liu SJ, et al. Outer Radial Glia-like Cancer Stem Cells Contribute to Heterogeneity of Glioblastoma. Cell Stem Cell. 2020; 26:48–63.e6. https://doi.org/10.1016/j.stem.2019.11.015 [PubMed]
  • 8. Braun K, Ahluwalia MS. Treatment of Glioblastoma in Older Adults. Curr Oncol Rep. 2017; 19:81. https://doi.org/10.1007/s11912-017-0644-z [PubMed]
  • 9. Lim M, Xia Y, Bettegowda C, Weller M. Current state of immunotherapy for glioblastoma. Nat Rev Clin Oncol. 2018; 15:422–42. https://doi.org/10.1038/s41571-018-0003-5 [PubMed]
  • 10. Bader JE, Voss K, Rathmell JC. Targeting Metabolism to Improve the Tumor Microenvironment for Cancer Immunotherapy. Mol Cell. 2020; 78:1019–33. https://doi.org/10.1016/j.molcel.2020.05.034 [PubMed]
  • 11. Frankel T, Lanfranca MP, Zou W. The Role of Tumor Microenvironment in Cancer Immunotherapy. Adv Exp Med Biol. 2017; 1036:51–64. https://doi.org/10.1007/978-3-319-67577-0_4 [PubMed]
  • 12. Mansisidor A, Molinar T Jr, Srivastava P, Dartis DD, Pino Delgado A, Blitzblau HG, Klein H, Hochwagen A. Genomic Copy-Number Loss Is Rescued by Self-Limiting Production of DNA Circles. Mol Cell. 2018; 72:583–93.e4. https://doi.org/10.1016/j.molcel.2018.08.036 [PubMed]
  • 13. Paolella BR, Gibson WJ, Urbanski LM, Alberta JA, Zack TI, Bandopadhayay P, Nichols CA, Agarwalla PK, Brown MS, Lamothe R, Yu Y, Choi PS, Obeng EA, et al. Copy-number and gene dependency analysis reveals partial copy loss of wild-type SF3B1 as a novel cancer vulnerability. Elife. 2017; 6:e23268. https://doi.org/10.7554/eLife.23268 [PubMed]
  • 14. Jones DC, Scanteianu A, DiStefano M, Bouhaddou M, Birtwistle MR. Analysis of copy number loss of the ErbB4 receptor tyrosine kinase in glioblastoma. PLoS One. 2018; 13:e0190664. https://doi.org/10.1371/journal.pone.0190664 [PubMed]
  • 15. He M, Wang Y, Zhang G, Cao K, Yang M, Liu H. The prognostic significance of tumor-infiltrating lymphocytes in cervical cancer. J Gynecol Oncol. 2021; 32:e32. https://doi.org/10.3802/jgo.2021.32.e32 [PubMed]
  • 16. Tian Y, Wei Y, Liu H, Shang H, Xu Y, Wu T, Liu W, Huang A, Dang Q, Sun Y. Significance of CD8+ T cell infiltration-related biomarkers and the corresponding prediction model for the prognosis of kidney renal clear cell carcinoma. Aging (Albany NY). 2021; 13:22912–33. https://doi.org/10.18632/aging.203584 [PubMed]
  • 17. Wu T, Dai Y. Tumor microenvironment and therapeutic response. Cancer Lett. 2017; 387:61–8. https://doi.org/10.1016/j.canlet.2016.01.043 [PubMed]
  • 18. Valkenburg KC, de Groot AE, Pienta KJ. Targeting the tumour stroma to improve cancer therapy. Nat Rev Clin Oncol. 2018; 15:366–81. https://doi.org/10.1038/s41571-018-0007-1 [PubMed]
  • 19. Matsumura Y. Cancer stromal targeting therapy to overcome the pitfall of EPR effect. Adv Drug Deliv Rev. 2020; 154–155:142–50. https://doi.org/10.1016/j.addr.2020.07.003 [PubMed]
  • 20. Khan I, Waqas M, Shamim MS. Prognostic significance of IDH 1 mutation in patients with glioblastoma multiforme. J Pak Med Assoc. 2017; 67:816–7. [PubMed]
  • 21. Yan H, Parsons DW, Jin G, McLendon R, Rasheed BA, Yuan W, Kos I, Batinic-Haberle I, Jones S, Riggins GJ, Friedman H, Friedman A, Reardon D, et al. IDH1 and IDH2 mutations in gliomas. N Engl J Med. 2009; 360:765–73. https://doi.org/10.1056/NEJMoa0808710 [PubMed]
  • 22. Li LC, Zhang M, Feng YK, Wang XJ. IDH1-R132H Suppresses Glioblastoma Malignancy through FAT1-ROS-HIF-1α Signaling. Neurol India. 2020; 68:1050–8. [PubMed]
  • 23. Muir M, Gopakumar S, Traylor J, Lee S, Rao G. Glioblastoma multiforme: novel therapeutic targets. Expert Opin Ther Targets. 2020; 24:605–14. https://doi.org/10.1080/14728222.2020.1762568 [PubMed]
  • 24. Lee-Chang C, Rashidi A, Miska J, Zhang P, Pituch KC, Hou D, Xiao T, Fischietti M, Kang SJ, Appin CL, Horbinski C, Platanias LC, Lopez-Rosas A, et al. Myeloid-Derived Suppressive Cells Promote B cell-Mediated Immunosuppression via Transfer of PD-L1 in Glioblastoma. Cancer Immunol Res. 2019; 7:1928–43. https://doi.org/10.1158/2326-6066.CIR-19-0240 [PubMed]
  • 25. Saha D, Martuza RL, Rabkin SD. Macrophage Polarization Contributes to Glioblastoma Eradication by Combination Immunovirotherapy and Immune Checkpoint Blockade. Cancer Cell. 2017; 32:253–67.e5. https://doi.org/10.1016/j.ccell.2017.07.006 [PubMed]
  • 26. Cui X, Morales RT, Qian W, Wang H, Gagner JP, Dolgalev I, Placantonakis D, Zagzag D, Cimmino L, Snuderl M, Lam RHW, Chen W. Hacking macrophage-associated immunosuppression for regulating glioblastoma angiogenesis. Biomaterials. 2018; 161:164–78. https://doi.org/10.1016/j.biomaterials.2018.01.053 [PubMed]
  • 27. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, Li Z, Traugh N, Bu X, Li B, Liu J, Freeman GJ, Brown MA, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018; 24:1550–8. https://doi.org/10.1038/s41591-018-0136-1 [PubMed]
  • 28. Chen X, Xu R, He D, Zhang Y, Chen H, Zhu Y, Cheng Y, Liu R, Zhu R, Gong L, Xiao M, Wang Z, Deng L, Cao K. CD8+ T effector and immune checkpoint signatures predict prognosis and responsiveness to immunotherapy in bladder cancer. Oncogene. 2021; 40:6223–34. https://doi.org/10.1038/s41388-021-02019-6 [PubMed]
  • 29. Wang Q, Li M, Yang M, Yang Y, Song F, Zhang W, Li X, Chen K. Analysis of immune-related signatures of lung adenocarcinoma identified two distinct subtypes: implications for immune checkpoint blockade therapy. Aging (Albany NY). 2020; 12:3312–39. https://doi.org/10.18632/aging.102814 [PubMed]
  • 30. Danaher P, Warren S, Lu R, Samayoa J, Sullivan A, Pekker I, Wallden B, Marincola FM, Cesano A. Pan-cancer adaptive immune resistance as defined by the Tumor Inflammation Signature (TIS): results from The Cancer Genome Atlas (TCGA). J Immunother Cancer. 2018; 6:63. https://doi.org/10.1186/s40425-018-0367-1 [PubMed]
  • 31. Chen Y, Li ZY, Zhou GQ, Sun Y. An Immune-Related Gene Prognostic Index for Head and Neck Squamous Cell Carcinoma. Clin Cancer Res. 2021; 27:330–41. https://doi.org/10.1158/1078-0432.CCR-20-2166 [PubMed]