Gene expression data for patients from the TCGA-BLCA project, including 41 normal and 487 tumor samples, were obtained using the TCGAbiolinks R package. A query was then constructed to retrieve Gene Expression Quantification data processed with the STAR - Counts workflow. FPKM-normalized expression values were obtained from the fpkm_unstranded assay provided within the GDC-generated SummarizedExperiment (Supplementary Data S1). Ensembl gene identifiers were converted to official gene symbols using a custom Perl script and a reference human GTF annotation file from the GENCODE project (https://www.gencodegenes.org/human/) (Supplementary Data S2). Duplicate gene symbols were resolved by retaining only the entry with the highest overall expression across all samples, resulting in a non-redundant, gene-level expression matrix suitable for downstream analyses.

For external validation, three independent bladder cancer cohorts were obtained from the GEO database: GSE13507 (GPL6102, n = 166), GSE32894 (GPL6947, n = 224), and GSE48276 (GPL14951, n = 73). Probe IDs from each platform were mapped to gene symbols using corresponding annotation files. Each dataset was independently processed using the limma package, with log₂ transformation applied where appropriate and quantile normalization was performed. Genes common across all datasets were identified, and batch effects were corrected using the ComBat() function from the sva package. The resulting normalized and integrated expression matrix was used for downstream validation analyses.

To establish a comprehensive set of PRGs, we integrated data from three well-recognized and previously validated sources: PeroxisomeDB 2.0 (comprising 100 genes), the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (83 genes), and a proteomic dataset specific to human liver peroxisomes (60 genes) (14–16). After consolidating these datasets and eliminating redundant entries, we curated a final list of 113 non-overlapping PRGs, as detailed in Supplementary Table S1.

To enable integrative analysis, genes common to both TCGA and GEO cohorts were identified. Expression matrices were merged, and batch effects between TCGA and GEO samples were corrected using the ComBat function from the sva package (Supplementary Data S3). Shared PRGs were extracted from the corrected matrices and saved for downstream prognostic modeling. To develop the PRG-based prognostic signature, univariate Cox regression was used to identify survival-associated genes, followed by multivariate Cox regression to select independent prognostic markers. Risk scores were calculated as the weighted sum of gene expression and regression coefficients. Patients were stratified into high- and low-risk groups based on the median risk score.

To identify prognostic biomarkers, univariate Cox proportional hazards regression was first performed to screen survival-associated genes (p < 0.05). Candidate genes were then further evaluated using a random survival forest model implemented with the randomForestSRC and randomSurvivalForest packages to assess prognostic importance. Genes with a relative importance score greater than 0.5 were considered prognostically relevant and retained for subsequent analyses. For diagnostic gene selection, a random forest classifier was applied, and feature importance was ranked using the Mean Decrease in Gini index, with top-ranking genes selected for downstream diagnostic evaluation.

Somatic mutation data for bladder cancer (TCGA-BLCA) were downloaded and processed using the TCGAbiolinks package. Mutation Annotation Format (MAF) files were curated and integrated with clinical risk information to stratify samples into high- and low-risk groups. Mutational profiles were analyzed and compared between groups using the maftools package. Mutations were further filtered to highlight selected signaling or metabolic pathways. Gene-specific mutation patterns for TP53 were visualized using lollipopPlot(), and co-occurrence or mutual exclusivity among frequently mutated genes was assessed using somaticInteractions() and coOncoplot(). Tumor mutation burden (TMB) was calculated using the tmb() function.

To assess the enrichment of hallmark cancer pathways (KEGG Hallmarks) and immune-related pathways (curated gene lists; Supplementary Table S2) across PRG subtypes, single-sample Gene Set Enrichment Analysis (ssGSEA) was performed using the GSVA package in R. Normalized gene expression data, including log2 transformation and batch correction where applicable, were used as input to compute ssGSEA enrichment scores with the gsva function (method = “ssgsea”, kcdf = “Gaussian”, abs.ranking = TRUE) (Supplementary Data S4). The resulting ssGSEA scores for each gene set were scaled to the 0–1 range using min–max normalization performed per gene set. Differences in ssGSEA scores between groups were evaluated using the Wilcoxon rank-sum test, and the normalized ssGSEA matrix was used for downstream statistical analysis and visualization. To evaluate Gene Ontology (GO) and KEGG pathway enrichment between PRG subtypes, we utilized the clusterProfiler package in R, along with org.Hs.eg.db, by converting differentially expressed genes into Entrez IDs and applying enrichKEGG with p- and q-value thresholds set at 0.05. (Supplementary Data S5).

To comprehensively characterize the tumor microenvironment (TME), we first utilized the ESTIMATE algorithm via the “ESTIMATE” R package to infer tumor purity and quantify stromal and immune cell infiltration based on gene expression data (17). To estimate the absolute abundance of tumor microenvironment cell populations, we applied the MCP-counter method, which quantifies ten major cell types, including eight immune and two stromal populations, using predefined transcriptomic marker signatures. These signatures correspond to a fixed set of 112 curated marker genes implemented as the default gene sets in the MCP-counter R package (Supplementary Data S6; Supplementary Table S3) (18). These complementary approaches enabled a robust and multidimensional assessment of the TME landscape.

Drug–gene interaction enrichment analysis was performed using the Drug Signature Database (DSigDB) (https://dsigdb.tanlab.org/DSigDBv1.0/). Drug–gene association data were analyzed using the clusterProfiler package with gene annotation from org.Hs.eg.db to identify compounds significantly enriched for high-risk PRGs. Statistically significant candidates were selected based on p < 0.05 and FDR < 0.05 thresholds. The full list of enriched compounds was retained for transparency, while the final candidate was selected based on (i) biological plausibility, (ii) drug-likeness and safety, and (iii) structural tractability for docking analyses.

To evaluate potential interactions between the lead compound and high-risk PRGs, molecular docking simulations were performed against three target proteins: ACOX2, IDI1, and PRDX1. The 3D structure of erythorbic acid was retrieved from the PubChem database (CID: 54675810) in SDF format. The 3D crystal structures of IDI1 (PDB ID: 2ICK) and PRDX1 (PDB ID: 7WET) were downloaded from the RCSB Protein Data Bank (PDB), while the structure of ACOX2 was obtained from the AlphaFold Protein Structure Database (AF-Q99424-F1-model_v6.pdb) due to the absence of an experimentally resolved model. All receptor structures were pre-processed by removing water molecules, co-crystallized ligands, and heteroatoms. Docking simulations were conducted using CB-Dock2 (https://cadd.labshare.cn/cb-dock2/index.php), which automatically detects binding cavities and ranks docking poses based on predicted binding affinity (Vina score). The top-ranked pocket for each protein was selected according to the lowest binding energy and the presence of catalytically relevant residues.

This study employed the human bladder cancer cell line T24 and the normal bladder epithelial cell line XVHUC1, both sourced from the Committee of Type Culture Collection, Chinese Academy of Sciences (Shanghai, China). Cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum (FBS), 100 U/ml penicillin, and 100 μg/ml streptomycin. All cultures were maintained in a humidified incubator at 37 °C with 5% CO₂. Cell line authenticity was routinely confirmed through morphological assessment and regular mycoplasma testing to ensure contamination-free conditions.

Total RNA was extracted and purified using TRIzol Reagent (Takara, Otsu, Japan), followed by reverse transcription to generate cDNA. Quantitative real-time PCR (qRT-PCR) was carried out using the SYBR Green PCR Kit (Takara, Otsu, Japan). GAPDH was used as an internal reference to normalize mRNA expression levels, and relative expression in the treated samples was compared to that of the control group. Primers were designed using Primer Premier 5.0 and Beacon Designer 7.8 software and synthesized by General Biosystems Co., Ltd. The primer sequences are provided in Table 1.

Cell viability was evaluated using the Cell Counting Kit-8 (CCK-8) assay in the human bladder cancer cell lines T24 and HT-1197, as well as the normal urothelial cell line XVHUC-1. Cells were harvested by washing with phosphate-buffered saline (PBS) followed by brief trypsinization (1–2 min). After neutralization with complete culture medium, cells were gently resuspended to obtain a single-cell suspension and counted using a hemocytometer. The cell density was adjusted to 5 × 10⁴ cells/mL, and 100 μL of the suspension was seeded into each well of a 96-well plate. After 24 h of incubation, the medium was replaced with 100 μL of fresh medium containing the drug at indicated concentrations (0, 0.4, 0.8, 1.2, and 1.6 mM). Each concentration was tested in triplicate wells. Cells were then incubated under standard culture conditions, and cell viability was assessed by measuring optical density (OD) values at 24 h and 48 h following drug treatment according to the manufacturer’s instructions. Absorbance values were normalized to untreated control wells to calculate percent cell viability.

Categorical variables were compared using the chi-square test, and continuous variables using the Student’s t-test or Wilcoxon rank-sum test, as appropriate. Kaplan-Meier survival curves were evaluated with the log-rank test, and correlations were assessed using Pearson or Spearman methods. A p-value < 0.05 was considered statistically significant. All statistical analyses were performed using R (v4.0.3).

The overall analytical workflow and key findings of this study are summarized in Figure 1. Using peroxisome-related genes (PRGs), bladder cancer patients were stratified into distinct molecular subtypes with different prognostic outcomes. These subtypes displayed marked differences in clinical features, mutation patterns, pathway activities, and tumor microenvironment characteristics, leading to the identification of representative PRG biomarkers and a candidate therapeutic compound with experimental validation.

To explore the prognostic role of PRGs in bladder cancer, we applied a two-step feature selection strategy involving univariate Cox regression followed by multivariate Cox regression. Transcriptomic data from the TCGA-BLCA cohort served as the training dataset (n=401), while a merged expression matrix compiled from three GEO bladder cancer cohorts (n=463) was used for independent validation. After integrating the datasets and retaining only overlapping genes, a total of 91 shared PRGs were identified for analysis. Univariate Cox regression analysis on the TCGA-BLCA cohort revealed 12 PRGs significantly associated with patient prognosis (Figure 2A). These were subsequently refined to 6 core prognostic genes using multivariate Cox regression (Figure 2B). Among them, four genes including Peroxiredoxin 1 (PRDX1), Isopentenyl-Diphosphate Delta Isomerase 1 (IDI1), Acyl-CoA Oxidase 2 (ACOX2), and Dehydrogenase/Reductase 4 (DHRS4) were classified as high-risk genes, as their elevated expression was associated with poor prognosis. In contrast, Acyl-CoA Synthetase Long Chain Family Member 5 (ACSL5) and Xanthine Dehydrogenase (XDH) were identified as protective, with higher expression linked to better survival in bladder cancer. Based on the expression levels and corresponding regression coefficients of these six genes, a risk score was calculated for each patient. Using the median risk score as a cutoff, patients in both the training and validation cohorts were stratified into high-risk and low-risk subgroups. This stratification revealed significant prognostic differences between the two subgroups in each dataset, supporting the ability of PRG-derived risk scores to define novel molecular subtypes of bladder cancer (Figures 2C, D). The expression profiles of the six PRGs across both datasets are shown in Figures 2E, F. Principal component analysis (PCA) confirmed the subgroup distinction, with the first principal component (PC1) accounting for 26.3% of the variance in the TCGA BLCA cohort and 22.3% in the validation dataset (Figures 2G, H).

We next assessed the clinical characteristics associated with the identified PRG subtypes. No significant differences in age or gender distribution were observed among the PRG-defined subtypes (Figure 3A). However, marked differences emerged with respect to tumor stage, primary tumor size and extent, and lymph node involvement. The high-risk subgroup showed a greater proportion of patients with advanced-stage disease (stage III/IV) and extensive nodal metastasis. Additionally, tumors classified as T1-T3 were more commonly observed in the high-risk group, while early-stage tumors (Ta, stage 0) were more prevalent in the low-risk subgroup. Similarly, patients with one to three positive lymph nodes were more frequent in the high-risk category (Figure 3A). Notably, the PRG-based risk model retained its prognostic value across different clinical strata, including age, gender, and stage groups (0-II/III-IV), demonstrating its potential as an independent and robust classifier (Figures 3B–D). In conclusion, PRG subtyping effectively stratifies bladder cancer patients by tumor aggressiveness and provides meaningful prognostic insights beyond conventional clinical parameters.

We conducted an in-depth analysis of genomic aberrations across the PRG subtypes to uncover key mutational patterns. The oncoplot visualizations highlight distinct mutation distributions between PRG subtypes (Figures 4A, B). The oncoplot of mutational pathways revealed three key altered pathways distinguishing bladder cancer risk groups. The genome integrity pathway, characterized by mutations in TP53 and STAG2, and the Chromatin SWI/SNF complex pathway, defined by ARID1A mutations, were observed across cohorts. Additionally, a transcription factor pathway driven by ELF3 and NFE2L2 mutations was exclusive to the high-risk subgroup, whereas the low-risk subgroup exhibited an “Other” pathway defined by mutations in SPTAN1, APOB, and CACNA1A. The lollipop plot illustrated a higher mutation rate in the high-risk subgroup (56.72%, max height 13) compared to the low-risk subgroup (41.09%, max height 5), with mutation hotspots clustering within key functional domains of TP53 (Figures 4C, D). Mutation burden analysis demonstrated that, despite the higher mutation rate in the high-risk group, the low-risk subgroup showed a slightly greater overall mutation burden (3.17 mut/Mb vs. 2.64 mut/Mb) (Figures 4E, F). This suggests that mutations in the low-risk group tend to affect larger genomic regions despite occurring in fewer patients. Somatic interaction analysis further revealed distinct co-occurrence patterns (Figures 4G, H). TP53 mutations frequently co-occurred with RB1, while TTN mutations predominantly co-occurred with MUC16. Notably, the low-risk subgroup exhibited more co-occurring partner genes for both TP53 and TTN than the high-risk subgroup, indicating a higher frequency and broader mutational diversity within low-risk tumors. Together, these results delineate distinct mutational landscapes and genomic interactions between bladder cancer risk groups, providing insights that may inform tailored therapeutic approaches.

The identified PRG subtypes exhibited distinct clinical, molecular, and mutational characteristics. To further explore the underlying biological mechanisms differentiating these subtypes, we conducted a transcriptomic analysis to uncover key functional pathways. Pathway enrichment analysis revealed broad activation of hallmark oncogenic and stress-response pathways in the high-risk subgroup (Figure 5A). These included key cell cycle regulators (E2F targets, G2M checkpoint, MYC targets V1, mitotic spindle), DNA damage repair, hypoxia, EMT, angiogenesis, apoptosis, and immune-related signatures such as inflammatory response, IL6-JAK-STAT3 signaling, and allograft rejection, suggesting a more aggressive tumor phenotype. In contrast, the peroxisome pathway was among the few significantly upregulated in the low-risk subgroup, indicating a more regulated metabolic state (Figures 5B, C). As ACSL5 and XDH were the only PRGs associated with a protective prognostic effect, and ACSL5 is a key enzyme in long-chain fatty acid metabolism, ACSL5 may serve as a peroxisome-related biomarker of favorable prognosis in bladder cancer. Conversely, PRDX1, a peroxisomal antioxidant enzyme strongly expressed in the high-risk subgroup, may reflect heightened oxidative stress and survival signaling. Together, these results highlight distinct metabolic and stress-adaptive signaling between risk subgroups, with peroxisome regulation—particularly via ACSL5 and PRDX1—contributing to the biological divergence and prognostic differences observed.

Separately, several immune-related pathways showed distinct enrichment patterns between subgroups (Figure 5D). The low-risk subgroup was characterized by enrichment of the Type II interferon response, whereas the high-risk subgroup exhibited enrichment of HLA, CCR signaling, inflammation-promoting responses, T cell co-inhibition and checkpoint pathways, and cytolytic activity. These differences highlight subgroup-specific variations in immune regulation, surveillance, and inflammatory signaling.

GO and KEGG enrichment analyses based on differentially expressed genes (DEGs = 27; logFC ≥ 0.5, adjusted p-value < 0.05) further supported the biological divergence between subgroups (Figures 5E–G). GO terms highlighted enrichment in receptor–ligand interactions, extracellular matrix organization, and serine-type peptidase activity, suggesting differences in cell adhesion, matrix remodeling, and intercellular communication. Terms related to muscle structure and immune cell migration also indicated subgroup-specific variation in tissue and immune processes (Figure 5F). KEGG analysis identified significant enrichment in PPAR signaling, a key regulator of lipid metabolism and inflammation, along with ECM–receptor interaction, cytokine signaling, and steroid hormone biosynthesis. These findings underscore subgroup differences in metabolic regulation, immune signaling, and hormonal activity, aligning with observed transcriptional and functional divergence (Figure 5G).

To further understand the tumor microenvironment (TME) differences between the two PRG subtypes, we evaluated the immune and stromal components using the ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) algorithm, which infers stromal and immune cell presence based on gene expression data. This analysis revealed a significantly enriched tumor microenvironment in the high-risk subgroup, characterized by higher immune- and stromal-related signature scores, as reflected by elevated ImmuneScore and StromalScore values derived from the ESTIMATE algorithm (Figures 6A–C).

To precisely characterize cellular composition, we applied the MCP-counter (Microenvironment Cell Populations counter) algorithm, which identified a notable increase in monocytic lineage cells, natural killer (NK) cells, dendritic cells of monocytic lineage, cytotoxic T cells, CD8+ T cells, and fibroblasts within the high-risk subgroup (Figure 6D). At the gene level, ACOX2 and PRDX1 exhibited positive correlations with immune and stromal cell abundance scores, whereas ACSL5 showed an overall negative correlation trend (Figure 5E). Notably, the most pronounced correlations were observed with monocytic lineage cells and fibroblasts, while associations with other immune cell types were comparatively modest. In addition, the PRG risk score showed positive correlations with the enrichment levels of multiple immune and stromal cell populations, particularly fibroblasts and monocytic lineage cells, as illustrated for representative cell types in Figures 6F–H.

To identify robust biomarkers that distinguish PRG-based subtypes in BLCA, we employed a random forest classifier. Among the diagnostic feature genes, PRDX1 and ACSL5 emerged as the top two candidates based on the MeanDecreaseGini index, indicating their strong contribution to subtype classification (Figures 7A, B). To evaluate their prognostic relevance, we conducted univariate Cox regression and random survival forest analyses. Notably, PRDX1 showed the highest variable importance score in the survival forest analysis, followed closely by IDI1 and ACSL5 (Figures 7C–E), underscoring their association with overall survival within PRG-defined risk groups. Next, we stratified patients according to the median expression levels of PRDX1 and ACSL5. Heatmap visualization revealed strong concordance between gene-based classification and established risk subgroups, as well as key clinical features (Figure 7F). PCA further supported this stratification, with the first PC1 explaining 55.9% of the total variance between subgroups (Figures 7G, H). Diagnostic performance was assessed using ROC analysis, where ACSL5 achieved an AUC of 0.799 and PRDX1 an AUC of 0.783, highlighting their effectiveness in differentiating PRG-based subtypes (Figure 7I). Moreover, Kaplan–Meier survival analysis revealed a significant difference in overall survival between high and low ACSL5 expression groups (p = 0.013), while PRDX1 showed even stronger prognostic value, with an inverse survival trend and highly significant difference (p < 0.001) (Figures 7J, K). Collectively, these results position PRDX1 and ACSL5 as representative diagnostic and prognostic biomarkers for PRG-based molecular stratification in BLCA, with potential utility in improving patient risk assessment and personalized treatment strategies.

We experimentally assessed the mRNA expression levels of candidate PRG risk genes in the bladder cancer cell line T24, compared to the normal bladder epithelial cell line XVHUC1. The analysis revealed that high risk PRG genes were significantly upregulated in T24 cells relative to XVHUC1 cells. ACSL5 showed significant downregulation while XDH did not show any significant change in expression (Figure 8).

Drug-gene interaction enrichment analysis using the DSigDB database identified several compounds significantly associated with the three high-risk peroxisome-related genes (ACOX2, IDI1, and PRDX1) (Figure 9A). Among the top statistically significant hits, 1-chloro-2,4-dinitrobenzene exhibited the lowest p-value (9.29 × 10^-6) and interacted with all three target genes, followed by DL-methionine, acrylamide, erythorbic acid, hydrogen, and L-cysteine, each associated with subsets of ACOX2 and PRDX1 (Figure 9B). However, given the known toxic or nonspecific reactivity of compounds such as 1-chloro-2,4-dinitrobenzene and acrylamide, erythorbic acid was prioritized as the most biologically relevant and pharmacologically safe candidate for subsequent molecular docking. Erythorbic acid demonstrated significant enrichment (p = 5.3 × 10^-4; FDR = 0.0107) and shared strong predicted associations with both ACOX2 and PRDX1, two enzymes upregulated in the high-risk molecular subgroup of bladder cancer. These results supported the selection of erythorbic acid for in-depth docking and structural interaction analysis.

To investigate the binding interaction of erythorbic acid with PRGs, we first retrieved the 3D chemical structure of erythorbic acid from the PubChem database as shown in Figure 9C. Protein 3D structures for ACOX2, IDI1 and PRDX1 were also obtained and were then used as input for docking analysis via CB-Dock2, a cavity-detection-based blind docking platform (Figures 9D–F). Molecular docking of erythorbic acid with ACOX2, IDI1, and PRDX1 revealed stable and energetically favorable binding within functionally important regions of all three proteins (Table 2; Figures 9D–L).

For ACOX2, erythorbic acid occupied the enzyme’s largest FAD-adjacent pocket (volume = 13,116 ų) with a binding energy of –6.2 kcal/mol (Table 2; Figures 7D–F). Key interacting residues included ARG50, TYR105, GLU437, and ASP628, consistent with the catalytic environment responsible for β-oxidation. This interaction suggests that erythorbic acid may modulate peroxisomal lipid oxidation and ROS production. In IDI1, erythorbic acid exhibited the strongest affinity (–6.6 kcal/mol) within the Mn²⁺-coordinated catalytic pocket, interacting with HIS41, HIS52, CYS87, and TYR137 (Table 2; Figures 9G–I). This site corresponds to the isomerase’s active domain involved in the mevalonate pathway, implicating erythorbic acid as a potential modulator of isoprenoid biosynthesis and redox coupling in tumor metabolism. For PRDX1, the compound bound stably near the peroxidatic cysteine region (Pocket 2; volume = 249 ų; –5.4 kcal/mol), contacting residues LYS37, TYR38, LEU69, ASN70, ASP134, and PHE165 adjacent to Cys52, the active thiol responsible for peroxide reduction (Table 2; Figures 7J-L). This interaction indicates potential interference with PRDX1’s antioxidant function, which may reduce the redox resilience of high-risk tumor cells.

Overall, erythorbic acid demonstrated multi-target binding across metabolic and antioxidant PRG hubs, supporting its potential as a broad-spectrum peroxisomal modulator for high-risk bladder cancer subtypes.

To evaluate the effect of erythorbic acid on bladder cancer cell viability, CCK-8 assays were performed in the bladder cancer cell lines T24 and HT-1197, with XV-HUC-1 normal urothelial cells included as a non-malignant comparator. Cells were treated with increasing concentrations of erythorbic acid and assessed at 24 h and 48 h. Erythorbic acid induced a dose- and time-dependent reduction in cell viability in both T24 and HT-1197 cells relative to untreated controls, with progressively greater decreases observed at higher concentrations and longer exposure durations. In contrast, XV-HUC-1 cells exhibited comparatively limited sensitivity, maintaining near-baseline viability at 0.4–0.8 mM, while higher concentrations resulted in only a moderate reduction, with viability remaining at or above approximately 80%, a level commonly regarded as indicative of minimal cytotoxic effects in in vitro assays of normal epithelial cells (19) (Figure 10).