An Exploration of the Tumor Microenvironment Identified a Novel Five-Gene Model for Predicting Outcomes in Bladder Cancer

Abstract

Bladder cancer (BC) is one of the top ten most common cancer types globally, accounting for approximately 7% of all male malignancies. In the last few decades, cancer research has focused on identifying oncogenes and tumor suppressors. Recent studies have revealed that the interplay between tumor cells and the tumor microenvironment (TME) plays an important role in the initiation and development of cancer. However, the current knowledge regarding its effect on BC is scarce. This study aims to explore how the TME influences the development of BC. We focused on immune and stromal components, which represent the major components of TME. We found that the proportion of immune and stromal components within the TME was associated with the prognosis of BC. Furthermore, based on the scores of immune and stromal components, 811 TME-related differentially expressed genes were identified. Three subclasses with distinct biological features were divided based on these TME-genes. Finally, five prognostic genes were identified and used to develop a prognostic prediction model for BC patients based on TME-related genes. Additionally, we validated the prognostic value of the five-gene model using three independent cohorts. By further analyzing features based on the five-gene signature, higher CD8+ T cells, higher tumor mutational burden, and higher chemosensitivity were found in the low-risk group, which presented a better prognosis. In conclusion, our exploration comprehensively analyzed the TME and identified TME-related prognostic genes for BC, providing new insights into potential therapeutic targets.

Keywords: bladder cancer; immune infiltration; prognosis; stroma; tumor microenvironment.

Figure 1

Correlation of scores with the survival of BC patients and clinicopathological staging characteristics. (A–C) Kaplan–Meier survival analysis for BC patients of high and low score in ImmuneScore, StromalScore and ESTIMATEScore. p = 0.57, 0.014 and 0.045 by Log Rank test. (D) Distribution of ImmuneScore in stage, T classification, N classification and M classification. The p = 0.14, 0.16, 0.99 and 0.044, respectively. (E) Distribution of StromalScore in stage, T classification, N classification and M classification. The p = 4.9e−09, 3.3e−07, 0.012, 0.55, respectively. (F) Distribution of ESTIMATEScore in stage, T classification, N classification and M classification. The p = 8.4e−05, 23e−04, 0.26 and 0.56, respectively.

Figure 2

Heatmaps, Venn plots, and enrichment analysis of GO and KEGG for DEGs. (A, B) Heatmap for DEGs generated by comparison of the high score group vs. the low score group in StromalScore and ImmuneScore. Row is the gene, and column name is the samples which not shown in plot. Differentially expressed genes were determined by Wilcoxon rank sum test with q <0.05 and log2foldchange >1.5 as the significance threshold. (C, D) Venn plots showing common up-regulated and down-regulated DEGs shared by StromalScore and ImmuneScore. (E, F) GO and KEGG enrichment analysis for 811 DEGs, terms with p and q <0.05 were considered significant.

Figure 3

Classification of BC patients. (A) Heatmap for the expression of top 200 high variable genes. Clinical information was annotated for subclasses. (B) PCA plot for the assigned subclasses. (C) Kaplan–Meier survival analysis for BC patients grouped into three subclasses. Log Rank test comparing overall survival between the class1 and class3 subtypes reached a p <0.001. For the comparison between class 2 and class 3, p value <0.01. For class 1 and class 2, p = 0.067.

Figure 4

Immune and stromal features of subclasses. (A) Immune and stromal scores among three subclasses. (B) Immune-related molecular biomarkers among the three subclasses. (C) Stromal-related molecular biomarkers among three subclasses. (D) The fractions of 22 immune cells calculated by CIBERSORT among three subclasses.

Figure 5

The biomarkers of other subtypes and mutations among three subclasses. (A) The biomarkers of other published subtypes among three subclasses. (B) Percent of Robertson’s subtypes. (C) Top mutations in class 1, class 2 and class 3.

Figure 6

Characteristics of the prognostic gene signature. The prognostic index was imputed as follows: (0.234 ∗ FPR1) + (0.051 ∗ TNFAIP6) + (0.145 ∗ GFPT2) + (0.106 ∗ IL-10) + (−0.172 ∗ ZNF683). (A) Identification of the optimal penalization coefficient lambda in the Lasso regression model. (B) LASSO Cox regression algorithm was used to identify the most robust prognostic genes. (C) Forest plots presenting the multivariate Cox proportional hazards regression analysis of prognostic selected genes in overall survival. (D) Kaplan–Meier curves for patients grouped by expression levels of selected genes.

Figure 7

Time-dependent ROC analysis and Kaplan–Meier analysis for the validation of prognostic model in TCGA, GSE13507, GSE31684 and GSE32894. Kaplan–Meier curve and time-dependent ROC analysis of risk score in (A) entire TCGA cohort, (B) GSE13507 cohort, (C) GSE31684 cohort, (D) GSE32894 cohort.

Figure 8

The relationship of five-gene risk groups with functional annotations, TME landscape, immune cell fractions, tumor mutation burden, PDL1 expression and chemotherapy response. (A) The Hallmark enrichment of high- and low-risk groups by GSEA method. (B) Landscape of immune and stromal microenvironment based on immune and stromal signature. (C) Boxplot of the distribution of 22 immune cells in the high- and low-risk groups. (D) PDL1 expression differences. (E) Tumor mutation burden difference. (F) Estimated IC50 indicates the efficiency of chemotherapy to high- and low-risk groups by Methotrexate.

Figure 9

CMap database analysis identifies potential cancidate small molecular drugs targeting the DEGs between high- and low-groups.