The Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) was used to mine the transcriptomic sequencing data in the GSE32894 (GPL6947) and GSE13507 (GPL6102) datasets related to BLCA (Sjodahl et al., 2012; Lee et al., 2010). The dataset GSE32894 consisted of a total of 308 BLCA tissue samples, of which 224 samples with survival information were specifically chosen for the validation of the prognostic model. A total of 67 control samples, including 58 paracancerous controls and 9 normal controls, along with 165 primary BLCA samples from GSE13507, were selected to validate prognostic gene expression. The gene expression matrix and clinical traits of The Cancer Genome Atlas (TCGA)-BLCA were retrieved from the University of California, Santa Cruz (UCSC) Xena (https://xena.ucsc.edu/). A total of 430 samples were acquired from TCGA-BLCA, of which 405 cancer tissue samples (400 cases with available survival information) and 19 adjacent tissue samples were selected as the training set. The MitoCarta3.0 database (http://www.broadinstitute.org/mitocarta) offered the acquisition of 1,136 mitochondrial-related genes (MRGs). We also identified a total of 18 distinct patterns of PCD and key regulatory genes from literature, encompassing 1,964 genes and including apoptosis, autophagy, and alkaliptosis. After eliminating duplicates, we ultimately acquired 1,548 PCD-related genes (PCD-RGs) (Qin et al., 2023).

To identify the differences in gene expression between BLCA tissues and control tissues, thereby laying a foundation for subsequent screening of prognostic genes, the differentially expressed genes (DEGs) between BLCA and control samples in the TCGA-BLCA dataset were identified by the DESeq2 package (version 1.38.3) (|log₂ FC| ≥ 0.5, p.adj <0.05) (Love et al., 2014). The results were visualized through a volcano plot and heatmap, created using the ggplot2 package (version 3.3.6) and the heatmap3 package (version 1.1.9), respectively (Gustavsson et al., 2022; Zhao et al., 2014).

The intersection genes of DEGs, MRGs, and PCD-RGs were obtained using the UpSetR package (version 1.4.0) and applied for subsequent analysis (Conway et al., 2017) to understand the biological functions and pathways related to intersection genes. The intersection genes were subjected to ID conversion using the clusterProfiler package (version 4.4.4). Then, functional annotations were made, focusing on the analysis of the Gene Ontology (GO) function and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (p.adj <0.05) (Wu et al., 2021). Visualization was enabled using the ggplot2 package (version 3.3.6). The STRING database (https://cn.string-db.org/) was utilized to explore the protein interactions among the intersecting genes, thereby constructing a protein–protein interaction (PPI) network (minimum required interaction score = 0.4).

The study aimed to identify prognostic genes significantly associated with the survival outcomes of patients with BLCA from the previously screened cross-genes and to construct a risk assessment model based on these genes. In TCGA-BLCA, the survival package (version 3.4-0) was utilized to conduct univariate Cox regression analysis on intersection genes, employing p < 0.05 and HR≠1 as criteria for identifying genes significantly associated with survival outcomes (https://CRAN.R-project.org/package=survival). Subsequently, a proportional hazards (PH) assumption test was performed on the results of univariate Cox analysis to confirm genes (p > 0.05). Stepwise regression method was employed to further identify prognostic genes, and multivariate Cox analysis was employed to construct risk models based on the prognostic genes using the formula: r i s k   s c o r e = ∑ i = 1 i c o e f i * e x p i . The “coefi” represents the prognostic gene coefficients, and “expi” indicates the expression situation of prognostic genes.

To evaluate the predictive performance of the BLCA prognostic model constructed based on the genes screened in the earlier stage and verify its generalizability, the BLCA samples from TCGA-BLCA were utilized for evaluating the prognostic model, while those from GSE32894 were employed for validating the prognostic model. The risk score of the samples was calculated using the formula, and based on its optimal threshold, the samples were divided into high- and low-risk categories. Then, the Kaplan–Meier (KM) curves for the two datasets were separately plotted by applying the survival package (version 3.4-0). The survival receiver operating characteristic (ROC) curves for 1–3 years were separately generated employing the timeROC package (version 0.4) in R (https://CRAN.R-project.org/package=survival, doi:10.1093). The area under the curve (AUC) >0.6 indicated that the model had a relatively good predictive performance for the survival status of the samples.

The study aimed to clarify the association between risk scores and clinical characteristics of BLCA, identify independent factors influencing patient prognosis, and construct a tool for intuitively predicting patients’ survival outcomes. In TCGA-BLCA, the proportion of risk scores and the significance of the differences between groups for the five clinical traits (age, gender, pathologic_T, grade, and stage) were analyzed using the Wilcox test for two groups and the Kruskal test for multiple groups. The traits pathologic_M and pathologic_N were excluded from the analysis due to a higher number of missing values. Therefore, the significance of the risk score and five clinical traits on prognosis was analyzed by implementing univariate Cox regression, followed by performing a PH hypothesis test on significant traits. Subsequently, independent prognostic factors were obtained by multivariate Cox regression analysis. Finally, the nomogram for predicting the survival rate of patients from 1 to 3 years was created using the survival and rms packages (version 6.3-0) (Xu et al., 2023). Additionally, an evaluation of the nomogram was conducted using calibration curves and decision curve analysis (DCA).

Finally, gene set enrichment analysis (GSEA) was conducted using the clusterProfiler package (version 4.4.4) to explore differences in relevant biological functions and involved signaling pathways between two risk categories (Wu et al., 2021; Zhou et al., 2024). First, a differential analysis was conducted on high- and low-risk categories. The DEGs between the two risk categories were then sorted and enriched utilizing the GSEA function of the R language (p.adj <0.05). The Hallmark and KEGG gene sets were downloaded as background sets implementing the msigdbr package (version 7.5.1) (https://CRAN.R-project.org/package=msigdbr).

In order to analyze the variation in immune infiltration between two risk categories, we implemented the ssGSEA algorithm provided by the GSVA package (version 1.44.5) to assess the enrichment scores of 28 immune cells in all samples within the TCGA-BLCA dataset, and then a Wilcoxon test was utilized to compare the variations in cell content between these two categories (p < 0.05) (Hanzelmann et al., 2013). Furthermore, Spearman correlation was applied to analyze the correlation between prognostic genes and risk scores with discrepant immune cell enrichment scores. In addition, the common immune checkpoint expressions of the two risk categories were compared, covering IDO1, PD-1, PD-L1, CD27, PD-L2, LAG3, TIM-3, ICOS, CTLA4, and LGALS9. The tumor immune dysfunction and exclusion (TIDE), dysfunction, and exclusion scores were computed for the samples in TCGA-BLCA, and their correlation with risk scores was computed to evaluate the immunotherapy efficacy of patients.

The between-group differences in IC₅₀ values of 198 drugs from Genomics of Drug Sensitivity in Cancer (GDSC, http://www.cancerrxgene.org/) were computed by the oncoPredict package (version 0.2), and the correlation between discrepant drugs and risk scores was assessed (p < 0.01) (Maeser DGruener RFHuang RS, 2021). The expressions of prognostic genes were definitively validated in both the TCGA-BLCA and GSE13507 datasets and in clinical samples. A total of 10 tissue samples were obtained from Shanxi Bethune Hospital, consisting of five control specimens and five BLCA specimens. Prior to sample collection, all participants provided signed informed consent forms, and the study received ethical approval from the Clinical Research Ethics Committee of Shanxi Bethune Hospital (SBQLL-2024-154).

The extraction and purification of total RNA were carried out using TRIzol (Ambion, Austin, United States), following the manufacturer’s instructions. Subsequently, cDNA was synthesized from the RNA using the SureScript First-Strand cDNA Synthesis Kit (ServiceBio, Wuhan, China). The resulting cDNA was then subjected to analysis using 2x Universal Blue SYBR Green qPCR Master Mix (ServiceBio) along with specific primer sequences (Additional file 1). Real-time quantitative polymerase chain reaction (RT-qPCR) amplification was performed under the following thermal cycling conditions: initial denaturation at 95 °C for 1 min, followed by 40 cycles of denaturation at 95 °C for 20 s, annealing at 55 °C for 20 s, and extension at 72 °C for 30 s. Data analysis employed the method known as 2^(-ΔΔCt), where GAPDH was used as an internal reference for normalization.

Data processing and analysis were carried out using R software (version 4.2.1), and difference analysis was conducted on the TCGA-BLCA dataset using the DESeq2 package (version 1.38.3). Enrichment analysis was performed using the clusterProfiler package (version 4.4.4), and Cox regression analysis was conducted using the survival package (version 3.4-0). ROC validation was conducted using the timeROC package (version 0.4), and nomograms related to independent prognostic factors were constructed using the survival and rms packages (version 6.3-0). GSVA analysis was performed using the GSVA package (version 1.44.5), expression analysis was conducted using the Wilcoxon test, and statistical significance was evaluated by setting a p-value threshold lower than 0.05.

A total of 5,686 BLCA-DEGs were identified in the TCGA-BLCA dataset, with 3,266 upregulated and 2,420 downregulated genes (Figures 1A,B). The crossing of 5,686 BLCA-DEGs, 1,548 PCD-RGs, and 1,136 MRGs resulted in the identification of 41 intersection genes (Figure 1C). These intersection genes might be involved in both the regulation of programmed cell death and the remodeling of metabolic abnormalities in bladder cancer cells, providing potential candidate targets for subsequent exploration of tumor progression mechanisms and targeted intervention strategies.

In GO functions, these intersection genes were enriched in the intrinsic apoptotic signaling pathway, apoptotic mitochondrial changes, release of cytochrome c from mitochondria, and so on (Figure 1D). Additionally, they were associated with KEGG pathways, such as apoptosis, platinum drug resistance, and apoptosis—multiple species (Figure 1E). Notably, intersection genes were enriched in mitochondrial and PCD-related functions and pathways, which was consistent with our screening expectations. After excluding 8 genes without any interactions, we built a PPI network comprising 33 nodes and 147 pairs of interactions, such as VDAC2 with BAX and CASP9 with VDAC1 (Figure 1F). The results of these enrichment analyses, together with the construction of the PPI network, revealed the potential molecular network underlying the regulation of programmed cell death and metabolic abnormalities in bladder cancer, providing clues for in-depth analysis of BLCA progression.

After univariate Cox analysis and PH assumption testing of 41 intersection genes, a total of 6 prognostic-related genes were retained (Figure 2A). These six genes underwent further screening through stepwise regression to identify five prognostic genes (POLB, FASN, CASP9, VDAC2, and RHOT2), and a multivariate Cox analysis was then employed to construct a prognostic model (Figure 2B). The prognostic model was formulated as follows: riskScore = POLB × (−0.41) + FASN × 0.35 + CASP9 × (−0.36) + VDAC2 × 0.23 + RHOT2 × (−0.27). The risk scores of BLCA samples in TCGA-BLCA were calculated and divided into high- and low-risk categories according to their optimal threshold (−0.4). The high-risk category exhibited a higher number of deceased samples than the low-risk category, indicating that the risk score could provide a certain degree of prognostic evaluation for patients (Figure 2C). The expression of prognostic genes in categories is shown in Figure 2D. The KM curve revealed significant disparities in survival rates between the two risk categories, with the higher one exhibiting poor survival outcomes (Figure 2E). The model demonstrated a robust predictive ability for sample survival, as evidenced by the high AUC values of 0.66, 0.67, and 0.71 for the 1-, 2-, and 3-year survival in ROC curves, respectively (Figure 2F). The construction of this prognostic model provided a quantitative tool for the individualized risk assessment of bladder cancer patients. The five key genes involved might influence BLCA progression by regulating programmed cell death and metabolic reprogramming, offering potential molecular targets for the development of targeted therapeutic strategies.

The generalizability of the prognostic model was assessed by performing the same analysis in GES32894. The optimal threshold for all samples in this dataset was determined to be −4.1. The high-risk sample exhibited a higher mortality rate and significantly poorer survival, reaffirming the results observed in the TCGA-BLCA (Figures 3A–C). The AUC values of the 1- (0.60), 2- (061), and 3-year (0.61) survival in ROC curves exceeded 0.60, further confirming the reliability of the risk model (Figure 3D).

The distribution of risk scores significantly differed among the age, grade, and stage groups, with risk scores increasing in correlation with age and disease severity in grade and stage (Additional file 2). In TCGA-BLCA, univariate Cox analysis, PH hypothesis test, and multivariate Cox analysis were implemented to determine the risk score, stage, and age as independent prognostic factors to construct the nomogram (Figures 4A–C). In the nomogram, as the total points increase, the patient’s survival probability decreases. The slope of the calibration curve approached 1, and the net proceeds of the model in the DCA were superior to those of the individual factors, indicating that the nomogram exhibited enhanced predictive efficacy (Figures 4D–G). By constructing a nomogram incorporating risk score, age, stage, and grade, the independent roles of these factors in the prognosis of bladder cancer patients were revealed. This provided a potential prognostic assessment tool for clinical practice and held certain clinical application value.

In the Hallmark background, epithelial–mesenchymal transition (EMT), MTORC1 signaling, and inflammatory response pathways were enriched in the high-risk patients (Figure 5A). KEGG pathways such as ECM–receptor interactions, cytokine–cytokine receptor interactions, and spot-of-adhesion were enriched in the high-risk category, while ribosome was enriched in the low-risk category (Figure 5B). The involvement of cytokine–cytokine receptor interactions, epithelial–mesenchymal transition, and inflammatory responses has been demonstrated in the pathogenesis of BLCA. This result revealed the potential roles of EMT, MTORC1 signaling pathway, and inflammatory response in the pathogenesis of BLCA, further expanding our understanding of the biological characteristics of BLCA.

The two risk groups displayed remarkable differences in the levels of 16 immune cells, such as natural killer (NK) T cells, mast cells, γδT cells, activated CD4⁺ T cells, and CD56dim NK cells (Figure 6A). Interestingly, with the exception of CD56dim natural killer cells and monocytes, the levels of the remaining 14 immune cells were found to be higher in the high-risk patients, and the risk score was positively correlated with these cells (Figure 6B). This indicated that the prognosis of BLCA could be markedly influenced by the immune system. We also found that POLB, FASN, CASP9, and RHOT2 were negatively correlated with most of the differential immune cells, while VDAC2 was the opposite. Further analysis of immune checkpoint expression showed that PDCD1LG2, HAVCR2, LAG3, and LGALS9 were notably differentially expressed between the groups (Figure 6C). The analysis of the immunotherapy effect revealed that both TIDE and exclusion scores exhibited higher values in the high-risk patients, showing an increased trend following the increase of the risk score (Figures 6D–F). The findings implied a greater likelihood of immune evasion in patients at high risk. These results revealed the role of the immune system in the prognosis of bladder cancer and provided certain information for the subsequent optimization of immunotherapeutic strategies. The results may offer theoretical support for the formulation of individualized immunotherapy regimens for high-risk patients.

A total of 29 drugs (including pictilisib, dactolisib, and erlotinib) exhibited prominent differences between the two risk categories, with a correlation p-value less than 0.01 for the association between drug IC₅₀ value and risk score (Additional file 3). Among them, Figures 7A,B illustrate the results of pictilisib and dactolisib. Their low IC₅₀ values indicated their potent inhibitory effect on tumor cells even at small drug concentrations. Except for CASP9 and POLB, the expression trends of the other three prognostic genes (FASN, VDAC2, and RHOT2) in the BLCA and control groups of the TCGA-BLCA and GSE13507 datasets were consistent, and the intergroup differences were noteworthy (Figures 7C–H). These results revealed the potential application value of specific drugs in the treatment of bladder cancer and provided a reference basis for individualized medication strategies based on risk scores.