High-throughput sequencing data and clinical information for bladder cancer (BLCA) were obtained from the UCSC Xena platform (https://xenabrowser.net/) (19).The paired single-cell datasets, including normal bladder tissues and cancer tissues from GSE222315, were acquired via the GEO database. To ensure data quality and standardization, batch effect correction and integration of multi-group gene expression data were performed using R packages such as “limma” and “sva”. For single-cell data, the Harmony integration algorithm was applied to eliminate potential experimental batch variations. Additionally, a taurine metabolism-related gene set was systematically screened from the GeneCard database, serving as the foundation for subsequent gene scoring and functional enrichment analyses.

Preprocessing of single-cell transcriptome data was conducted using the “Seurat” R package. Low-quality cells were first detected via the “PercentageFeatureSet” function, with single-cell gene expression set to a range of 200–10,000 and mitochondrial gene expression ratio below 20%. Data normalization was then performed using the “NormalizeData” function, converting processed data into Seurat objects for subsequent analysis. Highly variable genes were identified via the built-in “FindVariableFeatures” function, followed by scaling and principal component analysis (PCA) using the “RunPCA” tool. To reveal latent structural features, dimensionality reduction and visualization were further performed using UMAP and t-SNE algorithms. In cell subset analysis, differential expression analysis of distinct cell types was conducted via “FindAllMarkers” combined with Wilcoxon rank-sum tests and manual annotation based on data source annotations, with cell subset annotation completed using feature genes. For assessing purine metabolism characteristics, the aforementioned taurine metabolism-related gene set was integrated with five classic analytical methods—AddModuleScore, single-sample gene set enrichment analysis (ssGSEA), AUCell, UCell, and singscore—to score the endogenous tumor microenvironment of BLCA, with metabolic activity characteristics of each cell subset ultimately visualized via interactive bubble plots.

Pseudotime trajectory analysis of tumor cells was performed using the “Monocle 3” R package based on dynamic changes in single-cell gene expression data. Specifically, dimensionality reduction was first conducted via the DDRTree algorithm, followed by pseudotime analysis to reveal continuous changes in epithelial cells with high/low taurine expression from initial to terminal states. The “PLOT_CELL_TRACTURE” function was used to dynamically order distinct cell subsets along the pseudotime axis. For intercellular communication analysis, the “CellChat” software was employed to construct ligand-receptor regulatory networks across cell types. Integrating interaction data from CellPhoneDB and STRING databases via the netVisualDiffInteraction algorithm with cell-specific expression profile data, bidirectional regulatory mechanisms between high/low taurine-expressing epithelial cells and immune/stromal cells were revealed, providing a theoretical basis for targeted tumor microenvironment therapy.

The spatial transcriptome dataset GSE171351 from the GEO database was selected as an independent validation cohort. These datasets included fresh-frozen tumor tissue sections from which high-resolution spatial gene expression matrices were obtained via the 10x Genomics Visium platform. During data preprocessing, normalization and correction were performed using the “SCTransform” algorithm in “Seurat”, followed by dimensionality reduction of standardized UMI count matrices via the “RunPCA” function to identify spatial transcriptomic subsets with similar gene expression patterns. For spatial location information in spatial transcriptome data, deep mining was performed using “stLearn”, a tool developed based on Scanpy. This method integrated cell type annotation, gene co-expression networks, and spatial coordinate information to identify significantly enriched gene modules via spatial autocorrelation analysis. Cell type deconvolution of spatial transcriptome data was performed using a Bayesian hierarchical model based on single-cell reference atlases, effectively distinguishing spatial distribution characteristics of epithelial cells, immune cells, and stromal cells with high/low taurine metabolic features and revealing interaction patterns among cell types in the tumor microenvironment.

Based on the TCGA-BLCA cohort, the dataset was stratified and randomly sampled at a 7:3 ratio into training and validation sets. Single-factor Cox proportional hazards regression was first performed via survival analysis R packages to systematically screen genes significantly associated with overall survival (OS, p<0.05). To optimize model complexity, Lasso regression was used for feature dimensionality reduction, with the optimal λ value selected via 1000-fold cross-validation to further screen model-building genes. Combined with multi-factor stepwise regression analysis, key oncogenes were identified to construct a multivariate Cox risk model. Based on this model, a taurine metabolic score (TM score) was calculated for each patient. Survival analysis was performed via the Kaplan-Meier method, with intergroup differences evaluated via log-rank tests. Model predictive efficacy was assessed via time-dependent ROC curves. The final prognostic nomogram integrated three independent variables: TM score, age, and clinical stage. Model stability was validated via Bootstrap resampling, with calibration curves showing high consistency between predicted and observed survival rates.

The ESTIMATE algorithm was used to calculate stromal score, immune score, and tumor purity for each patient based on expression matrices of 141 stromal and 225 immune signature genes, quantifying tumor microenvironment characteristics by inferring proportions of tumor mesenchyme and immune cells from gene expression data. To systematically evaluate immune infiltration characteristics of BLCA patients stratified by TM score, analysis results from seven mainstream algorithms—including CIBERSORT, MCPcounter, and XCELL—were obtained from the TIMER 2.0 database (20).Activation states of 28 immune cell types were further quantified via single-sample gene set enrichment analysis (ssGSEA). For immune escape assessment, the TIDE computational framework (http://tide.dfci.harvard.edu/) was used to analyze immunotherapy response potential in tumor samples and evaluate immune checkpoint inhibitor efficacy. Multi-dimensional charts were generated using ggplot2 throughout the study.

Bladder cancer cell line T24 and human bladder tumor-associated fibroblasts hCAF-BT (both purchased from the Chinese Academy of Sciences Cell Bank) were used in this experimental section. T24 cells were cultured in high-glucose DMEM medium with 10% FBS, while hCAF-BT cells were maintained in F12K medium supplemented with 1% PS double-antibiotics. Cells were cultured in a 37°C, 5% CO₂ incubator and passaged every 2–3 days. For Transwell co-culture system construction, log-phase T24 cells were seeded in the upper chamber and hCAF-BT cells in the lower chamber. After 24-hour synchronization with serum-free medium, cultures were switched to 1% FBS medium and maintained for 72 hours.

A lentiviral vector expressing shRNA targeting the FAAH gene (designed by GenePharma) and a negative control sh-NC were used. T24 cells were seeded at 5×10⁵ cells/well in 6-well plates, and transfection was performed using Lipofectamine 3000 when cell confluency reached 60%. Forty-eight hours post-transfection, puromycin (2μg/mL) was added for resistance screening, with stable knockdown cell lines obtained after 14 days of continuous culture.

T24 cells were adjusted to a single-cell suspension at 5×10³ cells/mL and seeded at 100μL/well in 96-well plates with 6 replicates per group. At 1, 2, 3, 4, 5 and 6 days post-seeding, 10μL CCK-8 reagent (Dojindo) was added, and after 2-hour incubation in the dark, absorbance at 450nm was measured using a microplate reader. Cell growth curves were plotted with time on the x-axis and OD values on the y-axis, with growth trends fitted using a four-parameter Logistic equation.

All data calculations and statistical analyses were performed in a 64-bit R (4.2 version) environment, with major R packages used described in the methods section. Data analysis workflows strictly adhered to FAIR principles. Statistical significance was denoted by * for P<0.05, ** for P<0.01, and *** for P<0.001. For multiple comparison scenarios (e.g., GSEA analysis or intergroup difference tests), Bonferroni correction was used by default.

To characterize the BLCA microenvironment, paired samples from six BLCA patients were subjected to rigorous quality control and dimensionality reduction, eliminating noise and redundant data ( Figures 1A, B ). The six samples mentioned above were further reduced in dimension to form 18 cell clusters. Subsequently, single-cell annotation techniques were employed to reveal that single-cell clustering had identified distinct cell populations, including mononuclear macrophages, plasma cells, epithelial cells, T/NK cells, B cells, endothelial cells, and fibroblasts ( Figures 1C, D ). The balanced distribution of these cell types in BLCA suggests a dynamic equilibrium critical for TME maintenance and progression. Marker gene analysis revealed cell-type-specific expression signatures ( Figure 1E ). Leveraging five single-cell scoring algorithms (UCell, AddModuleScore, AUCell, GSVA, singscore), we quantified taurine metabolic activity across cell types. Endothelial cells, macrophages, and epithelial cells exhibited the highest taurine metabolic scores (TMscore), as visualized in interactive bubble plots ( Figures 1F, G ). Comparative analysis between tumor and normal tissues highlighted differential TMscore in epithelial cells, fibroblasts, and macrophages, indicating tumor-specific reprogramming of taurine metabolism ( Figure 1H ).

Using pseudotime analysis via Monocle 3, we reconstructed differentiation trajectories of TMscore-defined epithelial cell clusters. Eight distinct epithelial subpopulations were identified, with Cluster 6 displaying terminal differentiation markers at the trajectory endpoint ( Figure 2A ). Stratifying epithelial cells by TMscore median revealed enrichment of Notch signaling and EGFR tyrosine kinase inhibitor resistance pathways in high-TMscore cells ( Figure 2B ). Transcription factor binding motif analysis identified HEY2 and PPARD as key regulators of intercellular crosstalk, with their binding sites enriched in differentially expressed gene promoters ( Figure 2C ). Ligand-receptor interaction modeling using CellChat uncovered enhanced communication between high-TMscore epithelial cells and fibroblasts/macrophages. The CCL28-CCR10 and WNT7B-FZD1 signaling axes emerged as dominant pathways mediating this crosstalk ( Figures 2D, E ), serving as molecular hubs for intercellular signaling.

Integrating single-cell and spatial transcriptomics, we mapped 13 distinct cell subpopulations in BLCA tissue sections ( Figure 3A ). Spatial distribution of taurine metabolic genes (SLC27A5, CAV1) revealed pronounced heterogeneity, with tumor core regions exhibiting hyperactivation of taurine metabolism ( Figures 3B, C ). In the tumor microenvironment, metabolic crosstalk often has a systemic impact, where changes in one pathway can affect the entire network. To simultaneously analyze alterations in different metabolic pathways within the current spatial transcriptomics sections, we performed a metabolic pathway enrichment analysis via scMetabolism, which identified Cluster 11 as a metabolic hotspot ( Figure 3D ). RCTD-based spatial mapping showed high-TMscore epithelial cells preferentially localized to niches enriched in fibroblasts and macrophages ( Figures 3E–G ), validated by ligand-receptor co-expression and proximity ligation assays. PAR signaling pathway networks were identified as critical mediators of spatial cell–cell interactions in high-TMscore regions ( Figure 3H ), underscoring the role of metabolic reprogramming in shaping TME architecture.

A taurine metabolic dysregulation index (TMs) was developed through integrative analysis of differentially expressed genes and survival data. Using Lasso-Cox regression, a multi-gene prognostic model was constructed, dividing patients into high- and low-TMs groups based on median scores ( Figures 4A–C ). Survival analysis revealed significantly shorter overall survival (OS) in high-TMs cohorts across TCGA-BLCA datasets ( Figures 4D, E ). Time-dependent ROC curves demonstrated superior prognostic accuracy of TMs compared to conventional clinicopathological features (AUC: 0.78 vs. 0.62, Figures 4K, L ). Multivariate Cox analysis confirmed TMs as an independent prognostic factor (HR=1.87, P<0.001, Figures 4H, I ). A clinical nomogram integrating TMs, age, and stage showed high calibration accuracy for 1-, 3-, and 5-year OS prediction ( Figures 4G, M, N ). Genomic profiling revealed high-TMs tumors harbored increased missense mutations, particularly in TP53, TTN, and RYR2 ( Figures 5A, D ). Combined TMs-TMB analysis identified a high-TMs/low-TMB subgroup with poorest prognosis ( Figures 5B, C ), suggesting synergistic effects of metabolic dysregulation and immune evasion. No significant association was observed between TMs and MSI, emphasizing TMB as the primary genomic correlate ( Figure 5E ).

Emerging evidence has highlighted the intricate crosstalk between metabolic reprogramming and the TME in BLCA. Building upon our previous findings implicating taurine metabolic dysregulation in shaping TME composition, we systematically profiled immune cell infiltration across TMs subgroups. Employing a suite of deconvolution algorithms, including CIBERSORT, TIMER, and xCELL, we observed that high-TMs tumors were characterized by marked enrichment of immune populations such as natural killer (NK) cells, memory B cells, CD8⁺ T cells, regulatory T cells (Tregs), and M2-polarized macrophages ( Figure 6A ). The coexistence of both immunostimulatory and immunosuppressive subsets suggests a complex immune reprogramming landscape accompanying taurine metabolic alterations. To further delineate the immunological shifts associated with TMs status, we analyzed the expression of immune-related genes, revealing heightened immune gene activity in the high-TMs group ( Figure 6B ). Single-sample gene set enrichment analysis (ssGSEA) corroborated these findings, showing significant enrichment of macrophage- and helper T cell-associated signatures in the high-TMs cohort ( Figures 6C, D ). In parallel, GSVA of canonical tumor immune pathways revealed increased activity in antigen presentation, immune checkpoint signaling, and cytotoxic effector functions in high-TMs tumors ( Figure 6E ), supporting a metabolically driven immunological remodeling. Quantitative assessment using the ESTIMATE algorithm further confirmed a consistent reduction in stromal, immune, and overall ESTIMATE scores in high-TMs tumors ( Figure 6F ), reflecting a shift toward an immunosuppressive and metabolically reprogrammed TME. Collectively, these data implicate taurine metabolism as a pivotal modulator of immune contexture in BLCA.

Given the central role of the immune milieu in determining response to ICB, we next interrogated the impact of taurine metabolic imbalance on immunotherapy-related parameters. Expression profiling revealed that nearly all immune checkpoint genes were significantly upregulated in high-TMs tumors( Figure 7A ), consistent with an immunosuppressive, potentially therapy-resistant state.To further quantify immune evasion and dysfunction, we employed the Tumor Immune Dysfunction and Exclusion (TIDE) framework. High-TMs tumors displayed elevated immune exclusion scores and reduced dysfunction scores ( Figures 7B, C ), suggesting active immune escape despite preserved effector function. GSVA analyses demonstrated concurrent upregulation of multiple immune response pathways, reinforcing the presence of a metabolically reinforced, yet immunologically active, tumor state ( Figures 7D, E ). Stratification by the Immuno-pheno-score (IPS)—a composite predictor of ICB response—revealed that TMs risk groups exhibited distinct sensitivities to PD-1 and CTLA-4 blockade, with notable differences in CTLA-4 responsiveness ( Figures 7F–H ). Correlation analyses confirmed that the TMs-derived risk score positively associated with multiple features of the cancer-immunity cycle, including tumor antigen release, presentation, and immune infiltration ( Figure 7I ). These findings underscore the potential of TMs status as a biomarker for predicting ICB efficacy and guiding therapeutic stratification in BLCA.

Considering the value of TMs mentioned above for tumor immunity, we attempted to address this potential relationship through internal mechanisms. We constructed an ssGSEA score using taurine disorder criteria and evaluated the relationship between this criterion and common cytokines in the BLCA population. Cross-cohort correlation analyses demonstrated a robust association between TMs and the chemokine CCL15 ( Figure 8A ), a fibroblast-and tumor-derived factor implicated in immune recruitment and suppression. Given its high loading score in the taurine metabolic signature, FAAH emerged as a candidate effector of metabolic reprogramming in BLCA. These results suggest that FAAH may mediate immune modulation via paracrine signaling. In immunotherapy-treated cohorts, high FAAH expression was significantly associated with non-responder status (NR), and showed predictive value across clinical datasets ( Figures 8B–E ). The scatter plot reveals the correlation between FAAH and CCL15 in TCGA data ( Figure 8F ). Notably, conditioned media (CM) derived from CAFs exhibited significantly higher levels of CCL15 compared to that from normal fibroblasts ( Figure 8G ), implicating a FAAH-CCL15 axis in fibroblast-mediated immune modulation. To mechanistically validate these findings, we established a series of co-culture systems: control, T24 + CAF CM, T24 + CAF CM with FAAH knockdown (CAF-sh FAAH), and direct T24-CAF co-culture. CCK-8 demonstrated that CAF-derived CM enhanced T24 proliferation, an effect abrogated by FAAH knockdown ( Figures 8H, I ). Together, these data position FAAH as a functional driver of taurine metabolic reprogramming and an actionable target for reversing fibroblast-mediated immunosuppression in BLCA.