Your new experience awaits. Try the new design now and help us make it even better

ORIGINAL RESEARCH article

Front. Immunol., 17 October 2025

Sec. Cancer Immunity and Immunotherapy

Volume 16 - 2025 | https://doi.org/10.3389/fimmu.2025.1600254

This article is part of the Research TopicImmune Predictive and Prognostic Biomarkers in Immuno-Oncology: Refining the Immunological Landscape of CancerView all 46 articles

Comprehensive analysis of single-cell and bulk transcriptomes reveals key B-cell genes and immune microenvironment regulation in bladder cancer

Lijun Wang&#x;Lijun Wang1†Juan Yang,&#x;Juan Yang1,2†Zhangxiao Xu,&#x;Zhangxiao Xu1,2†Bo TaoBo Tao1Yunpeng HeYunpeng He1Yuan ZhaoYuan Zhao1Jian WuJian Wu1Yiran MaYiran Ma1Zitao Zhong*Zitao Zhong2*Lin Ye*Lin Ye1*
  • 1Faculty of Life Science and Technology & The Affiliated Anning First People’s Hospital, Kunming University of Science and Technology, Kunming, China
  • 2Faculty of Life Science and Technology, Kunming University of Science and Technology, Kunming, Yunnan, China

Background: Bladder cancer is a significant malignancy, for which prognostic prediction and understanding of the tumor immune microenvironment are crucial. B cells play a key regulatory role in this environment, making their study essential for advancing bladder cancer research.

Methods: In this study, a multi-omics analysis strategy combining single-cell RNA sequencing (scRNA-seq) and bulk RNA-seq was used to establish single-cell transcriptome profiles of tumor tissues from bladder cancer patients, focusing on B-cell populations and their interactions with other cell types in the tumor microenvironment. Large public databases were used to screen for key prognostic genes associated with bladder cancer B cells, and their biomarker expression was verified by in vitro experiments.

Results: Based on tumor samples from eight patients with bladder cancer and four normal samples, we selected 84, 967 cells for single-cell sequencing analysis. From these, we identified 10, 967 B cells and identified 508 key genes associated with B cells in bladder cancer from five different B cell subtypes. By integrating a large amount of RNA sequencing data, we identified VCL, FLNA, TAGLN, ACTA2, COL6A2, and CALD1 as potential biomarkers for B-cell-associated bladder cancer, and experimentally verified that these markers were significantly lower in bladder cancer patients than in the normal group, and were effective in predicting the survival rate of the patients and the status of the tumor immune microenvironment.

Conclusions: Using a combination of transcriptomic and experimental validation at single-cell and batch levels, this study provides insights into the key gene signatures of B cells from patients with bladder cancer and their roles in regulating the tumor immune microenvironment, providing new biomarkers and potential therapeutic targets for predicting patient’ prognosis and immunotherapy response.

1 Introduction

Bladder cancer (BCa) is a malignant tumor originating from the epithelium of the urinary tract, which continues to increase in morbidity and mortality, and is one of the most common malignant tumors (1). Epidemiological estimates for 2025 indicate that bladder cancer is the fourth most common cancer among men in the United States (n = 65,080), with a total of 84,870 new cases nationwide. The incidence in men is 3.3 times higher than that in women, showing a significant difference (2). Bladder cancer is a highly heterogeneous disease, with many challenges in its classification, staging, and grading (3, 4). Bladder cancer is characterized by high rates of recurrence and metastasis. Based on the depth of tumor infiltration, they can be divided into non-muscle invasive bladder cancer (NMIBC) and muscle invasive bladder cancer (MIBC) (5). Of these, 75% of uroepithelial carcinomas of the bladder present as non-muscle invasive bladder cancer (6), and the other 25% of urothelial carcinomas of the bladder present as muscle-invasive bladder cancer with a risk of metastatic spread, with a 5-year survival rate of approximately 40%-50% (7), approximately 50% of bladder cancer cases with muscle invasive bladder cancer eventually metastasize, and the 5-year survival rate for muscle invasive bladder cancer with distant metastases is only 10% (8), which is the main cause of death in bladder cancer patients. Currently, radical cystectomy (RC) combined with chemotherapy or immunotherapy is the first-line treatment for BCa (9). Although chemotherapy and immunotherapy can improve survival to some extent, a subset of patients respond poorly to these therapies, resulting in missed opportunities for RC and reduced survival (10, 11). However, in some cases, it may not be effective in preventing cancer recurrence and metastasis (12). Therefore, recurrence, metastasis, and spread of cancer have become one of the biggest resistances to cancer treatment, and it is extremely important to study the mechanism of bladder cancer development, metastasis, and spread, and to find new targets for bladder cancer treatment.

The tumor microenvironment (TME) has been a hotspot in cancer biology research and is a relevant therapeutic target for drug discovery. Notably bladder cancer is one of the cancers with the least immune infiltration (13), This may account for the poor response to anti-PD1 therapy. Molecular subtypes of bladder cancer show different cell type-specific expression patterns (14), suggesting that the heterogeneity of BCa is, at least in part, due to the different cell type components of the microenvironment. However, until recently, relevant studies have been scarce. Previous studies have shown that the abundance of B cells in cancer is positively correlated with a favorable clinical outcome, whereas others have shown that they promote tumors, implying that the biological function of B cells is a complex landscape (15). Relatively little is known about comprehensive molecular analyses of B cells in bladder cancer. Therefore, the biological functions of B cells in bladder cancer remain to be explored.

The rise in single-cell RNA sequencing (scRNA-seq) technology has provided an unprecedented opportunity to resolve the molecular signatures of different immune cell populations in the TME (16). The scRNA-seq provides more precise and detailed analyses at the single-cell level than traditional bulk RNA-sequencing methods (17). In this study, we innovatively combined bulk RNA seq and scRNA-seq data to systematically reveal the infiltration pattern of B cells in the microenvironment of bladder cancer, and successfully identified seven bladder cancer markers with potential clinical applications. In this study, the molecular features of tumor-infiltrating B cells were explored in depth and their specific marker genes were identified through single-cell RNA sequencing analysis of bladder cancer samples. RNA sequencing data from large databases were integrated and analyzed to screen and validate potential key genes associated with cancer immune responses. The workflow of this study is illustrated in Figure 1.

Figure 1
Flowchart illustrating analysis workflows. At the top, “Transcriptome Analysis” involves GSE13507 to WGCNA, and TCGA to Difference Analysis, leading to Intersection Analysis. On the right, “Construction of the protein-protein interaction (PPI) network” follows, combined with Survival Analysis. Below, “Single-cell Analysis” starts with GSE22231, progressing through data integration, dimensionality reduction, data clustering, cell annotation, DEGs, marker gene identification, and enrichment analysis. Analysis of diagnostic significance, ROC analysis, TCGA, GSEGSE37815, and Wilcoxon rank-sum test-RT-qPCR is shown. Relevant genes listed include VCL, FLNA, TAGLN, ACTA2, COL6A2, CALD1.

Figure 1. Research design.

2 Methods

2.1 Data source and preprocessing

In this study, a multi-omics integration analysis strategy was used to integrate four independent datasets from The Cancer Genome Atlas (TCGA) (https://www.cancer.gov/ccg/research/genome-sequencing/tcga) bladder cancer dataset, and the Comprehensive Gene Expression Database (GEO) (https://www.ncbi.nlm.nih.gov/) (Table 1). Eight bladder cancer tissue samples and four normal tissue samples from nine patients in the GSE222315 dataset were selected for the scRNA-seq study. The selection criteria were as follows:(a) each sample should contain no less than 300 and no more than 7000 cells; (b) each cell should express more than 300 (250) genes; (c) each gene should be expressed in at least three cells; and (d) the mitochondrial RNA content in each cell should be less than 20% (18). Finally, we screened using the R package Seurat and obtained 84, 967 cells for subsequent analysis. Differentially expressed genes (DEGs) were identified from TCGA bladder cancer dataset using the DESeq, limma, and Wilcoxon function packages of the R package. Specifically, we applied a threshold of |log2 FoldChange| > 2 with an adjusted P < 0.05 to select DEGs for subsequent analysis (19). The GSE37815 dataset was used for marker gene validation.

Table 1
www.frontiersin.org

Table 1. Grouping information.

2.2 Data integration

After data integration and filtering, the Seurat package was used for data normalization and integration analyses. First, the gene expression values were normalized using the ‘NormalizeData’ function. Specifically, the expression value of each gene was multiplied by the total gene expression of the cell, and the scaling factor (10,000) and the natural logarithmic transformation (log(x+1)) were performed to avoid a zero logarithm. Next, 2000–3000 highly variable genes were identified by the ‘FindVariableFeatures’ function and centered using the ‘ScaleData’ function. To eliminate technical differences between batches of samples (due to the inability to sequence samples simultaneously), we used the ‘FindIntegrationAnchors’ or ‘RunHarmony’ function to set up 2000 anchors for data integration. This anchor-based integration strategy can effectively map homologous cell types from different datasets to a small number of key anchors, thereby reducing the batch effect and improving the comparability and reliability of data.

2.3 Dimensionality reduction and data clustering

Dimensionality Reduction and Data Clustering: Considering that each gene in the sample exists as a dimension, high-dimensional data becomes difficult to visualize. Therefore, dimensionality reduction techniques must be employed to represent the true data structure by using a reduced number of dimensions (20). First, the highly variable set of genes obtained was downscaled using the ‘RunPCA’ function in the Seurat package. Subsequently, data integration was performed by setting 2000 anchors using the ‘FindIntegrationAnchors’ (or ‘RunHarmony’) function. This ‘anchor’ based integration strategy effectively converts the same cell types from different datasets into a small number of key anchors, thus mitigating batch effects. Next, the ‘FindNeighbors’ and ‘FindClusters’ functions in the Seurat package were used to perform cluster analysis of the downscaled data. To enhance the visual presentation of the clustering results, the t-distributed stochastic nearest neighbor embedding (t-SNE) method was applied. Specifically, a clear visualization of the cell population was achieved by the ‘RunTSNE’ function of the Seurat package with the parameters dims=1:15 and resolution=1.0.

2.4 Cell annotation, DEGs, and marker genes identification

Cellular annotation, DEG and marker gene identification were performed by reviewing relevant literature and websites and manually searching for cellular marker genes for cell-type annotation. This approach effectively correlates gene expression in different cell types with that in cells. The weakest correlation for each cell type was eliminated through iterations and the corresponding cell types were identified (21). In this study, samples were manually annotated for cell types using literature and web resources. Subsequently, DEG were identified using the ‘FindAllMarkers’ function in the Seurat software package. A threshold value of |log2FoldChange| > 2 (0.25 between cells and between groups) was set and adjusted to P < 0.05. Finally, marker genes specific to each cell type were identified using the ‘FindAllMarkers’ function of the Seurat software package. This series of analyses provides an important basis for precise annotation and functional study of cell types.

2.5 GO enrichment analysis and KEGG pathway enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the David and Metascape databases (clusterProfler packages, version 3.14.3). The analysis revealed functional enrichment of genes and related signalling pathways. All the results were visualized using the ggplot2 software package (22, 23).

2.6 Cell trajectory analysis

Cells constantly undergo dynamic changes, transforming from one cell type to another, which leads to changes in gene expression and functional state (24). Each cell was arranged along a corresponding cell track, representing a pseudotemporal order, and cells were grouped into different differentiated states by employing gene expression profiles. We used the Monocle2 package for the pseudotemporal analysis of B-cell subtypes (25). To explore the differentiation trajectories of B cell subtypes and related genes in different states, we used the ‘plot_cell_trajectory’ function to sort the cells according to their pseudotimes. The ‘BEAM’ function was used to identify genes responsible for cell branching and differentiation. The ‘plot_genes_branched_heatmap’ function was used to visualize the results.

2.7 Cell-cell communication analysis

The analysis of cell-cell communication is based on the interaction of ligands and receptors on the cell surface, and this intercellular communication plays a key role in a variety of biological processes. To investigate the interaction patterns between different cell types, the CellChat software package was used to construct cellular communication networks. The software simulates the communication between cells by constructing an interaction network of ligands, receptors, and their related factors. Based on the gene expression profiles of ligands and receptors in different cell types, CellChat was able to infer the strength of cellular interactions and reveal the rich patterns of ligand-receptor interactions between different cell types. This analysis provides an important basis for a deeper understanding of intercellular signalling and functional regulation (26).

2.8 Construction of weighted gene co-expression network

Weighted gene co-expression networks were constructed using the R package ‘WGCNA’ (27). Sample analyses were initiated by screening for outlier samples using clustering methods. Subsequently, the gene expression matrix was converted into a similarity matrix by calculating the Pearson correlation value (cmn) and adjacency (amn) between genes, where parameter β served as a soft threshold that could modulate the correlation between genes. Based on the assessment of scale independence and average connectivity, the soft threshold power was set to 8 (β = 8, R2 = 0.9). The neighbor-joining matrix was further transformed into a topological overlap matrix so that the gene interactions conformed maximally to the scale-free distribution characteristics. Hierarchical clustering was performed using the dynamic tree-cutting algorithm with the unsigned network construction type, which was set to two to determine the module division, and the minimum module size of the gene dendrograms was set to 30. The module difference degree (MES) of the module dendrograms was calculated, and modules with a difference degree of less than 0.3 were merged to construct the final network structure.

2.9 Differentially expressed genes for bladder cancer-related genes

412 bladder cancer gene expression data and 19 paraneoplastic tissue gene datasets were downloaded from TCGA database. The downloaded gene datasets were integrated into a matrix using Perl (command line) software. Screening criteria for identifying DEGs were established using the R packages DESeq, limma, and Wilcoxon. Specifically, we applied a threshold of |log2 fold change | > 2 with an adjusted P < 0.05. The DEGs were selected for subsequent analyses. A Venn diagram was used to analyze the intersection of three sets of gene sets: sc-seq tagged genes from bladder cancer-associated B cells, high ESTIMATEScore-related genes obtained from WGCNA analysis, and differentially expressed genes screened by the TCGA database using the three methods to identify core overlapping genes.

2.10 Integrated survival, ROC, and PPI network analyses for the identification of bladder cancer biomarkers

To select the core genes with prognostic value, Kaplan-Meier survival analysis was performed on these key genes using the R software package ‘survivalPath’ (28), and survival differences were evaluated using the log-rank test. Genes with P < 0.05 were considered to have prognostic significance and were subjected to further in-depth analysis. Assessment of the diagnostic significance of hub genes. To investigate whether prognostic genes could distinguish tumor samples from non-tumor samples, we used the R package ‘survival ROC’ (29) and ROC analyses were performed on the hub genes. To elucidate the interaction network of infiltrating immune cell-associated genes at the protein level, protein-protein interaction (PPI) network mapping was constructed using the STRING database (https://string-db.org/). The interactions were visualized using the Cytoscape 3.7.1 platform to clearly demonstrate the functional association patterns among the genes.

2.11 Core gene database and protein validation

The expression levels of the core genes in tumor and non-tumor samples were compared in TCGA and GSE37815 datasets. Additionally, the protein expression of these genes in bladder cancer was validated using the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/).

2.12 Validation of core genes in cells and clinical samples by RT-qPCR

A normal human uroepithelial cell line (SV-HUC-1) and three human bladder cancer cell lines (5637, T24, and HT1376) were obtained from Yunnan Tengyue Biotech Co. E2112), 5637 cells, and HT1376 cells were cultured in RPMI-1640 medium (Gibco, C11875500BT), T24 cells were cultured in McCoy’s 5A medium (EvaCell, E2110), and the above cells were cultured in medium supplemented with 10% FBS (Gibco, A5256701) and 1% penicillin-streptomycin (Gibco,15140122). Cells were cultured at 37 °C, 5% CO2 cell culture incubator. Total RNA was extracted using the RNAfast200 kit (Shanghai Feijie Biotechnology Co., Ltd.), and mRNA was reverse-transcribed into cDNA using the Evo M-MLV Reverse Transcription Reagent Pre-mix (Hunan Acres Bio, AG11706). RT-qPCR was performed using the SYBR Green Pro Taq Hs Pre-mixed qPCR Kit (Hunan Acres Bio, AG11701). RT-qPCR was performed to verify the expression of the biomarkers in normal human uroepithelial cells (SV-HUC-1) and human bladder cancer cell lines (5637, T24, and HT1376). To further validate the experimental results, we also collected tumor tissues and adjacent normal tissues from seven bladder cancer patients at the First People’s Hospital of Anning, Affiliated Hospital of Kunming University of Science and Technology. All participants provided written informed consent, and the study was approved by the Ethics Committee of the First People’s Hospital of Anning, Affiliated Hospital of Kunming University of Science and Technology. The sequences of all primers are listed in Table 2.

Table 2
www.frontiersin.org

Table 2. Primers of hub genes.

2.13 Data processing

Statistical analysis was performed using the Statistical Analysis R software (version 4.4.3) and GraphPad Prism 9.0 (GraphPad Software, USA). Differences between samples were assessed using the Wilcoxon rank-sum test, and statistical significance was determined using Student’s t-test. Experimental results are expressed as the mean ± standard deviation (mean ± SD), with all data derived from at least three independent experiments. Significance was denoted as follows: NS: p > 0.05, *: p < 0.05, **: p < 0.01, ***: p < 0.001, and ****: p < 0.0001.

3 Results

3.1 Single-cell transcriptome profiles of bladder cancer

Twelve samples (including tumors and paired normal tissues) from nine patients in the GSE222315 dataset were analyzed in this study. After rigorous quality control screening, 84,967 high-quality cells were obtained for subsequent analysis. After eliminating the batch effect using the ‘anchor’ integration strategy, the data were normalized, centered and downscaled by principal component analysis (PCA) to retain the top 30 principal components. Subsequent clustering visualization using the t-SNE algorithm classified the cell population into 23 cell clusters (Figure 2A), which were grouped into seven major cell types by manual annotation (Figure 2B). Figure 2C shows the single-cell transcriptome profiles of the normal and tumor samples. Figure 2D shows single-cell transcriptome profiles of different samples using manual cellular annotation revealing seven cell types (Figure 2E): CD3D, CD3E, CD2, GNLY, KLRD1 high-expressing T cells, EPCAM high-expressing epithelial cells, CD79A, MZB1, MS4A1 high-expressing B cells, FCGR3A, TYROBP-overexpressing NK cells, COL3A1, COL1A1, DCN, C1R-overexpressing Fibroblast cells, PECAM1, CD34, CDH5, VWF-overexpressing Endothelial cells, and KIT, TPSAB1, and TPSB2-overexpressing Mast cells. Differential gene expression heatmaps (Figure 2F) and UMAP maps (Figure 2G) were used to visualize the transcriptional features and marker gene expression patterns of the seven cell types. In addition, we analyzed the proportional distribution of each cell type in the normal and tumor tissues (Figure 2H).

Figure 2
Seven-panel scientific figure. Panel A: UMAP plot with 22 clusters differentiated by color. Panel B: UMAP plot labeling seven cell types. Panel C: UMAP plot showing normal versus tumor groupings. Panel D: Series of UMAP plots displaying cell distributions across samples p1 to p9. Panel E: Dot plot showing gene features by cell identity, indicating percent and average expression. Panel F: Heatmap showing gene expression levels across different cell identities. Panel G: Seven individual UMAP plots for specific genes with expression intensity. Panel H: Stacked bar chart showing the proportion of cell clusters in normal and tumor samples.

Figure 2. Overview of single-cell mapping of normal and tumor samples from bladder cancer cancer patients. (A) Umap plot depicting single cell samples clustered into 23 clusters. (B) Identification of 7 cell types based on marker gene expression. (C) Umap clustering diagram comparing normal and tumor tissues. (D) Umap clustering diagram of 12 samples. (E) Bubble plots showing marker gene expression of 7 cell types across cells. (F) Heatmap showing differential gene expression across 7 cell types. (G) Umap plot highlighting expression patterns of marker genes across 7 cell types. (H) Cell type sections of normal and tumor samples in the scRNA-seq dataset.

3.2 Single-cell transcriptome mapping of B cells in scRNA-seq

ScRNA-seq analysis of bladder cancer revealed a high abundance of B cells in all samples. In total, 10,967 B cells were extracted for clustering and annotation. Secondary clustering divided these cells into four clusters (Figure 3A), which were manually annotated based on marker gene expression patterns (Figure 3B). A comparison of UMAP clustering between normal and tumor tissues revealed significant differences (Figure 3C). Figure 3D illustrates the transcriptomic features of B cells across samples, which were classified into three identifiable B cell subtypes and one unannotated subtype: plasma cells with high CD38 expression, mature B cells with high CD19/CD22 expression, and memory B cells with high CD27 expression (Figure 3E). Further analysis demonstrated that these five B cell subtypes exhibited significant differences in gene expression between normal and tumor tissues (Figure 3F), and the UMAP plots visually depicted the expression of marker genes across the subtypes (Figure 3G). Overall, this single-cell transcriptomic study provides a comprehensive characterization of abundant B cell populations in the bladder tumor microenvironment, offering insights into their potential roles in tumor immune regulation.

Figure 3
UMAP plots and heatmaps visualize cell clustering and gene expression in normal and tumor samples. Panels A-C show B cells, plasma cells, and macrophages, with color-coded identities. Panel D illustrates sample distribution across UMAP dimensions. Panel E displays feature expression levels like CD38 and CD19 across cell types. Panel F is a heatmap showing gene expression with identities. Panel G presents expression of genes like CD83 and SDF2L1. Panel H compares cluster proportions in normal versus tumor states.

Figure 3. Single-cell transcriptome mapping of B cells in scRNA-seq. (A) Umap plot of B cells in scRNA-seq, clustered into 5 different clusters. (B) Sequence-based identification of B based on 5 cell types. (C) Umap clustering plot comparing normal and tumor tissues. (D) Umap clustering plot of 12 samples. (E) Bubble map showing marker genes for 5 B-cell subtypes. (F) Heatmap showing differential gene expression for 5 cell types. (G) Umap plot highlighting the expression patterns of marker genes for the 5 cell types. (H) River bar stack plot of cell types for normal and tumor samples in the scRNA-seq dataset.

3.3 Enrichment analyses associated with B cells

To investigate specific gene expression among the different cell subtypes, we first performed a volcano plot analysis of the seven cell subtypes (Figure 4A). To elucidate the biological functions of the marker genes in B cells, KEGG and GO enrichment analyses were performed. KEGG results showed that marker genes were mainly enriched in pathways activated to function in B cells, including ribosome, coronavirus disease-COVID-19, intestinal immune network for lgA production, B cell receptor signalling pathway, hematopoietic cell lineage, Leishmaniasis, Antigen processing and presentation, Toxoplasmosis, Allograft rejection, and Type I diabetes mellitus (Figure 4B). GO enrichment results showed that biological processes (BP) were mainly associated with cytoplasmic translation, immunoglobulin production, production of molecular mediators of immune response, B cell-mediated immunity, and immunoglobulin-mediated immune response. Cellular Component (CC) is mainly related to immunoglobulin complexes, cytosolic ribosomes, ribosomal subunits, ribosomes, and cytosolic large ribosomal subunits. Molecular function (MF) was associated with antigen-binding and structural constituents of ribosomes (Figures 4C, D). We also mapped the differential gene scatter plots (Figure 4E) and volcano plots (Figure 4F) of B cells between tumor and normal tissues and performed differential gene enrichment analyses (Figure 4G) to further explore the changing characteristics of B cells in the tumor microenvironment.

Figure 4
Multiple graphs display data analysis results. Graph A shows cell cluster changes with average log2 fold change on the vertical axis. Graphs B to D illustrate pathway enrichment, highlighting KEGG and GO pathways with varying p-values and counts. Graph E is a volcano plot showing gene expression data, and Graph F compares normal versus tumor expression levels, highlighting specific genes like RFX7 and IGHG. Graph G presents a bar chart of enriched pathways related to extracellular matrix and ribosomal components.

Figure 4. Enrichment analyses between cells of specific subpopulations and between groups. (A) Overall cell subpopulation intercellular differential analysis volcano plot (red dots are genes specific to cell clusters relative to other cell clusters). (B) B-cell intercellular differential gene KEGG functional enrichment map. (C) Functional enrichment of the B cell intercellular differential gene GO. (D) B-cell intercellular differential gene GO functional enrichment map (faceted to show the enrichment of BP, MF, and CC in GO entries). (E) Scatter plot of B-cell intergroup difference analysis. (F) Volcano plot of B-cell intergroup difference analysis. (G) Bar graph of differential gene enrichment between B-cell groups.

3.4 pseudotemporal analysis

Pseudotemporal analysis, also known as cell trajectory analysis, simulates developmental trajectories of different cells based on the expression patterns of temporal genes in single-cell samples. B-cell maturation trajectories were analyzed in real time using the Monocle software package for tumor-associated B cells. We extracted B cells to demonstrate their developmental trajectories, revealing five branches of B cell subtype (Figure 5A). As shown in Figures 5B, C, B cells exhibited five differentiation states during development. Figure 5B depicts the chronological order of cell subtype differentiation, with darker colored cells gradually transforming into lighter colored ones. The results showed progression from Branch 3 to Branch 5, Branch 1, and Branch 2. There were differences in the expression of model genes at different developmental stages of eutrophication. During the development of tumor-associated B cells, the expression of DERL3, FKBP11, FKBP2, ITM2C, MZB1, PRDX4, RRBP1, and XBP1 gradually increases. Subsequently, the temporal gene expression of the four branches was presented as a heatmap using the BEAM function. BPs were explored by performing GO enrichment analysis, which showed that branch 1 was mainly associated with the immune response-activating signalling pathway, B cell activation, B cell receptor signalling pathway, B cell proliferation, B cell differentiation, and B cell-mediated immunity. Branch 2 is mainly involved in B-cell activation, homeostasis, proliferation, and differentiation. Branch 3 is mainly related to the cellular response to tumor necrosis factor, the response to tumor necrosis factor, and positive regulation of the humoral immune response. Finally, branch 4 was mainly related to the cellular response to tumor necrosis factor, response to tumor necrosis factor, and positive regulation of the humoral immune response (Figure 5E). These pathways are mainly involved in B cell differentiation and immune responses.

Figure 5
Charts illustrate immune cell data and gene expression analysis. Panels A-C display component plots differentiating cell states or types and pseudotime. Panel D highlights clusters linked with immune responses, ribosome biogenesis, and more, with a heat map showing expression levels. Panel E presents relative expression of specific genes over pseudotime for different cell stimuli, color-coded by cell type or cluster state.

Figure 5. Pseudotime analysis of specific subpopulations. (A) Position of different subpopulations in the pseudotime (stim). (B) Differentiation of different subpopulations in the pseudotime (Pseudotime). (C) Position of cells in different pseudotime stages in the pseudotime (State). (D) Heatmap of distribution of genes with differences in the pseudotime as well as the enrichment situation (Heatmap of genes with changes in differences in the pseudotime is shown in four clusters; text in the left bar is the enrichment pathway). (E) Scatterplot of gene changes in the pseudotime.

3.5 Analysis of cell-cell interactions associated with B cells

To explore the communication properties of B cells with other cell types, the Cellchat function was used to detect ligand-receptor and molecular interactions between different cells. The results showed that the interactions between the seven cell types (T cells, B cells, NK cells, fibroblast cells, endothelial cells, and mast cells) were stronger in bladder cancer tissues than in normal tissues (Figures 6A, B). The FN1 signalling pathway was upregulated in tumor samples and was widely present in ligand-receptor interactions in various cell types (Figures 6C-E). Among these, the ligand-receptor pair FN1-CD44 contributed the most to the FN1 pathway (Figure 6F), and SCD1 and ITGB1 showed higher expression levels in multiple cell types (Figures 6G, H).

Figure 6
Diagram comparing cellular interactions and signaling pathways between normal and tumor tissues. Panel A shows interaction networks; B presents bar graphs of interaction counts and strengths; C highlights relative information flow; D is a bubble plot of communication probabilities; E displays heatmaps of FN1 signaling; F shows a bar chart of ligand-receptor pair contributions; G illustrates chord diagrams of cell interactions; H contains violin plots of relative contributions by cell type.

Figure 6. Cell–cell communication analysis. (A) Circos plots of intercellular interactions in normal and tumor groups. (B) Comparison of the number and strength of cell–cell communications between tumor and normal tissues. (C) Differential signaling pathway-related communications between tumor and normal tissues. (D) Communication probabilities of macrophages with other cell types mediated by ligand–receptor pairs. (E) Heatmap of FN1 signaling interactions among different cell types in tumor and normal groups. (F) Relative contribution of ligand–receptor pairs within the FN1 pathway in tumors. (G) Cell–cell communication mediated by the FN1CD44 ligand–receptor pair. (H) Expressiondifferences of FN1 pathway molecules between tumor and normal samples.

3.6 Identification of the hub gene in bladder cancer

WGCNA-based identification of infiltrative immune cell-related genes to screen infiltrative immune cell-related genes, we performed WGCNA assays on the GSE13507 dataset. First, sample clustering analysis showed that all 197 samples in the GSE13507 dataset could be used to construct weighted gene co-expression networks (Figure 7A). Subsequently, a weighted gene co-expression network containing eight modules was constructed when β was set to 5 (R2 = 0.85), and modules with module feature gene dissimilarity less than 0.3 were merged (Figures 7B, C). The associations between these eight modules and the four algorithms (StromalScore, ImmuneScore, ESTIMATEScore, and TumorPurity) are shown in Figure 7D. We found that the MEred module was associated with the ImmuneScore (Figure 7D). In addition, we calculated the correlation between the genes in the MEred module and ImmuneScore and identified 1390 infiltrating immune cell-associated genes by setting the correlation coefficient.

Figure 7
A collage of several data visualizations depicting various analyses in a research study. (A) A sample dendrogram and trait heatmap show hierarchical clustering and score distribution. (B) Two line charts illustrate scale independence and mean connectivity for soft threshold power. (C) A cluster dendrogram with dynamic tree cut and color-coded clusters. (D) A heatmap displays module-trait relationships with values and significance. (E) A Venn diagram shows overlaps among five datasets, highlighting shared elements. (F) A network diagram presents gene interactions with nodes in varying sizes and colors indicating centrality and importance.

Figure 7. Identification of modules associated with bladder cancer using WGCNA. (A) Sample clustering plot; (B) Scale-free index and average connectivity analysis for different soft threshold powers. (C) Dendrogram of all differentially expressed genes clustered according to the heterogeneity measure (1-TOM). Ribbons show results from automated single-block analysis; (D) Heatmap of correlation between module signature genes and bladder cancer immune-related traits. (E) Venn plot showing the results of the intersection analysis of b-cell related differentially expressed genes in single cell samples and bladder cancer TCGA dataset, showing a total of 43 immune-related bladder cancer b-cell key genes. (F) Protein interactions network graph constructed by the identified 43 key genes, in which the red nodes marked the core genes with the top ten BC values. TOM: topological overlap matrix, ME: module eigengene; BC: betweenness centrality.

Gene expression datasets of 412 bladder cancer tissues and 19 paracancerous tissues were obtained from TCGA database. Differential expression analysis between bladder cancer and paracancerous tissues was performed to screen for genes related to bladder carcinogenesis and progression. Three methods, R package DESeq, limma, and Wilcoxon, were used to identify DEGs, and the screening criteria were set as |log2 fold change | > 2 and a corrected P < 0.05. A total of 3104, 2563, and 2836 differentially expressed genes were identified using these three methods, respectively (Figure 7E). This study investigated bladder cancer B-cell marker genes, immune-related genes screened using WGCNA, and bladder cancer-related genes obtained by differential analysis (Figure 7E). A total of 43 shared genes were identified. To explore the interaction network of core genes at the protein level, the PPI network was constructed using the STRING database and visualized using Cytoscape 3.7.1 (Figure 7F) Ten core genes were obtained by screening based on BC values: VCL, DUSP1, FLNA, THBS1, TAGLN, ACTA2, COL6A2, VWF, CXCR4, and CALD1.

3.7 Evaluating the prognostic and diagnostic value of key genes

To assess the association between core genes and bladder cancer prognosis, we systematically evaluated ten candidate genes using Kaplan-Meier survival analysis combined with the log-rank test. The results of the analysis showed that seven genes (VCL, FLNA, THBS1, TAGLN, ACTA2, COL6A2, and CALD1) were significantly correlated (P < 0.05) with the prognosis of bladder cancer patients, and these genes were included in the subsequent in-depth study (Figure 8A). To further validate the diagnostic value of these core genes in distinguishing between rejected and non-rejected samples, we performed subject work characteristic (ROC) curve analyses in the TCGA training and GSE37815 validation sets. The results showed that the area under the curve (AUC) for all seven core genes exceeded 0.7 in both independent datasets (Figure 8B), strongly suggesting that these genes may serve as potential biomarkers for bladder cancer diagnosis.

Figure 8
Graphs displaying survival probability and ROC curves for various markers. Part A shows Kaplan-Meier survival plots for VCL, FLNA, THBS1, TAGLN, ACTA2, COL6A2, and CALD1 with p-values indicating statistical significance. Part B illustrates ROC curves for each marker, with AUC values: VCL 0.771, FLNA 0.809, THBS1 0.786, TAGLN 0.894, ACTA2 0.911, COL6A2 0.808, and CALD1 0.788.

Figure 8. Assessment of the value of hub gene diagnosis. (A) Survival curves of 8 genes with prognostic value (B) Receiver Operating Characteristic (ROC) curves of hub genes in the TCGA dataset.

3.8 Key gene database and protein expression results

We then performed an in-depth analysis of the expression patterns of these key genes in the TCGA and GSE37815 datasets, focusing on comparing the expression differences between tumor tissues and paired non-tumor tissue samples. The results showed that, except for THBS1, the remaining six genes exhibited statistically significant differential expression in both the training and validation sets (P < 0.05); therefore, THBS1 was excluded from subsequent analyses (Figures 9A, B). Ultimately, six genes—VCL, FLNA, TAGLN, ACTA2, COL6A2, and CALD1—were identified as potential marker genes. To investigate protein-level interactions between these six candidate markers and FN1, a protein-protein interaction (PPI) network was constructed. As shown in Supplementary Figure 1, FN1 displayed strong interactions with all the six marker proteins. As shown in Figure 9C, compared to normal tissues, VCL, FLNA, TAGLN, ACTA2, COL6A2, and CALD1 were significantly downregulated in BLCA tissues.

Figure 9
Bar graphs and box plots comparing gene expression levels of VCL, FLNA, THBS1, TAGLN, ACTA2, COL6A2, and CALD1 between control and treated groups with significant differences noted. Tissue section images below show immunohistochemical staining for these genes in normal and bladder cancer (BLCA) tissues, indicating protein localization and expression differences between the tissue types.

Figure 9. (A) TCGA dataset showed that the expression of VCL, FLNA, THBS1, TAGLN, ACTA2, COL6A2, and CALD1 were down-regulated in the tumor samples as compared to normal samples. (B) Analysis of the GSE37815 dataset revealed that VCL, FLNA, TAGLN, ACTA2, COL6A2, and CALD1 was downregulated in tumor samples compared to normal tissues, whereas THBS1 showed no significant difference. The Wilcoxon rank-sum test was used to assess discrepancies between the different samples. Each p-value is written above the box plots (NS: p > 0.05, *: p < 0.05, **: p < 0.01, ***: p < 0.001, and ****: p < 0.0001). (C) Validation of the expression of the six genes in the Human Protein Atlas (HPA) database showed that their expression levels were significantly lower in the tumor samples than in the normal group.

3.9 RT-qPCR results of key genes in cell lines and clinical samples

To validate the expression differences of the six potential biomarkers, further experiments were performed on bladder cancer cells and normal bladder cells. RT-qPCR results showed that compared with normal bladder epithelial cells SV-HUC-1, the expression levels of six biomarkers, namely VCL, FLNA, TAGLN, ACTA2, COL6A2 and CALD1, the expression levels of all six biomarkers were significantly down-regulated (P < 0.05), a result that was highly consistent with the results of our pre-biomarker analysis (Figure 10). Further validation was performed in tumor tissues and adjacent normal tissues from seven bladder cancer patients, and the results showed that the expression changes of the six downregulated genes were consistent with bioinformatics analysis and cell experiments (Figure 11). Overall, Alterations in the expression of VCL, FLNA, TAGLN, ACTA2, COL6A2, and CALD1 may play important roles in the development of bladder cancer, suggesting their potential value in early detection, molecular subtyping, and targeted therapy.

Figure 10
Bar charts comparing mRNA relative expression levels of VCL, FLNA, TAGLN, CALD1, ACTA2, and COL6A2 across different cell lines (SV, 5637, T24, HT1376). Each graph shows a decrease in expression from SV to HT1376. Significant differences are marked with asterisks, indicating statistical significance.

Figure 10. RT-qPCR validation of six key genes (VCL, FLNA, TAGLN, CALD1, ACTA2, and COL6A2) was performed in cell lines. The results showed that these genes were significantly downregulated in bladder cancer cells. Statistical significance is indicated above the bars in the histogram (NS: p > 0.05, *: p < 0.05, **: p < 0.01, ***: p < 0.001, and ****: p < 0.0001).

Figure 11
Graphs comparing mRNA relative expression levels of six genes (VCL, FLNA, TAGLN, CALD1, ACTA2, COL6A2) between paracancerous and cancerous tissues. Each graph shows a significant decrease in expression in cancerous tissues, indicated by an asterisk and horizontal bar over the data points. Data points are scattered with error bars representing variability.

Figure 11. RT-qPCR validation of six key genes (VCL, FLNA, TAGLN, CALD1, ACTA2, and COL6A2) was performed in tumor tissues and adjacent normal tissues from bladder cancer patients. The histogram shows that all six genes were significantly downregulated in tumor tissues. Differences between tumor and adjacent normal tissues were assessed using the Wilcoxon rank-sum test, and statistical significance is indicated above the bars in the histogram. (NS:p > 0.05;:p < 0.05;:p < 0.01;:p < 0.001;****:p < 0.0001).

4 Discussions

Despite advances in preclinical and clinical studies, the prognosis of bladder cancer remains unclear. scRNA-seq technology can reveal cellular features at single-cell resolution and elucidate the regulatory mechanisms of RNA expression in pathological states (30, 31), offering the possibility of accurately identifying prognostic features associated with germinal centres. As important participants in the antitumor immune response, B cells play a key role in the presentation of tumor antigens to T cells (32). In this study, we explored the B cell-associated differentially expressed gene profiles based on scRNA-seq technology and preliminarily explored the potential clinical application value of these genes in early bladder cancer screening and molecularly targeted therapies by integrating TCGA and GEO datasets.

In this study, 52,721 single-cell transcriptomes from bladder cancer and normal samples were analyzed to construct a complete map of the TME of bladder cancer. Seven cell types were identified by annotation: T cells, B cells, endothelial cells, epithelial cells, fibroblasts, NK cells, and mast cells. Cell type identification revealed the presence of 10,967 B cells. Notably, the analysis showed that the expression level of B cells was significantly lower in tumor tissues than in normal tissues. It has been shown that B cells perform differentiated functions in different cancer types (33). B cells secrete immunoglobulins, which in turn inhibit tumor growth (34). Kroeger et al. found clonal B-cell expansion and increased plasmoblastoid infiltration in the tumor tissues of patients with bladder cancer, and B cells, as antigen-presenting cells, activated tumor-specific T-cell responses (35). These studies revealed an immune-mediated relationship between B cells and the pathogenesis of bladder cancer and that downregulation of B cell expression in the tumor microenvironment may lead to a reduction in the production of immune proteins, which in turn exacerbates tumor progression. The association between downregulated B cell expression in bladder cancer and increased bladder cancer risk observed in this study is consistent with the above findings. This study provides new insights into the future development of tumor immunotherapy and highlights the importance of fully considering the role of B cells in tumors.

In addition, GO enrichment analysis showed that B cell marker genes were significantly enriched in five immune-related pathways: immunoglobulin production, production of molecular mediators of the immune response, B cell-mediated immunity, B cell-mediated immunity, immunoglobulin-mediated immune response, and immunoglobulin complex. KEGG results showed that marker genes were significantly enriched in the B-cell receptor signalling pathway, a pathway in which B cells are activated to function. B-cell pseudotemporal trajectory analyses revealed dynamic changes in gene expression over time. These genes were classified into four functional clusters involved in lymphocyte-mediated immunity, protein folding, immunoglobulin, and cytoplasmic translation processes. Intercellular communication network analysis revealed a close association between ligands and receptors in different cells, in which cytokines mainly interact through the FN1 (Fibronectin 1) signalling pathway, a finding that highlights the value of FN1 as a potential biomarker for bladder cancer. Heat map analysis further showed that the expression of genes related to the FN1 signalling pathway was mainly reflected in the interaction between fibroblasts and epithelial cells.

Evidence suggests that FN1 positively regulates the proliferation and migration of a variety of cancer cells and has been identified as a biomarker for a number of cancers, including gastric cancer (36) and cervical cancer (37). Targeted inhibition significantly reduces the proliferation, invasive ability, and metastatic potential of cancer cells. Zhang et al. revealed that FN1 may be involved in the development and progression of bladder cancer and has potential as a prognostic marker and therapeutic target for bladder cancer (38). In this study, we confirmed the critical role of FN1 in bladder cancer through in-depth analyses, identified it as a potential oncogenic pathway in bladder cancer, and elucidated its functional mechanism in fibroblasts. Additionally, our study found that FN1 exhibits strong protein-protein interactions with six proteins, including VCL, FLNA, TAGLN, ACTA2, COL6A2, and CALD1. Previous studies have also reported an interaction between FN1 and VCL, consistent with our findings (39). Based on the above findings, the protein-protein interactions between FN1 and these six identified biomarkers require further experimental validation. These findings not only open a new direction for FN1-targeted therapy but also show important research and clinical application prospects in enhancing the effectiveness of cancer immunotherapy and deepening the knowledge of the dynamic role and molecular mechanisms of FN1 in cancer.

Next, the TCGA bladder cancer dataset was analyzed for differences. Subsequently, immune-related genes were screened using WGCNA in the GEO dataset. Finally, six B-cell marker genes (VCL, FLNA, TAGLN, ACTA2, COL6A2, and CALD1) that were downregulated in bladder cancer were screened using the PPI network, survival curves, ROC curve analysis, and expression validation. Previous studies have shown that Vinculin (VCL) expression is downregulated in bladder cancer with tumor suppressor gene properties (40), which is consistent with the results of the present study. Wu et al. further confirmed that VCL can be used as a potential protein marker for bladder cancer (41) and provided strong support for the results of this study. Filament protein A (FLNA) plays a key role in the development of blood vessels, heart, and brain organs, particularly in the formation of intercellular contacts and adherent junctions (42). It has been demonstrated that the expression level of FLNA is significantly reduced in bladder cancer tissues and that FLNA overexpression can effectively inhibit the proliferation, invasion, and metastasis of bladder cancer cells (43), corroborating the results of the present study. Transglutamine protein (TAGLN) is an important actin-associated protein, which has been found to be expressed at significantly reduced levels in a variety of cancers, including prostate cancer. Based on this feature, TAGLN is widely considered to have tumor suppressor effects (44), and further studies have shown that TAGLN may not only serve as a potential diagnostic biomarker for bladder cancer (BCa), but may also be a promising therapeutic target (45). The actin alpha 2 (ACTA2, actin alpha 2) gene is primarily responsible for encoding smooth muscle alpha-actin, an essential component of the cytoskeleton, which plays a central role in key biological processes, such as cell contraction, migration, and signal transduction (46). Previous studies have shown that ACTA2 is one of the pivotal genes closely associated with the immune response in muscle invasive bladder cancer (MIBC) (47), which not only reveals the important role of ACTA2 in the development of MIBC but also provides a new research direction for an in-depth understanding of the mechanism of tumor immune microenvironment regulation. The collagen type VI alpha 2 chain (COL6A2) gene is an important component of the three alpha chains encoding collagen type VI. Recent studies have shown that COL6A2 plays a key role in many tumors, especially lung adenocarcinoma, and its expression level can be used as a prognostic predictor and a potential marker for targeted therapies (48). Notably, through an in-depth study of bladder cancer tissues, Piao et al. found that the mRNA expression level of COL6A2 was significantly downregulated in both non-muscle-invasive bladder cancer (NMIBC) and muscle-invasive bladder cancer (MIBC) tissue samples compared to that in normal bladder tissues (49), which is consistent with our findings. Caldesmon 1 (CALD1) is an important cytoskeleton-associated protein that affects cell morphology and migration capacity by regulating the dynamic balance of actin filaments (50). Recent studies have revealed that CALD1 not only serves as a prognostic biomarker for bladder cancer but may also contribute to its progression of bladder cancer by participating in the remodelling process of the tumor microenvironment. Further studies have shown that the expression of CALD1 is significantly correlated with immune cell infiltration in bladder cancer, which provides a new perspective for understanding the regulation of the tumor immune microenvironment (51).

Although the detection of downregulated genes in clinical practice presents certain challenges, they can still serve as reliable biomarkers, and their clinical feasibility has been supported by multiple lines of evidence. First, the principle of “loss indicates abnormality” suggests that the reduced expression of tumor suppressor genes itself can indicate pathological changes. This mechanism is exemplified by classical tumor suppressors such as TP53, the “guardian of the genome,” which regulates the cell cycle, DNA damage repair, and apoptosis (52), and whose dysfunction or downregulation is implicated in over 50% of human tumors (53). Although the downregulated genes in our study were not classical tumor suppressors, their reduced expression may similarly contribute to bladder cancer development, supporting their potential as diagnostic biomarkers or therapeutic targets. Consistently, previous studies have shown that downregulated genes can serve as prognostic biomarkers and immunotherapy targets. For example, SOX7 is downregulated in colorectal cancer, and its low expression correlates with poor prognosis (54). Second, the detection of downregulated genes can be achieved through several approaches, including: (1) highly sensitive RNA quantification techniques (e.g., RT-qPCR), which allow accurate discrimination even at low expression levels; (2) protein-level assessment of functional loss, such as immunohistochemistry (IHC) detection of PTEN loss, which is widely used as a prognostic indicator in prostate cancer (55, 56); and (3) epigenetic markers such as DNA methylation, which provide stable signals, exemplified by the FDA-approved colorectal cancer screening test Epi proColon®, based on SEPT9 promoter methylation (57). Collectively, these findings suggest that although downregulation itself may be difficult to capture directly, its associated transcriptomic, proteomic, and epigenetic alterations can be transformed into stable and detectable biological signals. Therefore, the six genes identified in our study as negatively associated with bladder tumors may represent promising diagnostic biomarkers and potential therapeutic targets for clinical application.

This study has made significant progress in constructing a B cell-based prognostic model for bladder cancer and its tumor microenvironment. However, this study had limitations. First, although the B cell prognostic model developed using public datasets demonstrated good predictive performance, its clinical utility requires validation through large-scale prospective clinical studies. Second, due to limitations in data availability, the sample size included in this study was relatively small; future multi-center, large-cohort studies are needed to further assess the robustness of the model. Finally, although we confirmed the expression levels of key prognostic genes and their associations with patient outcomes using public clinical data and RT-qPCR experiments, in-depth functional studies on these genes have not yet been conducted. Therefore, in future studies, we plan to employ systematic gene editing and functional assays to elucidate the molecular mechanisms of these key genes in bladder cancer progression, providing a stronger theoretical foundation for clinical applications.

5 Conclusion

In this study, we constructed single-cell transcriptome profiles of bladder cancer using scRNA-seq, and explored the role of B cells in cell trajectories, transcription factor regulatory networks, and intercellular communication mechanisms. By integrating single-cell samples with TCGA and GEO data, six immune cell-associated genes closely related to B cells, with significant differences, were successfully identified. These findings suggest six potential molecular targets for bladder cancer treatment. Although the specific mechanisms of action of these genes need to be validated by further studies, these results are important for an in-depth understanding of the molecular mechanisms of bladder cancer pathogenesis.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Ethics statement

The studies involving humans were approved by the Ethics Committee of the First People’s Hospital of Anning, Affiliated to Kunming University of Science and Technology. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Author contributions

LW: Conceptualization, Funding acquisition, Writing – original draft. JY: Data curation, Formal Analysis, Validation, Writing – original draft. ZX: Data curation, Formal Analysis, Validation, Writing – original draft. BT: Project administration, Writing – original draft. YH: Software, Writing – original draft. YZ: Software, Writing – original draft. JW: Formal Analysis, Writing – original draft. YM: Writing – original draft. ZZ: Conceptualization, Writing – review & editing. LY: Conceptualization, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research and/or publication of this article. This study was supported by Yunnan Fundamental Research Projects (KUST-AN2023009Y), Yunnan Fundamental Research Projects(202301BE070001-019), Medical Technology Centre Construction Project of Kunming Health Science and Technology Talent Cultivation Project(2022-SW(technology)-020), and Kunming Health Science and Technology Talent Cultivation Project Medical Science and Technology Discipline Leader Training Program (2023-SW (with head)-008).

Acknowledgments

We thank the Department of Urology and the Research Office of the First People’s Hospital of Anning, Affiliated to Kunming University of Science and Technology, for their support of this work. We also thank the server provided by RIZHIMANG for its assistance during the data analysis process.

Conflict of interest

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1600254/full#supplementary-material

References

1. Babjuk M, Burger M, Capoun O, Cohen D, Compérat EM, Dominguez Escrig JL, et al. European association of urology guidelines on non-muscle-invasive bladder cancer (Ta, T1, and carcinoma in situ). Eur Urol. (2022) 81:75–94. doi: 10.1016/j.eururo.2021.08.010

PubMed Abstract | Crossref Full Text | Google Scholar

2. Siegel RL, Kratzer TB, Giaquinto AN, Sung H, and Jemal A. Cancer statistics, 2025. CA: Cancer J Clin. (2025) 75:10–45. doi: 10.3322/caac.21871

PubMed Abstract | Crossref Full Text | Google Scholar

3. Genitsch V, Kollár A, Vandekerkhove G, Blarer J, Furrer M, Annala M, et al. Morphologic and genomic characterization of urothelial to sarcomatoid transition in muscle-invasive bladder cancer. Urol Oncol. (2019) 37:573.e19–.e29. doi: 10.1016/j.urolonc.2019.09.025

PubMed Abstract | Crossref Full Text | Google Scholar

4. Mohanty SK, Lobo A, Mishra SK, and Cheng L. Precision medicine in bladder cancer: present challenges and future directions. J Personalized Med. (2023) 13(5):756. doi: 10.3390/jpm13050756

PubMed Abstract | Crossref Full Text | Google Scholar

5. Minoli M, Kiener M, Thalmann GN, Kruithof-de Julio M, and Seiler R. Evolution of urothelial bladder cancer in the context of molecular classifications. Int J Mol Sci. (2020) 21(16):5670. doi: 10.3390/ijms21165670

PubMed Abstract | Crossref Full Text | Google Scholar

6. Ottley EC, Pell R, Brazier B, Hollidge J, Kartsonaki C, Browning L, et al. Greater utility of molecular subtype rather than epithelial-to-mesenchymal transition (EMT) markers for prognosis in high-risk non-muscle-invasive (HGT1) bladder cancer. J Pathol Clin Res. (2020) 6:238–51. doi: 10.1002/cjp2.167

PubMed Abstract | Crossref Full Text | Google Scholar

7. Wilson F, Joseph N, and Choudhury A. Biomarkers in muscle invasive bladder cancer. Adv Clin Chem. (2022) 107:265–97. doi: 10.1016/bs.acc.2021.07.005

PubMed Abstract | Crossref Full Text | Google Scholar

8. Kamoun A, de Reyniès A, Allory Y, Sjödahl G, Robertson AG, Seiler R, et al. A consensus molecular classification of muscle-invasive bladder cancer. Eur Urol. (2020) 77:420–33. doi: 10.1016/j.eururo.2019.09.006

PubMed Abstract | Crossref Full Text | Google Scholar

9. Lopez-Beltran A, Cookson MS, Guercio BJ, and Cheng L. Advances in diagnosis and treatment of bladder cancer. BMJ. (2024) 384:e076743. doi: 10.1136/bmj-2023-076743

PubMed Abstract | Crossref Full Text | Google Scholar

10. Compérat E, Amin MB, Cathomas R, Choudhury A, De Santis M, Kamat A, et al. Current best practice for bladder cancer: a narrative review of diagnostics and treatments. Lancet (London England). (2022) 400:1712–21. doi: 10.1016/S0140-6736(22)01188-6

PubMed Abstract | Crossref Full Text | Google Scholar

11. Lenis AT, Lec PM, and Chamie K. Bladder cancer. JAMA. (2020) 324:2006. doi: 10.1001/jama.2020.17601

PubMed Abstract | Crossref Full Text | Google Scholar

12. Siegel RL, Miller KD, Fuchs HE, and Jemal A. Cancer statistics, 2022. CA: Cancer J Clin. (2022) 72:7–33. doi: 10.3322/caac.21708

PubMed Abstract | Crossref Full Text | Google Scholar

13. Chen YP, Zhang Y, Lv JW, Li YQ, Wang YQ, He QM, et al. Genomic analysis of tumor microenvironment immune types across 14 solid cancer types: immunotherapeutic implications. Theranostics. (2017) 7:3585–94. doi: 10.7150/thno.21471

PubMed Abstract | Crossref Full Text | Google Scholar

14. Robertson AG, Kim J, Al-Ahmadie H, Bellmunt J, Guo G, Cherniack AD, et al. Comprehensive molecular characterization of muscle-invasive bladder cancer. Cell. (2017) 171:540–56.e25. doi: 10.1016/j.cell.2017.09.007

PubMed Abstract | Crossref Full Text | Google Scholar

15. Zhang E, Ding C, Li S, Zhou X, Aikemu B, Fan X, et al. Roles and mechanisms of tumour-infiltrating B cells in human cancer: a new force in immunotherapy. biomark Res. (2023) 11:28. doi: 10.1186/s40364-023-00460-1

PubMed Abstract | Crossref Full Text | Google Scholar

16. Chen H, Ye F, and Guo G. Revolutionizing immunology with single-cell RNA sequencing. Cell Mol Immunol. (2019) 16:242–9. doi: 10.1038/s41423-019-0214-4

PubMed Abstract | Crossref Full Text | Google Scholar

17. Ma KY, Schonnesen AA, Brock A, Van Den Berg C, Eckhardt SG, Liu Z, et al. Single-cell RNA sequencing of lung adenocarcinoma reveals heterogeneity of immune response-related genes. JCI Insight. (2019) 4(4):e121387. doi: 10.1172/jci.insight.121387

PubMed Abstract | Crossref Full Text | Google Scholar

18. Butler A, Hoffman P, Smibert P, Papalexi E, and Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. (2018) 36:411–20. doi: 10.1038/nbt.4096

PubMed Abstract | Crossref Full Text | Google Scholar

19. Love MI, Huber W, and Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | Crossref Full Text | Google Scholar

20. Stuart T and Satija R. Integrative single-cell analysis. Nat Rev Genet. (2019) 20:257–72. doi: 10.1038/s41576-019-0093-7

PubMed Abstract | Crossref Full Text | Google Scholar

21. Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol. (2019) 20:163–72. doi: 10.1038/s41590-018-0276-y

PubMed Abstract | Crossref Full Text | Google Scholar

22. Sherman BT, Hao M, Qiu J, Jiao X, Baseler MW, Lane HC, et al. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic Acids Res. (2022) 50:W216–w21. doi: 10.1093/nar/gkac194

PubMed Abstract | Crossref Full Text | Google Scholar

23. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. (2019) 10:1523. doi: 10.1038/s41467-019-09234-6

PubMed Abstract | Crossref Full Text | Google Scholar

24. Haghverdi L, Büttner M, Wolf FA, Buettner F, and Theis FJ. Diffusion pseudotime robustly reconstructs lineage branching. Nat Methods. (2016) 13:845–8. doi: 10.1038/nmeth.3971

PubMed Abstract | Crossref Full Text | Google Scholar

25. Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. (2014) 32:381–6. doi: 10.1038/nbt.2859

PubMed Abstract | Crossref Full Text | Google Scholar

26. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. (2021) 12:1088. doi: 10.1038/s41467-021-21246-9

PubMed Abstract | Crossref Full Text | Google Scholar

27. Langfelder P and Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf. (2008) 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | Crossref Full Text | Google Scholar

28. Shen L, Mo J, Yang C, Jiang Y, Ke L, Hou D, et al. SurvivalPath:A R package for conducting personalized survival path mapping based on time-series survival data. PloS Comput Biol. (2023) 19:e1010830. doi: 10.1371/journal.pcbi.1010830

PubMed Abstract | Crossref Full Text | Google Scholar

29. Heagerty PJ and Zheng Y. Survival model predictive accuracy and ROC curves. Biometrics. (2005) 61:92–105. doi: 10.1111/j.0006-341X.2005.030814.x

PubMed Abstract | Crossref Full Text | Google Scholar

30. Segerstolpe Å, Palasantza A, Eliasson P, Andersson EM, Andréasson AC, Sun X, et al. Single-cell transcriptome profiling of human pancreatic islets in health and type 2 diabetes. Cell Metab. (2016) 24:593–607. doi: 10.1016/j.cmet.2016.08.020

PubMed Abstract | Crossref Full Text | Google Scholar

31. Boxer E, Feigin N, Tschernichovsky R, Darnell NG, Greenwald AR, Hoefflin R, et al. Emerging clinical applications of single-cell RNA sequencing in oncology. Nat Rev Clin Oncol. (2025) 22:315–26. doi: 10.1038/s41571-025-01003-3

PubMed Abstract | Crossref Full Text | Google Scholar

32. Llaó-Cid L. Single-cell insights into the role of T cells in B-cell Malignancies. FEBS Lett. (2025). doi: 10.1002/1873-3468.70099

PubMed Abstract | Crossref Full Text | Google Scholar

33. Liu M, Sun Q, Wang J, Wei F, Yang L, and Ren X. A new perspective: Exploring future therapeutic strategies for cancer by understanding the dual role of B lymphocytes in tumor immunity. Int J Cancer. (2019) 144:2909–17. doi: 10.1002/ijc.31850

PubMed Abstract | Crossref Full Text | Google Scholar

34. Li Q, Teitz-Tennenbaum S, Donald EJ, Li M, and Chang AE. In vivo sensitized and in vitro activated B cells mediate tumor regression in cancer adoptive immunotherapy. J Immunol. (2009) 183:3195–203. doi: 10.4049/jimmunol.0803773

PubMed Abstract | Crossref Full Text | Google Scholar

35. Zirakzadeh AA, Marits P, Sherif A, and Winqvist O. Multiplex B cell characterization in blood, lymph nodes, and tumors from patients with Malignancies. J Immunol. (2013) 190:5847–55. doi: 10.4049/jimmunol.1203279

PubMed Abstract | Crossref Full Text | Google Scholar

36. Jiang K, Liu H, Xie D, and Xiao Q. Differentially expressed genes ASPN, COL1A1, FN1, VCAN and MUC5AC are potential prognostic biomarkers for gastric cancer. Oncol Lett. (2019) 17:3191–202. doi: 10.3892/ol.2019.9952

PubMed Abstract | Crossref Full Text | Google Scholar

37. Wang S, Gao B, Yang H, Liu X, Wu X, and Wang W. MicroRNA-432 is downregulated in cervical cancer and directly targets FN1 to inhibit cell proliferation and invasion. Oncol Lett. (2019) 18:1475–82. doi: 10.3892/ol.2019.10403

PubMed Abstract | Crossref Full Text | Google Scholar

38. Zhang L, Wang Y, Song M, Chang A, Zhuo W, and Zhu Y. Fibronectin 1 as a key gene in the genesis and progression of cadmium-related bladder cancer. Biol Trace Element Res. (2023) 201:4349–59. doi: 10.1007/s12011-022-03510-1

PubMed Abstract | Crossref Full Text | Google Scholar

39. Sahana J, Corydon TJ, Wehland M, Krüger M, Kopp S, Melnik D, et al. Alterations of growth and focal adhesion molecules in human breast cancer cells exposed to the random positioning machine. Front Cell Dev Biol. (2021) 9:672098. doi: 10.3389/fcell.2021.672098

PubMed Abstract | Crossref Full Text | Google Scholar

40. Wu H, Jiang W, Ji G, Xu R, Zhou G, and Yu H. Exploring microRNA target genes and identifying hub genes in bladder cancer based on bioinformatic analysis. BMC Urol. (2021) 21:90. doi: 10.1186/s12894-021-00857-w

PubMed Abstract | Crossref Full Text | Google Scholar

41. Wu Z, Wang S, Jiang F, Li Q, Wang C, Wang H, et al. Mass spectrometric detection combined with bioinformatic analysis identified possible protein markers and key pathways associated with bladder cancer. Gene. (2017) 626:407–13. doi: 10.1016/j.gene.2017.05.054

PubMed Abstract | Crossref Full Text | Google Scholar

42. Song M, He Q, Berk BA, Hartwig JH, Stossel TP, and Nakamura F. An adventitious interaction of filamin A with RhoGDI2(Tyr153Glu). Biochem Biophys Res Commun. (2016) 469:659–64. doi: 10.1016/j.bbrc.2015.12.044

PubMed Abstract | Crossref Full Text | Google Scholar

43. Wang Z, Li C, Jiang M, Chen J, Yang M, and Pu J. Filamin A (FLNA) regulates autophagy of bladder carcinoma cell and affects its proliferation, invasion and metastasis. Int Urol Nephrol. (2018) 50:263–73. doi: 10.1007/s11255-017-1772-y

PubMed Abstract | Crossref Full Text | Google Scholar

44. Wen F, Sun X, Sun C, Dong Z, Jia G, Bao W, et al. TAGLN is downregulated by TRAF6-mediated proteasomal degradation in prostate cancer cells. Mol Cancer Res: MCR. (2021) 19:1113–22. doi: 10.1158/1541-7786.MCR-20-0513

PubMed Abstract | Crossref Full Text | Google Scholar

45. Liu Y, Wu X, Wang G, Hu S, Zhang Y, and Zhao S. CALD1, CNN1, and TAGLN identified as potential prognostic molecular markers of bladder cancer by bioinformatics analysis. Medicine. (2019) 98:e13847. doi: 10.1097/MD.0000000000013847

PubMed Abstract | Crossref Full Text | Google Scholar

46. Guo DC, Papke CL, Tran-Fadulu V, Regalado ES, Avidan N, Johnson RJ, et al. Mutations in smooth muscle alpha-actin (ACTA2) cause coronary artery disease, stroke, and Moyamoya disease, along with thoracic aortic disease. Am J Hum Genet. (2009) 84:617–27. doi: 10.1016/j.ajhg.2009.04.007

PubMed Abstract | Crossref Full Text | Google Scholar

47. Chen Z, Liu G, Liu G, Bolkov MA, Shinwari K, Tuzankina IA, et al. Defining muscle-invasive bladder cancer immunotypes by introducing tumor mutation burden, CD8+ T cells, and molecular subtypes. Hereditas. (2021) 158:1. doi: 10.1186/s41065-020-00165-7

PubMed Abstract | Crossref Full Text | Google Scholar

48. Kong F, Lu Z, Xiong Y, Zhou L, and Ye Q. A novel cancer-associated fibroblasts risk score model predict survival and immunotherapy in lung adenocarcinoma. Mol Genet Genom: MGG. (2024) 299:70. doi: 10.1007/s00438-024-02156-z

PubMed Abstract | Crossref Full Text | Google Scholar

49. Piao XM, Hwang B, Jeong P, Byun YJ, Kang HW, Seo SP, et al. Collagen type VI−α1 and 2 repress the proliferation, migration and invasion of bladder cancer cells. Int J Oncol. (2021) 59(1):37. doi: 10.3892/ijo.2021.5217

PubMed Abstract | Crossref Full Text | Google Scholar

50. Nie S, Kee Y, and Bronner-Fraser M. Caldesmon regulates actin dynamics to influence cranial neural crest migration in Xenopus. Mol Biol Cell. (2011) 22:3355–65. doi: 10.1091/mbc.e11-02-0165

PubMed Abstract | Crossref Full Text | Google Scholar

51. Du Y, Jiang X, Wang B, Cao J, Wang Y, Yu J, et al. The cancer-associated fibroblasts related gene CALD1 is a prognostic biomarker and correlated with immune infiltration in bladder cancer. Cancer Cell Int. (2021) 21:283. doi: 10.1186/s12935-021-01896-x

PubMed Abstract | Crossref Full Text | Google Scholar

52. Aubrey BJ, Kelly GL, Janic A, Herold MJ, and Strasser A. How does p53 induce apoptosis and how does this relate to p53-mediated tumour suppression? Cell Death Differ. (2018) 25:104–13. doi: 10.1038/cdd.2017.169

PubMed Abstract | Crossref Full Text | Google Scholar

53. Guan YS, He Q, and Zou Q. Status quo of p53 in the treatment of tumors. Anti-cancer Drugs. (2016) 27:811–8. doi: 10.1097/CAD.0000000000000397

PubMed Abstract | Crossref Full Text | Google Scholar

54. Zhang J, Sun Z, Li G, Ding L, Wang Z, and Liu M. Discovering biomarkers associated with infiltration of CD8(+) T cells and tumor-associated fibrosis in colon adenocarcinoma using single-cell RNA sequencing and gene co-expression network. Front Immunol. (2025) 16:1496640. doi: 10.3389/fimmu.2025.1496640

PubMed Abstract | Crossref Full Text | Google Scholar

55. Lotan TL, Wei W, Ludkovski O, Morais CL, Guedes LB, Jamaspishvili T, et al. Analytic validation of a clinical-grade PTEN immunohistochemistry assay in prostate cancer by comparison with PTEN FISH. Modern Pathol. (2016) 29:904–14. doi: 10.1038/modpathol.2016.88

PubMed Abstract | Crossref Full Text | Google Scholar

56. Lotan TL, Heumann A, Rico SD, Hicks J, Lecksell K, Koop C, et al. PTEN loss detection in prostate cancer: comparison of PTEN immunohistochemistry and PTEN FISH in a large retrospective prostatectomy cohort. Oncotarget. (2017) 8:65566–76. doi: 10.18632/oncotarget.19217

PubMed Abstract | Crossref Full Text | Google Scholar

57. Song L, Jia J, Peng X, Xiao W, and Li Y. The performance of the SEPT9 gene methylation assay and a comparison with other CRC screening tests: A meta-analysis. Sci Rep. (2017) 7:3032. doi: 10.1038/s41598-017-03321-8

PubMed Abstract | Crossref Full Text | Google Scholar

Keywords: bladder cancer, B cells, immune microenvironment, single-cell RNA-seq analysis, biomarkers, RT-qPCR

Citation: Wang L, Yang J, Xu Z, Tao B, He Y, Zhao Y, Wu J, Ma Y, Zhong Z and Ye L (2025) Comprehensive analysis of single-cell and bulk transcriptomes reveals key B-cell genes and immune microenvironment regulation in bladder cancer. Front. Immunol. 16:1600254. doi: 10.3389/fimmu.2025.1600254

Received: 26 March 2025; Accepted: 25 September 2025;
Published: 17 October 2025.

Edited by:

Ana María Hurtado López, Universidad Católica San Antonio de Murcia (UCAM), Spain

Reviewed by:

Peng Xu, Yale University, United States
Xin Mu, Tianjin University, China

Copyright © 2025 Wang, Yang, Xu, Tao, He, Zhao, Wu, Ma, Zhong and Ye. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Zitao Zhong, enppdGFvQGt1c3QuZWR1LmNu; Lin Ye, NzkzNzU4OUBxcS5jb20=

These authors have contributed equally to this work

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.