A four-gene signature-derived risk score for glioblastoma: prospects for prognostic and response predictive analyses

Abstract

Objective: Glioblastoma (GBM) is the most common primary malignant brain tumor regulated by numerous genes, with poor survival outcomes and unsatisfactory response to therapy. Therefore, a robust, multi-gene signature-derived model is required to predict the prognosis and treatment response in GBM.

Methods: Gene expression data of GBM from TCGA and GEO datasets were used to identify differentially expressed genes (DEGs) through DESeq2 or LIMMA methods. The DEGs were then overlapped and used for survival analysis by univariate and multivariate COX regression. Based on the gene signature of multiple survival-associated DEGs, a risk score model was established, and its prognostic and predictive role was estimated through Kaplan-Meier analysis and log-rank test. Gene set enrichment analysis (GSEA) was conducted to explore high-risk score-associated pathways. Western blot was used for protein detection.

Results: Four survival-associated DEGs of GBM were identified: OSMR, HOXC10, SCARA3, and SLC39A10. The four-gene signature-derived risk score was higher in GBM than in normal brain tissues. GBM patients with a high-risk score had poor survival outcomes. The high-risk group treated with temozolomide chemotherapy or radiotherapy survived for a shorter duration than the low-risk group. GSEA showed that the high-risk score was enriched with pathways such as vasculature development and cell adhesion. Western blot confirmed that the proteins of these four genes were differentially expressed in GBM cells.

Conclusions: The four-gene signature-derived risk score functions well in predicting the prognosis and treatment response in GBM and will be useful for guiding therapeutic strategies for GBM patients.

Keywords: Differentially expressed genes; gene set enrichment analysis; glioblastoma prognosis; radiotherapy; temozolomide chemotherapy.

1

Differentially expressed genes between GBM and normal brain tissues. (A) Volcano plots showing the log₂ (fold change) of mRNA in GBM compared to normal brain tissues, and the corresponding–log₁₀ (adjusted P value) in TCGA, GSE4290 and GSE16011 datasets. Genes with adjusted P value below 0.05 and fold change above 2 (below -2) were marked with red (blue) dots. (B) Venn diagrams showing the gene numbers of the upregulated genes (left) and the downregulated genes (right) of GBM in TCGA, GSE4290 and GSE16011 cohorts. (C) Heatmaps of the overlapped genes in TCGA, GSE4290 and GSE16011 datasets. The overlapped genes that are from Figure 1B include 483 upregulated genes and 765 downregulated genes.

2

The four-gene signature-derived risk score is high in GBM. (A) The calculation formula and the value of risk score (top), the corresponding expression of four key genes (middle), and the associated clinicopathological parameters (bottom). “β” indicates the regression coefficient derived from multivariate COX stepwise regression in TCGA GBM cohort; “E” represents the expression value of the corresponding gene. The P-value indicates the correlation between the risk score and the clinicopathological parameters (Supplementary Table S5). KPS, Karnofsky performance status; IDH, isocitrate dehydrogenase; MGMT, O(6)-methylguanine-DNA methyltransferase; NA, not available. (B-D) Risk scores in GBM and normal brain tissues (B), in IDH-wt GBM and IDH-mut GBM (C), and in mesenchymal GBM and non-mesenchymal (proneural and classical) GBM (D) in TCGA GBM cohort. IDH, isocitrate dehydrogenase; IDH-wt, IDH-wild type; IDH-mut, IDH-mutation. Data are shown as mean ± SEM, *P < 0.05, *** P < 0.001, ns, not significant.

S1

Survival analysis of GBM patients. (A-B) Kaplan-Meier overall survival analysis of GBM patients stratified by IDH status (A), and expression subtypes (B). IDH, isocitrate dehydrogenase; IDH-wt, IDH-wild type; IDH-mut, IDH-mutation. Mes, mesenchymal; Non-Mes, non-mesenchymal, including proneural and classical subtypes.

3

Risk score predicts the prognosis and treatment response in TCGA GBM cohort. (A-C) Kaplan–Meier overall survival analysis among TCGA GBM patients stratified by risk score only (A), combined with IDH status (B), and combined with expression subtypes (C). Mes, mesenchymal; Non-Mes, non-mesenchymal, including proneural and classical subtypes. (D-E) Kaplan–Meier overall survival analysis of TCGA GBM patients with TMZ chemotherapy (D), or radiotherapy (E) according to the risk score. TMZ, temozolomide.

4

Performance of risk score in predicting the survival and treatment response in independent GEO GBM cohorts. (A) Comparison of risk scores between GBM and normal brain tissues in GSE4290, GSE16011, GSE59612 and GSE90604. (B) Kaplan-Meier overall survival analysis of high-risk and low-risk GBM patients in GSE16011, GSE43378 and GSE83300. (C) Kaplan-Meier overall survival analysis of GBM patients with radiotherapy in GSE16011 according to the risk score. (D) Risk scores in mesenchymal GBM and non-mesenchymal (proneural and classical) GBM (up), and Kaplan–Meier overall survival analysis of GBM patients stratified by risk scores combined with expression subtypes (down). Mes, mesenchymal; Non-Mes, non-mesenchymal, including proneural and classical subtypes. (E) Risk scores in IDH-wt GBM and IDH-mut GBM (up), and Kaplan–Meier overall survival analysis of GBM patients stratified by risk scores combined with IDH status (down). IDH, isocitrate dehydrogenase; IDH-wt, IDH-wild type; IDH-mut, IDH-mutation. Data are shown as mean ± SEM, *P < 0.05, ** P < 0.01, *** P < 0.001, ns, not significant.

5

Four-gene signature associated biological pathways in high-risk GBM. (A) Enrichment map showing the pathways enriched in high-risk GBM through GSEA analysis. Nodes represent enriched gene sets with P-value below 0.05. Node with deep red color correlates with small P-value. Node size corresponds to the number of genes within gene set. Edge thickness corresponds to the number of shared genes between gene sets. (B) Representative enriched pathways in high-risk GBM through GSEA analysis. NES, normalized enrichment score.

6

Validation of the differential expression of the four genes in GBM cells. (A and B) Western blot images (A) and the relevant quantification (B) of OSMR, HOXC10, SCARA3 and SLC39A10 in GBM cell line LN229, primary GBM cells (GBM-1 and GBM-2), and normal glial cell line HEB. The relative expression of target proteins is quantified in comparison with β-Actin and normalized to the corresponding expression in HEB cells. Data are shown as mean ± SEM from three independent experiments, *P < 0.05, ** P < 0.01, *** P < 0.001.