The data for BLCA patients was obtained from the TCGA and encompasses transcriptomic and clinical information. A total of 431 files (comprising 406 cases) were acquired, consisting of 412 tumor files and 19 normal files. To facilitate analysis, the TPM format was employed for the data. Owing to the incomplete clinical data, certain information was omitted from the clinical study. The 277 PANoptosis genes were discerned from the existing literature (29).

We computed the t-statistics, LogFC, and P value using the “eBays” function. The comparative limma analysis (version 4.3.3) identified 1.5-fold differently expressed transcripts (adj.P<0.05) in BLCA, indicating PANoptosis-related molecular signatures in the TCGA cohort (30). We intersected the two gene sets to produce a collection of PANoptosis-associated BLCA genes (BLCA-PANs) for subsequent investigation.

Enrichment analysis for the BLCA-PANs was performed using the “org.Hs.eg.db” and “clusterProfiler” R packages (version 4.3.3) (31). All P values were less than 0.05. Based on protein-protein interaction (PPI) analysis, we identified communications and many key genes within the BLCA-PANs (STRING: functional protein association networks (string-db.org)).

Unsupervised consensus clustering utilizing the K-means algorithm was executed with ConsensusClusterPlus to identify novel BLCA molecular subtypes based on characteristic gene profiles (32). The empirical cumulative distribution function (CDF) was utilized to ascertain the appropriate number of clusters (33). We subsequently evaluated prognostic variations using the R package “survival” and BLCA-PANs signatures in connection to clinical outcomes and immune infiltration.

We conducted univariate Cox regression analysis to identify genes with P-values of less than 0.05. Following data preprocessing, the raw data were randomly divided into training and testing sets (1:1 ratio) using the “caret” package. A prognostic model was developed using LASSO stepwise regression (34). Utilizing 10-fold cross-validation, the λ value associated with the smallest mean squared error and its standard error (SE) was identified as the stable solution, and regularization methods were employed to reduce the hazards of overfitting (35). This procedure discovered features with non-zero coefficients and produced coefficient path visualizations and cross-validation error curves. Hazard ratios (HR) and their 95% confidence intervals (CI) were derived using the model gene coefficients, with findings displayed in a forest plot format.

Compute the AUC (Area Under Curve) and the P-value for survival analysis. The threshold for the training set is established at P < 0.01, whereas the threshold for the test set is set at P < 0.05. The training set AUC exceeds 0.65, while the test set AUC surpasses 0.63 (36). Feature selection and model training are conducted solely on the training set, with the test set used only for final validation and verification. The formula for calculating the risk score is as follows: RiskScore=∑i=1n (βi×Expressioni). Risk stratification thresholds are established by the predetermined median risk score, categorizing patients into high-risk and low-risk groups for a comparative survival study, including overall survival (OS) and progression free survival (PFS). Graph the C-index, AUC curve, and decision curve analysis (DCA) to assess the correlation between the model risk score and clinical baseline variations (37). Develop a nomogram utilizing clinical parameters and generate the calibration curve. The “maftools” R package was used to evaluate tumor mutational burden (TMB) (38).

We conducted several analyses based on risk stratification, encompassing RNA stemness score (RNAss), immunological subtypes, Gene Set Enrichment Analysis (GSEA), and Single Sample Gene Set Enrichment Analysis (ssGSEA). RNAss is a score system derived from transcriptome data that evaluates stem cell characteristics, primarily utilized to analyze the stemness aspects of cells in tumor or other tissue samples. 1000 permutation tests determined levels of significance (P < 0.05, FDR<0.25) to guarantee robust statistical inference (39, 40). Comparisons between high-risk and low-risk groups revealed immunophenotypic difference across four dimensions: effector cells, signaling pathways, functional annotations, and the repertoire of Immune checkpoint (IC) molecules. Immune cell infiltration analysis primarily relies on the CIBERSORTx algorithm. Furthermore, we computed Tumor Microenvironment (TME) scores derived from the stromal score, immune score, ESTIMATE score, and tumor purity to evaluate variations in the tumor microenvironment (41). Additionally, the IMvigor 210 dataset from the immunotherapy cohort was utilized for relevant assessment, encompassing Tumor Immune Dysfunction and Exclusion (TIDE) and Microsatellite Instability (MSI), which can elucidate immune evasion and immunotherapy for high-risk and low-risk individuals (42). Ultimately, we conducted a drug susceptibility prediction study using the “oncoPredict” R package, which is grounded in the prognostic model.

We utilized four machine learning techniques to identify significant feature genes for the model. Boruta does a top-down feature relevance analysis by systematically comparing the significance of characteristics with that of shadow attributes generated through the random permutation of the original qualities (43). It assesses significance by its permuted equivalents and systematically removes extraneous aspects to stabilize the evaluation. Support Vector Machine Recursive Feature Elimination (SVM-RFE) was applied to a dataset subjected to 10-fold cross-validation, with the number of folds set at 10. This produced indices for the training and testing sets. Following the application of the SVM-RFE algorithm to each training fold, features were prioritized according to their average rank (44). The Random Forest model is trained utilizing the “randomForest” package, with a specification of 2000 trees. A graph illustrates the Out-of-Bag (OOB) error rate of the Random Forest in relation to the number of trees. The Random Forest model is reconfigured using the ideal tree count, and feature significance is assessed via the “importance” function (45). Additionally, we cross-referenced the model genes with genes exhibiting differential expression identified using multi-omic analysis of BLCA from our previous publication, which analyzed urinary specimens from five BLCA cases compared to five healthy donors (46) and supplemented by additional proteomics from Zhang et al. (47). Ultimately we identified PSMB5 as the primary gene of interest.

We performed extensive analyses on PSMB5, encompassing gene expression profiling, assessment of survival probability, evaluation of progression-free survival, and clinical correlation studies. Additionally, we performed an extensive analysis of the immune landscape and tumor mutational burden to clarify immune-related characteristics and investigate possible implications for immunotherapy.

To investigate the causal link between PSMB5 and BLCA, we used MR using Wald ratio methods (48). Exposure data comprised three single nucleotide polymorphisms (SNPs) from the eqtl-a-ENSG00000100804, filtered using a clumping window size of 10,000 Kb, R2<0.001, and F>10 (Supplementary Table 1). The outcome data was obtained from FinnGen (https://www.finngen.fi/en).

We employed single-cell RNA sequencing data from the GEO dataset GSE222315, which includes 9 BLCA cases and 4 surrounding normal tissue samples. Raw scRNA-seq data were converted into Seurat objects and underwent quality control according to specified thresholds to preserve high-quality cells: (1) Detection of 200–5,000 genes per cell; (2) Mitochondrial gene content not exceeding 15%; (3) Red blood cell gene expression rate surpassing 3%. Following normalization, batch effects were corrected using the Harmony integration method. Data was subjected to log-normalization and subsequently scaled using linear regression (49). Dimensionality reduction was performed using principal component analysis, followed by graph-based clustering via the “Find-clusters” algorithm (50). Visualization was conducted using UMAP, and the expression informed the annotation of various cell populations of classical marker genes (51). We examined the disparities in PSMB5 expression across different cell types and between the negative and positive groups. The correlation between PSMB5 expression and lymph node metastases was investigated in particular cell lines.

All cell lines were acquired from the Gansu Province Clinical Research Center for Urinary System Diseases. SV-HUC-1 urethral epithelial cells were cultivated in Ham’s F12K medium, while BLCA cell (T24, UMUC-3, J82, 5637) were sustained in RPMI-1640 (Shanghai Yuanpei Biotechnology). Both media included 10% fetal bovine serum (FBS) from PAN Biotech and 1% penicillin-streptomycin at a concentration of 100 U/mL-100 μg/mL from Solarbio. Standard incubation conditions of 37 °C, 5% CO₂, and humidity were maintained consistently.

The two small interfering RNAs (siRNAs) directed against PSMB5 were procured from Tsingke Biological, and the transfection reagent was sourced from Shanghai GenePharma Biotechnology (si1: 5’-CGAAAUGCUUCAUGGAACA-3’; si2: 5’-GGCAAUGUCGAAUCUAUGA-3’; si-NC: UUCUCCGAACGUGUCACGUTT). The efficacy of the knockdown was validated using western blot (WB) analysis at 48 hours post-transfection. Moreover, concurrent phenotypic experiments were conducted using the same procedure.

The whole coding sequence of human PSMB5 was inserted into the pLV3-CMV-3×FLAG-mCherry-Puro vector (Miaoling Bio, China). HEK293T cells were co-transfected with psPAX2 and pMD2.G vectors, and the viral supernatant was harvested to infect J82 cells. Following puromycin selection, stable cell lines exhibiting PSMB5 overexpression were established. Cell transfection was performed using Polybrene (Solaibao, China) according to the manufacturer’s instructions.

Total protein was extracted utilizing RIPA buffer (P0013B, Beyotime, China) augmented with protease inhibitors. Protein concentrations were measured via the Bicinchoninic Acid assay. After separation by SDS-PAGE electrophoresis, proteins were transferred to PVDF membranes. Membranes for immunoblotting were blocked with 6% non-fat dry milk and then treated with primary antibodies at 4 °C overnight. Protein bands were identified utilizing the Odyssey imaging system in conjunction with the appropriate secondary antibody (926-32211, Li-Cor, USA) for visualization. This work utilized the following antibodies: β-actin (cat#66009-1-Ig, Proteintech) and PSMB5 (cat#19178-1-AP, Proteintech).

The Cell Counting Kit-8 (CCK8) was utilized to evaluate the proliferation. In accordance with the guidelines, cells (2 × 10³/well) were inoculated in 100 µL of media using 96-well plates, with three replicate plates established for various time points. CCK-8 reagent (AbMole BioScience) was applied at 10 µL per well at intervals of 0 to 96 hours. Following a 2-hour incubation, the optical density at 450 nm was assessed via a BioTek plate reader.

For clonogenic tests, 6-well plates were inoculated with 1 × 10³ cells per well in 2 mL of medium. Following an 8–10 day cultivation at 37°C with 5% CO₂, colonies were fixed with 4% PFA (Biosharp #BL539A), stained with 0.1% crystal violet (Solarbio #G1063), and subsequently photographed and quantified.

Transfected cells (6×10⁵) attained confluence 48 hours after transfection. Monolayers were scraped with sterile 200 μL tips, rinsed with PBS, and subsequently treated with serum-free media. Migration was evaluated by photographing wounds at 0 and 24 hours using inverted microscopy, with closure rates measured using ImageJ.

BLCA cells (1×10⁵ in 200 μL of serum-free media) were inoculated into LABSELECT chambers (8 μm holes; #14342). The lower chambers had 600 μL of RPMI-1640 enriched with 20% FBS as a chemoattractant. After 24–48 hours of incubation at 37 °C with 5% CO₂, the transmigrated cells were subjected to methanol fixation (4%), crystal violet staining (0.1%; Solarbio #G1063), and subsequent microscopic counting.

Apoptosis was evaluated utilizing the Annexin V-FITC/PI kit (Multi Sciences #AP101) in accordance with the manufacturer’s specifications. Flow cytometric analysis (Beckman CytoFLEX S) was used to assess overall apoptosis by aggregating early and late apoptotic populations.

To investigate the therapeutic potential of PSMB5 as a target, we employed the Coremine medical ontology information retrieval tool (www.coremine.com/medical/) to delineate PSMB5. Additionally, to obtain the target protein result files, the structures of TCMs were retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/), while the structure of PSMB5 (PDB ID: 5l5w) was acquired from the PDB database (https://www.rcsb.org/). The requisite alterations to the receptor proteins, encompassing hydrogenation and charge equilibrium, were executed utilizing AutoDockTools 1.5.7 software. AutoDock Vina 1.1.2 was subsequently employed to mimic molecular docking between the pharmaceuticals and PSMB5 (52). The molecular docking results were visualized using PyMOL 3.1.5.1, focusing on high-affinity complexes.

Statistical analyses were conducted using R (v4.3.3) and GraphPad Prism (v9.0). Continuous variables were compared between groups using either Student’s t-test (parametric) or Wilcoxon rank-sum test (non-parametric), based on normality assessment. Categorical variables were evaluated with the χ² test or Fisher’s exact test, chosen based on anticipated cell frequencies. Survival outcomes were evaluated with Kaplan-Meier estimation and log-rank testing for group comparisons, augmented by multivariate Cox proportional hazards regression. All experimental techniques were conducted in three biological replicates, with data presented as mean ± standard deviation (SD). Statistical significance was determined at p < 0.05, with asterisk notation indicating non-significant (n.s.); *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001; ****p ≤ 0.0001. P values below 0.05 were considered statistically significant (53).

We found 4,968 differentially expressed genes in BLCA. These were compared with 277 PANoptosis-related genes sourced from the literature, resulting in the identification of 98 BLCA-PANs (Figures 1A, B). Gene Ontology (GO) analysis revealed the abundance and enrichment significance of BLCA-PANs across various levels (Figure 1C). The biological process (BP) exhibited significant enrichment in proteasome-mediated ubiquitin-dependent protein degradation (Figure 1D). The cellular component (CC) revealed that the principal enrichment functions of BLCA-PANs were endopeptidase, peptidase, and proteasome complexes (Figure 1E). The molecular function (MF) exhibited significant enrichment in DNA-binding transcription factor interactions and ubiquitin-related ligase interactions (Figure 1F). The Kyoto Encyclopedia of Genes and Genomes (KEGG) indicated that BLCA-PANs were predominantly abundant in the proteasome and apoptotic pathways (Figure 1G). Furthermore, we conducted a PPI analysis to demonstrate the interactions of BLCA-PANs and identified several key genes primarily associated with the proteasome subunit family (Figures 1H, I).

Consensus clustering analysis was performed on BLCA-PANs expression patterns to categorize patients into two distinct subtypes: Cluster A (n = 246) and Cluster B (n = 160) (Figures 2A, B). A heatmap was later generated to depict the differential expression of BLCA-PANs concerning molecular subtypes (Cluster A/B) and clinicopathological characteristics, including gender, age, and staging factors (T, N, M) (Figure 2C). OS analysis indicated that Cluster A demonstrated a markedly inferior overall survival probability relative to Cluster B (p = 0.033; Figure 2D). The study revealed that most ICs were considerably overexpressed in Cluster A, while only a select few showed elevated expression in Cluster B (p < 0.05; Figure 2E). Immune infiltration with ssGSEA indicated that Cluster A exhibited a statistically significant prevalence of γδ T cells, while Cluster B had a predominant infiltration of activated CD8 T cells and CD56bright natural killer cells (p < 0.05; Figure 2F). Waterfall charts of TMB indicated that Cluster A displayed elevated gene mutation rates compared to Cluster B (Figures 2G, H). Furthermore, the TMB score values indicated a significant difference between the two clusters (p < 0.05; Figures 2I-L).

Employing Cox regression studies, we developed a prediction model comprising four BLCA-PANs by LASSO regression (Figures 3A-D). The model exhibited enhanced predictive accuracy relative to the clinical baseline, as evidenced by C-index, AUC curve, and DCA analyses (Supplementary Figures 1A-C). The Receiver Operating Characteristic (ROC) curves demonstrated the model’s prognostic capability (Figure 3E). Based on the computed risk scores in BLCA (BLCA-Riskscore), patients were classified into High- and Low-risk categories (Figure 3F). Marked enhancements in OS (p < 0.001) and PFS (p = 0.007) were noted in the Low-risk group (Figures 3G, H). As a result, we constructed a nomogram (Figure 3I) and a standard curve (Supplementary Figure 1D).

To validate the model’s credibility and consistency, we partitioned the TCGA database into training and testing sets (Figures 4A, B). The operating system results for the two sets were consistent with the prior findings (p-value for test set = 0.038, p-value for train set = 0.001; (Figures 4D, E). Utilizing the BLCA-Riskscore to evaluate Clusters A and B, Cluster A had markedly higher risk scores compared to Cluster B, correlating with inferior overall survival rates in this cohort (Figure 4C). The Sankey diagram was used to illustrate the relationship between the groups and clinical features (Figure 4F). Furthermore, we conducted GSEA analyses for high-risk and low-risk groups based on GO and KEGG. The high-risk group was primarily characterized by the chemotaxis and migration of granulocytes and neutrophils, as well as the interaction with extracellular matrix receptors and the activation of the JAK-STAT signaling pathway. The low-risk group was linked to the epoxygenase P450 pathway, arachidonic acid epoxygenase or monooxygenase activity, and so on(Figures 4G-J).

The cascade charts illustrated the disparity in mutational landscapes between high-risk and low-risk groups (Figures 5A, B). The association investigation indicated a small inverse correlation between RNAss and risk score, implying diminished stemness characteristics (Figure 5C). We identified statistically significant dysregulations in pathways, notably impacting the KRAS cascade, NF-κB-mediated TNF-α signaling, β-catenin-dependent WNT pathway, TGF-β transduction, IL-6-JAK-STAT3 axis, and the PI3K-AKT-mTOR network (Figures 5D, E).

According to ssGSEA, plasmacytoid dendritic cells (pDC) exhibited substantial immunological infiltration in the high-risk cohort. Moreover, CD8⁺ T cells and CD56^brilliant natural killer (NK) cells were significantly infiltrated in the low-risk cohort (Figure 5F). The high-risk cohort had significant expression of M0, M1, and M2 macrophages, corresponding with specific immunological activity patterns. Conversely, the low-risk group exhibited a predominance of immunosuppressive regulatory T cells (Tregs), antibody-secreting plasma cells, monocytic lineage cells, and activated dendritic cell populations (Figure 5G). The examination of immune function indicated that immune responses were predominantly heightened in the high-risk group, encompassing APC co-inhibition, APC co-stimulation, MHC-I, neutrophils, para-inflammation, T cell co-inhibition, T cell co-stimulation, Th1 cells, Th2 cells, tumor-infiltrating lymphocytes (TILs), Tregs, and Type-I interferon response. Only the Type-II IFN response is considerably elevated in the low-risk group (Figure 5H). A cohort of 371 BLCA patients was categorized into four clusters (C1, C2, C3, C4) and classified as high-risk or low-risk based on established BLCA-Riskscore thresholds. A statistically significant difference was noted between risk strata using the chi-square test (Figure 5I). The comparative study of IC expression profiles differentiated the high-risk group from the low-risk group. Most ICs exhibited considerable differential expression among cohorts (Figure 5J).

Analysis of immune cell correlations indicated that seven immune cell types were significantly associated with the BLCA-Riskscore, comprising three positively correlated (M0, M2, Neutrophils) and four negatively correlated [Dendritic cells (activated), Monocytes, T follicular helper cells (Tfh), Tregs] cell types (Figures 6A-J). The TME score exhibited notable disparities (Figure 6K). Furthermore, assessments of TIDE and MSI indicated that the TIDE score, MSI status, and Dysfunction score, excluding the Exclusion score, exhibited considerable variance specific to the cohort (Figure 6L). Ultimately, drug susceptibility analysis revealed that 29 medications exhibited a significant correlation (Supplementary Figure 2).

To advance research on BLCA-PANs, we employed four machine learning techniques and integrated the results of Zhang et al. with our urine proteomics data to identify critical feature genes (Figure 7A). The Boruta algorithm demonstrates that PSMB5 attained the highest score (Supplementary Figures 3A, B). SVM-RFE indicates that PSMB5 is the nearest to the scatter point, exhibiting the highest average ranking (Supplementary Figure 3C). PSMB5 demonstrates the most excellent absolute coefficient value in the Lasso regression model (Supplementary Figure 3D). Random Forest demonstrates the most excellent Mean Decrease Gini score (Supplementary Figures 3E, F). The essential gene PSMB5 was ultimately acquired.

Differential expression analysis revealed that PSMB5 was markedly overexpressed in BLCA patients (Figure 7B). The studies of OS (p = 0.025) and PFS (p < 0.001) demonstrated that PSMB5 was significantly associated with clinical prognosis (Figures 7C, D). A heatmap illustrating the relationship between PSMB5 and clinical characteristics was generated (Figure 7E). Subsequently, we examined the connection between PSMB5 and other BLCA genes (Figure 7F). Analysis of immune cells yielded results consistent with the high-risk and low-risk groups, correlating with M0, M2, Monocytes, and Tregs (Figure 7G). Furthermore, we identified several ICs that were statistically significant with PSMB5 (Figure 7H). The TMB exhibited a positive correlation with PSMB5 expression (p = 0.036; Figure 7I).

Furthermore, we identified three SNPs in PSMB5 (rs12590429, rs117058979, rs11543947) to conduct MR. Results identified rs117058979 as a causative variant for BLCA [OR = 2.267 (1.008, 5.097), p = 0.048] (Supplementary Table 2).

Upon normalizing and annotating the single-cell database, the BLCA group and the normal group predominantly clustered into nine categories of cell lines (Figure 8A). PSMB5 exhibited markedly elevated expression in the BLCA cohort, predominantly among Endothelial cells, Epithelial cells, and Fibroblasts (Figures 8B, C). Consequently, we meticulously analyzed the relationship between PSMB5 expression and lymph node metastasis in all three cell lines, discovering substantial statistical differences for endothelial cells (p < 0.0001) and fibroblasts (p < 0.0001) (Figure 8D).

WB analysis revealed distinct expression profiles of PSMB5 across various BLCA cell lines, with significantly increased expression levels observed in T24 and UMUC-3 cells (Figure 9A). In these two cell lines, siRNA transfection resulted in a knockdown efficiency of about 60% for PSMB5. Subsequently, comprehensive in vitro functional studies were performed. (Figure 9B). The CCK-8 proliferation assay and analysis of colony formation consistently indicated that PSMB5 depletion markedly reduced cellular proliferation compared to the negative control (NC) groups (Figures 9C, D). Additionally, both wound-healing and transwell migration experiments demonstrated significantly reduced migratory ability in PSMB5-knockdown cells compared to controls (Figures 9E, F). Flow cytometric examination of apoptosis revealed that silencing PSMB5 markedly increased apoptotic rates compared to the NC groups (Figure 8G). After overexpressing PSMB5 in J82 with an overexpression efficiency of about 40%, the opposite biological behavior was displayed (Supplementary Figure 4A). The abilities of proliferation and migration are enhanced, and cell apoptosis is significantly reduced (Supplementary Figures 4B-F).

We identified five TCMs related to PSMB5 from the Coremine dataset: Chuanxiong Rhizoma (Chuan Xiong in Chinese), Ligusticum sinense Oliv. (Gao Ben in Chinese), Fuxiong Rhizome (Fu Xiong in Chinese), Tripterygium wilfordii Hook. f. (Lei Gong Teng in Chinese), and Scutellaria baicalensis Georgi (Huang Qin in Chinese) (Figure 10A). Subsequently, molecular dockings were performed, revealing that binding energies below -5.0 kcal·mol^-¹ indicated increased molecular affinity (Figures 10B-F). These interactions offer a potential pathway for further investigation into the use of TCMs in the treatment of BLCA.