A novel chemokine-based signature for prediction of prognosis and therapeutic response in glioma

Abstract

Aims: Gliomas are the primary malignant brain tumor and characterized as the striking cellular heterogeneity and intricate tumor microenvironment (TME), where chemokines regulate immune cell trafficking by shaping local networks. This study aimed to construct a chemokine-based gene signature to evaluate the prognosis and therapeutic response in glioma.

Methods: In this study, 1024 patients (699 from TCGA and 325 from CGGA database) with clinicopathological information and mRNA sequencing data were enrolled. A chemokine gene signature was constructed by combining LASSO and SVM-RFE algorithm. GO, KEGG, and GSVA analyses were performed for function annotations of the chemokine signature. Candidate mRNAs were subsequently verified through qRT-PCR in an independent cohort including 28 glioma samples. Then, through immunohistochemical staining (IHC), we detected the expression of immunosuppressive markers and explore the role of this gene signature in immunotherapy for glioma. Lastly, the Genomics of Drug Sensitivity in Cancer (GDSC) were leveraged to predict the potential drug related to the gene signature in glioma.

Results: A constructed chemokine gene signature was significantly associated with poorer survival, especially in glioblastoma, IDH wildtype. It also played an independent prognostic factor in both datasets. Moreover, biological function annotations of the predictive signature indicated the gene signature was positively associated with immune-relevant pathways, and the immunosuppressive protein expressions (PD-L1, IBA1, TMEM119, CD68, CSF1R, and TGFB1) were enriched in the high-risk group. In an immunotherapy of glioblastoma cohort, we confirmed the chemokine signature showed a good predictor for patients' response. Lastly, we predicted twelve potential agents for glioma patients with higher riskscore.

Conclusion: In all, our results highlighted a potential 4-chemokine signature for predicting prognosis in glioma and reflected the intricate immune landscape in glioma. It also threw light on integrating tailored risk stratification with precision therapy for glioblastoma.

Keywords: chemokine; glioma; immunotherapy; prognostic signature; tumor microenvironment.

FIGURE 1

Feature selection of a chemokine‐based gene signature. (A) LASSO coefficient profiles of the remained 27 chemokines. (B) The accuracy of the estimate generation for the SVM‐RFE algorithm. (C) The intersection feature selection between LASSO and SVM‐RFE algorithms. (D) The hazard ratio and p‐value of genes involved in multivariate Cox regression in the TCGA dataset

FIGURE 2

Landscape of clinical and molecular characteristics associated with the gene signature in gliomas. (A) TCGA (top) and CGGA (bottom) were arranged in an increasing order of the 4‐chemokine‐based riskscore. The relationship between the riskscore and patients' clinic pathological characteristics was evaluated (A, Spearman correlation tests; B, One‐way ANOVA test; C, Wilcoxon test). (B and C) The riskscore distribution between MGMT promoter methylated group and MGMT un‐methylated group in both TCGA (B) and CGGA (C) dataset. (D–F) The riskscore distribution in the different subgroup of glioma in a 28 independent cohort (Wilcoxon test and One‐way ANOVA test)

FIGURE 3

Kaplan–Meier survival analysis stratified by 4‐chemokine gene signature and the nomogram for survival prediction in TCGA dataset. (A) The riskscore of 4‐chemokine‐based signature distribution and survival status distribution for glioma patients. (B–F) Kaplan–Meier survival curves were plotted to estimate the overall survival probabilities in all grade's gliomas (B), GBM (C), IDH wildtype (D), IDHmut/1p19q intact (E) and IDHmut/1p19q codel (F). (G) The nomogram prediction of glioma patients for 1‐, 3‐, and 5‐year OS combining the signatures with clinic pathological features. (H) Calibration curves used to compare the predicted nomogram, the dashed diagonal line represents the ideal nomogram. (I) qRT‐PCR analysis of the four chemokine genes (CCL2, CCL5, CCL18, and CXCL16) expression between high‐ and low‐risk group in an independent validation cohort, 18S was used as an internal reference. (J) Kaplan–Meier survival curves were plotted to estimate the overall survival probabilities in an independent validation gliomas group

FIGURE 4

Distinct genomic alterations between high‐ and low‐risk group. (A) Differential somatic mutations were detected between high‐ and low‐risk group. (B) Top 20 significantly differential mutational genes (Fisher's exact test). (C) TMB between high‐ and low‐risk group (Wilcoxon test). (D and E) Distinct CNA profiles between gliomas in high‐ and low‐risk groups

FIGURE 5

Biological processes and signal pathways associated with the 4‐chemokine signature in TCGA dataset. (A) Correlation between chemokines‐based prognostic signature and transcriptomic expression profiles. (B) Biological processes enrichment of positively associated genes in TCGA dataset. (C) Enriched gene sets in HALLMARK collection in TCGA dataset. (D) Heatmap of signaling pathway activity scores by PROGENy (Wilcoxon test)

FIGURE 6

Immune cell infiltration and inflammatory profiles of the gene signature in TCGA dataset. (A) Heatmap of adaptive and innate immune cell types in high‐ and low‐risk group. (B) Heatmap of the MHC‐, costimulation‐, and inflammatory‐related genes expression in glioma patients from high‐ and low‐risk group. (C) The representative GAM related gene expression level between high‐ and low‐risk group in the TCGA dataset. (D) The representative IHC images of IBA1 and the correlation plot between riskscore and IBA1 protein expression. (E) The representative IHC images of TMEM119 and the correlation plot between riskscore and TMEM119 protein expression. (F) Sankey plot shows the relationship between glioma patients stratified by riskscore and 6 immune subtypes defiend by Thorsson et al. (G) The relationship between riskscore and inflammatory activity in glioma. (Wilcoxon test)

FIGURE 7

Distribution of immune‐response markers and the prognostic value of riskscore in patients with anti‐PD‐1 therapy. (A–D) Distribution of estimate score, stroma score, immune score and tumor purity in the TCGA dataset. (E) Distribution of different immune checkpoint mRNA expression in high‐ and low‐risk groups. (F) The correlation of riskscore and PD‐L1 protein expression in TCPA. (G) The distribution of PD‐L1 protein expression in high‐ and low‐risk groups by IHC staining. Scale bar, 50 μm. (H–J) ROC analysis of riskscore (H), TIS (I), and TIDE (J) score on overall survival at 6‐, 12‐, and 18‐month follow‐up in GBM patients receiving anti‐PD‐1 therapy

FIGURE 8

Identification of candidate agents with higher drug sensitivity in high‐risk group patients. (A and C) The dot chart showed the correlation coefficient between riskscore and the estimate half maximal inhibitory concentration (IC50) of 12 candidate drugs in TCGA (A) and CGGA (C) datasets. (B and D) The estimated IC50 of 12 candidate drugs were compared between high‐ and low‐risk groups in TCGA (B) and CGGA (D) datasets (Wilcoxon test)