Development and validation of a glycolysis-associated gene signature for predicting the prognosis, immune landscape, and drug sensitivity in bladder cancer

Abstract

Background: Bladder cancer (BCa) is one of the most common malignancies worldwide, and its prognostication and treatment remains challenging. The fast growth of various cancer cells requires reprogramming of its energy metabolism using aerobic glycolysis as a major energy source. However, the prognostic and therapeutic value of glycolysis-related genes in BCa remains to be determined.

Methods: The fused merge dateset from TCGA, GSE13507 and GSE31684 were used for the analysis of glycolysis-related genes expression or subtyping; and corresponding clinical data of these BCa patients were also collected. In the merge cohort, we constructed a 18 multigene signature using the least absolute shrinkage and selection operator (LASSO) Cox regression model. The four external cohorts (i.e., IMvigor210, GSE32894, GSE48276 and GSE48075) of BCa patients were used to validate the accuracy. We evaluated immune infiltration using seven published algorithms: CIBERSORT, QUANTISEQ, XCELL, TIMER, CIBERSORT-ABS, EPIC, and MCPCOUNTER. Subsequently, in order to analyze the correlation between risk groups(scores) and overall survival, recognised immunoregolatory cells or common chemotherapeutic agents, clinicopathological data and immune checkpoint-related genes of BCa patients, Wilcox rank test, chi-square test, cox regression and spearman's correlation were performed.

Results: Conspicuously, we could see that CD8+ T, cancer associated fibroblast, macrophage M2, NK, endothelial cells and so on were significantly dysregulated between the two risk groups. In addition, compared with the low-risk group, high-risk group predicted poor prognosis and relatively weak sensitivity of chemotherapy. Additionally, we also found that the expression level of partial genes in the model was significantly correlated with objective responses to anti-PD-1 or anti-PD-L1 treatment in the IMvigor210, GSE111636, GSE176307, GSE78220 or GSE67501 cohort; and its expression level was also varied in different objective response cases receiving tislelizumab combined with low-dose nab-paclitaxel therapy based on our mRNA sequencing (TRUCE-01). According to "GSEA" algorithm of R package "clusterProfiler", the most significantly enriched HALLMARK, KEGG pathway and GO term was separately the 'Epithelial Mesenchymal Transition', 'Ecm Receptor Interaction' and 'MF_Extracellular_matrix_structural_constitunet' in the high- vs. low-risk group. Subsequently, we verified the protein and mRNA expression of interested model-related genes from the Human Protein Atlas (HPA) and 10 paired BCa tissues collected by us. Furthermore, in vitro functional experiments demonstrated that FASN was a functional oncogene in BCa cells through promoting cell proliferation, migration, and invasion abilities.

Conclusion: In summary, the glycolysis-associated gene signature established by us exhibited a high predictive performance for the prognosis, immunotherapeutic responsiveness, and chemotherapeutic sensitivity of BCa. And, The model also might function as a chemotherapy and immune checkpoint inhibitor (ICI) treatment guidance.

Keywords: bladder cancer; gene signature; glycolysis; immune status; prognosis.

Figure 1

Flow diagram of the study.

Figure 2

Identification the glycolysis-associated molecular subtypes of BCa. (A) Cumulative distribution function (CDF) curves of the consensus scores for different subtype numbers (k = 2-9). (B) CDF delta area curve with k = 2 to 10. (C) Heatmaps of consensus matrices for k = 4, 5, 6 based on CDF curve estimation. (D) Principal component analysis (PCA) for six molecular typing. A single point represents each sample, and each subtype is represented by a different color. (E-H) KM survival analysis of OS, DSS, PFS, and DFS for the six subtypes.

Figure 3

Heatmap of 34 glycolysis-related prognosis DEGs, and the chi-square analysis between corresponding the six molecular subtypes and commonly clinical traits. *p < 0.05, ***p < 0.001.

Figure 4

Stratification of overall, training and self-validation set from the fused merge dateset. Kaplan-Meier survival plots, ROC curves for 1-, 3-, 5-years along with the c-index, calibration curve, the distribution of risk score ranking and survival status, model gene expression heatmaps with order of increasing risk score, and PCA analyses in the all (A), training (B) and self-validation set (C).

Figure 5

Establishment and evaluation of the nomogram. (A, B) The risk score was independent risk factors for BCa by univariate and multivariate Cox regression analyses of overall survival. (C) Construction of a glycolysis-related gene signature combined with clinical features nomogram for predicting the 1-, 3- and 5-year OS rates. (D) Calibration curve of nomogram. (E) A multi-index time-dependent ROC analysis was employed to evaluate the predictive accuracy of the glycolysis-related gene signature or our nomogram, and compare it with other clinical traits. (F) DCA of 1-, 3- and 5-year was applied to render clinical validity to the constructed gene signature or nomograms. (G–I) The most significant HALLMARK, KEGG pathways and GO functional enrichment in the high risk and low risk groups by GSEA method are displayed. *p < 0.01, ***p < 0.001.

Figure 6

External validation of the 18 model genes signature in the IMvigor210 cohort by different aspects such as (A) survival outcomes, (B) ROC curves, (C) calibration plots, (D) the distribution of risk scores, (E) the distribution of overall outcomes and (F) the gene expression heatmap analysis. (G, H) Correlation between riskscore group and clinicopathological data of BCa patients via Wilcox rank test or Chi-square test. CR, complete response; PR, partial response; SD, stable disease; PD, progressive disease.

Figure 7

External identification of the prognostic model for bladder cancer based on (A) GSE32894, (B) GSE48276, and (C) GSE48075 datasets.

Figure 8

(A) Sankey diagram of glycolysis-related gene cluster distribution in groups with different risk score and survival outcomes. (B) The boxplot showing the relationship between the 6 clusters and risk score. (C) Kaplan–Meier curves for patients stratified by both TMB and risk score. (D) The difference in the expression of immune cell between high and low risk score. (E, F) The difference in the expression of immune checkpoint-related genes and HLA genes between high-risk and low-risk groups. (G) The modeled potential application values of ICBs for BCa. *p < 0.05, **p < 0.01, ***p < 0.001.

Figure 9

The expression level of the modeled genes in groups with a different immunotherapy response status. (A-F) These modeled genes expression in different cohorts, including GSE111636, Imvigor210, GSE176307, GSE78220, our mRNA sequencing (TRUCE-01), and GSE67501.

Figure 10

The changes of partially modeled genes expression before and after tislelizumab combined with low-dose nab-paclitaxel (TRUCE-01) of bladder cancer (A-J).

Figure 11

Associations between glycolysis-related gene signature and immune-cell infiltration evaluated by seven different approaches. (A) There were significant differences in some immune infiltrate components between low- and high-risk groups. (B, C) Correlation analysis between the riskscore and infiltrating immune cells abundance. *P < 0.05, **P < 0.01 and ***P < 0.001.

Figure 12

Validation of model-associated gene expression. (A) The protein expression levels of the glycolysis-related partial model genes were confirmed by the HPA database. (B) qRT-PCR was utilized to detect the expression of SPINK4, SPINK5, DMRTA1, SLC1A6, and FASN in 10 paired tumor tissues. *P < 0.05.

Figure 13

Loss‐of‐function experiments were conducted to explore the biological function of FASN in vitro. (A) WB assay to examine the FASN protein expression in BCa tissues (T) and normal bladder tissues (N), β-actin protein served as control. (B) FASN small-interfering RNA (siRNA) transfection efficiency was assessed by WB in T24 and UM-UC-3 cells. (C) The wound healing assay demonstrated the capacity of migration in T24 and UM-UC-3 cells. (D) The transwell assay demonstrated the capacity of invasion in T24 and UM-UC-3 cells. (E) The colony-forming assay detected the proliferation ability of tumor cells. (F) CCK-8 assays were utilized to detect cell proliferation ability in T24 and UM-UC-3 cells. ****p < 0.0001.