Integrated analysis of bulk and single-cell RNA-seq data reveals cell differentiation-related subtypes and a scoring system in bladder cancer

Abstract

Bladder cancer (BLCA) exhibits notable molecular heterogeneity, influencing diverse clinical outcomes. However, the molecular subtypes associated with cell differentiation-related genes (CDR) and their prognostic implications remain unexplored. Analysing two GEO single-cell datasets, we identified genes linked to cell differentiation. Utilizing these genes, we explored distinct molecular subtypes. WGCNA analysis further identified CDR-associated genes, and the CDR score system, constructed using Lasso and Cox regression, was developed. Clinical prognosis and variations in immune-related factors among patient groups were assessed. Core genes were selected and confirmed through in vitro experiments. Two BLCA subtypes related to cell differentiation were identified: Subtype B demonstrated a favourable prognosis, while Subtype A exhibited significant immune cell infiltration. The CDR score system of nine genes revealed a positive correlation between higher scores and a poorer prognosis. The comprehensive analysis uncovered a positive association between CDR genes and M2 macrophages and unresponsiveness to immune therapy. Functional experiments validated that ANXA5 downregulation influences tumour cell migration without affecting proliferation. Our study reveals distinct cell differentiation-related molecular subtypes and introduces the CDR scoring system in BLCA. ANXA5 emerges as a potential therapeutic target, offering promising avenues for personalized treatment strategies.

Keywords: WGCNA; bladder cancer; cell differentiation; single‐cell data; subtypes.

FIGURE 1

The entire process of this study.

FIGURE 2

Analysing single‐cell RNA sequencing datasets derived from GSE190888 and GSE135337. Clusters annotation and cell types identification via tSNE in GSE190888 (A, B) and GSE135337 (D, E). Cell trajectory analysis for the identified epithelial cells in GSE190888 and GSE135337 (C, F). (G)Venn diagrams illustrate the overlap of CDR‐related genes identified in the two datasets, GSE190888 and GSE135337. (H). Gene Ontology (GO) analysis of CDR‐related genes was conducted to investigate the biological processes, molecular functions and cellular components associated with these genes.

FIGURE 3

Cell Differentiation‐related (CDR) Molecular Subtype Identification and WGCNA analysis in bladder cancer. (A) Consensus cluster analysis indicated that the optimal number of clusters was k = 2, dividing patients into two distinct subgroups. (B) Principal Component Analysis (PCA) analysis of two distinct subtypes. (C) Kaplan–Meier survival analysis between subtype A and B. (D) Analysis of immune cell infiltration in the two CDR subtypes was performed using ssGSEA (single‐sample Gene Set Enrichment Analysis). Statistical significance was indicated as follows: *p < 0.05; **p < 0.01; ***p < 0.001. (E, F) The heatmap demonstrates the relationship between various gene modules and clinical characteristics in GSE13507 and TCGA.

FIGURE 4

Clinical analysis of CDR scoring system. RNA sequencing data from 594 bladder cancer patients (TCGA and GSE13507) were combined, split into training and validation sets (1:1 ratio), and patients were categorized into high and low‐scoring groups based on the median score. K‐M survival analysis of CDR score in entire cohort (A), train cohort (B) and test cohort (C). ROC curves depicting the prediction of overall survival at 1, 3 and 5 years in the complete cohort (D), training subset (E) and test subset (F). Correlation between CDR subtypes and CDR scoring systems (G and H). Relationship between the scoring system and clinical pathological features (I‐L).

FIGURE 5

Correlation Between the Scoring System and the Immune Microenvironment. Three violin plots were used to visualize the expression of 9 genes related to the nine scoring system genes across the datasets GSE190888, GSE135337 and the external validation dataset GSE145281 (A, B, C). Correlation analysis between the scoring system genes and immune cells is shown in (D). Statistical significance was indicated as follows: *p < 0.05; **p < 0.01; ***p < 0.001. Analysis of immunotherapy sensitivity in both high and low‐scoring groups within the TCIA dataset is presented in (E). isp_ctla4_pos_pd1_neg refers to patients with positive CTLA4 expression and negative PD1 expression. isp_ctla4_neg_pd1_pos denotes patients with negative CTLA4 expression and positive PD1 expression. isp_ctla4_neg_pd1_neg represents patients with both CTLA4 and PD1 expression negative. isp_ctla4_pos_pd1_pos indicates patients with both CTLA4 and PD1 expression positive. Higher scores on the vertical axis represent greater sensitivity to immunotherapy.

FIGURE 6

Functional Verification of Knockdown ANXA5. ANXA5 expression in bladder cancer tissues and cells was detected using qPCR, and the interference efficiency of small interfering RNA (siRNA) was also validated (A, B and C). Evaluation of cell migration ability using wound healing and transwell assays was conducted to assess the effect of ANXA5 knockdown on bladder cancer cell migration (D and E). Statistical significance was indicated as follows: *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.