An Integrated Strategy for Predicting Cancer Aggressiveness and Therapeutic Vulnerability in Muscle Invasive Bladder Cancer

1. Introduction

According to the GLOBOCAN Cancer Statistics 2020, bladder cancer ranks as the tenth most common cancer worldwide, with respective incidence and mortality rates of 11.9 and 4.2 per 100,000 persons [1]. Men are approximately 4 times more commonly diagnosed than women, and thus bladder cancer is the 9th leading cause of cancer death among men. Tobacco smoking, occupational exposures in the aluminum, painting or rubber industries and drinking of arsenic contamination water are major risk factors of bladder cancer in several populations [2]. Histologically, the majority of bladder cancer is characterized as urothelial carcinoma (95% incidence), which could be further classified into non-muscle invasive (NMIBC) and muscle invasive (MIBC) bladder cancer. To determine the treatment and predict the outcomes, the tumor aggressiveness in a patient is usually staged and graded according to the extent and spread of the cancer cells by using transurethral resection of bladder tumor (TURBT) and radiologic imaging such as CT and MRI.

For patients with NMIBC, the lesions can be easily ‘‘shaved off’’ by TURBT. The follow-up treatments like adjuvant intravesical chemotherapy or Bacillus-Calmette-Guerin (BCG) instillation have much improved the clinical outcomes [3], and the mortality rates have been declining in most developed settings. Nevertheless, such treatment is insufficient for cases with MIBC due to the limitation that TURBT may not screen the whole organ, leading to a false diagnosis. Most patients at this stage eventually experience recurrence easily after receiving curative-intent treatment and start to develop metastatic disease [4]. The 5-year overall survival (OS) rate for patients with metastatic tumor is only 6%, in comparison with 90% in patients with NMIBC [1,5]. One standard treatment for MIBC is radical cystectomy to remove part of or the whole bladder. For patients in selected situations, chemotherapy may be applied. Although pivotal trials published decades ago have set first-line platinum-based combinations as the standard of care for patients at advanced stages, the median OS still remains 12-16 months, even treated with a second-line chemotherapy [6]. In recent years, immune checkpoint inhibitors (ICIs) in the neoadjuvant setting have emerged as another antitumor strategy for treating MIBC with some encouraging outcomes reported [7,8,9]. Nevertheless, only around 15% to 25% of patients with advanced BCs showed response to ICI treatment [10,11]. Therefore, it is highly desirable to find a good diagnostic strategy to predict aggressive MIBC with druggable targets.

Several attempts by genome-wide screening have advanced our understanding of molecular mechanisms involved in bladder cancer progression [12,13,14]. Overexpression of certain genes have been frequently linked with MIBC, such as EGFR, ERBB2/3, FGFR3, PPAR, E2F1, EZH2 and YWHAZ [12,13,14,15,16,17]. Although those findings provide opportunities for clinical applications, most of the suggested cancer-associated genes were defined based on candidate-oriented studies. Some other independent biomarkers which show higher clinical relevance may be still missing, awaiting further investigation. To this end, we have previously developed a novel platform, known as HDMAC (High-Dimensional Analysis of Molecular Alterations in Cancer), to facilitate multivariate regression study with a huge number of altered genes during cancer development [18]. Through HDMAC platform, several penalized regression models like ridge, lasso and adaptive lasso are offered for high-dimensional data analysis with Cox regression for survival and logistic regression for binary outcomes. In addition to the increased statistical strength, HDMAC also provides a method to control the false discovery rate (FDR) as the first-step screening. Newly defined biomarkers could be found for predicting aggressive bladder cancer through high-dimensional analyses on alterations of entire annotated genes.

In this study, we aimed to define a new panel of novel biomarkers of which overexpression can predict poor clinical outcomes in patients with MIBC by using HDMAC platform. The prognostic values of those biomarkers can be validated in two independent study cohorts. Furthermore, additive effects can be found that patients with more biomarkers positive show much shorter survival times than those with only one or none. Although no direct molecular interactions were defined, overexpression of those genes contribute to several similar down-stream pathways, suggesting common features involved in aggressiveness of MIBC. Therapeutics based on those target genes provide another opportunity to develop better treatments for patients with aggressive MIBC. Through matching the molecular features of drug responses in cancer cells with gene dependent cell vulnerabilities [19,20,21], several novel drugs were also defined for possible applications against aggressive MIBC.

2. Results

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

2.1. Potential Biomarkers for Aggressive MIBC Based on HDMAC Analysis

To find new biomarkers associated with shorter overall survival (OS) times in MIBC patients, we performed high-dimensional regression of alterations in gene expression (z-scores > 2.0 or <-2.0) in MIBC (Cell 2017, TCGA; n = 367) [14] by using HDMAC platform [18]. The data processing included univariate analysis with Cox-regression (1279 genes showed p-value < 0.05), followed by FDR correction (p-value < 1.0e-7) to delete false discovery univariates, resulting in 145 genes selected (Figure 1A). Then, those selected candidates were further filtered by using penalized regression models, lasso (L) (from 145 to 54 genes) and adaptive lasso (AL) (from 54 to 26 genes), to delete unrelated covariates (Figure 1B). Among the 26 selected genes, the expression of TCEANC, DOK7, C6ORF62, CTU1, BCL2L14, PDCL3 in MIBC showed negative AL coefficients, indicating their possible roles in tumor suppression [28,29,30]. Positive AL coefficients were detected for other 20 genes, suggesting their associations with worse clinical outcomes (Figure 1C). Such conclusion was further confirmed by comparing the tendency between AL coefficients with univariate coefficients (Figure 1D). We also noticed that overexpression of other known bladder cancer genes were highly associated with the expression of these selected genes (Supplementary Figure S1), whereas they were eliminated after FDR correction and penalized regression. Our data therefore suggest diagnostic potentials of overexpression in 26 genes in MIBC (Supplementary Table S1).

2.2. A Six-Gene Panel Serves as Novel Biomarkers for Predicting Poor Clinical Outcomes in Patients with MIBC

To further evaluate the clinical relevance of the selected 26 genes, we analyzed their associations with different clinical features. As shown in Figure 2A, overexpression of the selected genes were found frequently associated with tumor stage, subtype, histological grade and nodal status, suggesting their potent roles in promoting aggressive MIBC. To validate prognostic potentials of those genes, gene expression data from a Korean cohort study (GSE13507; n = 62) [12] were utilized to correlate gene expression levels with the hazard ratios by using a generalized additive model (GAM). As compared to the TCGA cohort, 18 out of the selected 24 genes (75.0 %) (PCLO and TCEANC were not annotated in the GSE13507 study) showed consistent tendency with patient overall survivals in the Korean cohort (Figure 2B, Supplementary Figure S2). Considering possible false conclusion for data with zero inflated, a panel of six genes (ARID3A, ARMH4, P4HB, PPT2, PRPF19 and SLC1A6) was suggested for further investigation. Although PPT2 overexpression did not show any association with the defined clinical features, it was still selected due to its significant p-values of GAM analyses in both cohorts (p = 0.0001 in TCGA, 0.0296 in Korean).

Based on the expression levels, we can successfully predict the significant differences in OS outcomes in both cohorts (Supplementary Figure S3). Strikingly, this six-gene panel also exhibits strong additive effects that patients negative to anyone of those biomarkers show much longer OS times than the ones positive to the defined biomarkers (Figure 3A). Furthermore, patients positive to more biomarkers (> 2) have shorter OS times (with median OS times of 8.94 and 6.87 months in TCGA and Korean cohorts, respectively) than the ones positive to only one or two biomarkers (19.78 months in TCGA and 15.40 months in Korean) (Figure 3B). Those results indicate high prognostic values of this six-gene panel for predicting outcomes of MIBC patients.

2.3. The Six-Gene Panel Shares Common Functional Features/Pathways

To confirm the key pathways governed by the six-gene panel in MIBC, gene-set enrichment analysis (GSEA) was performed using the transcriptome data (Supplementary Tables S2–S7). Interestingly, several functional features/pathways were commonly shared by different biomarkers, even though HDMAC selected candidates usually do not show direct molecular connection after high-dimensional regression (Supplementary Figure S4A). In addition, patients with positive scores to each biomarker seldom overlap with each other (Supplementary Figure S4B). Based on the top-20 enriched functions/pathways shared at least by two biomarkers, they can be classified into five major categories, including energy/nutrition homeostasis, ER-Gogi transportation & ER stress, glutamate & neuronal signaling, cancer-related signaling, and glycome alterations for tumor micro-environment (TME) remodeling (Figure 4A). Most of biomarkers participate in energy/nutrition homeostasis, cancer-related signaling, and glycome alterations for TME remodeling, which were also the top-3 categories activated in MIBC with poor clinical outcomes (Figure 4B). Of note, alterations in N-glycosylation profiles have been emerging as an important topic for developing diagnostic and therapeutic agents [31,32]. Glutamate & neuronal signaling is another interesting finding for MIBC, suggesting the neuronal features of aggressive MIBC [13,14,33]. Using the six-gene panel as the indicator, patients who show positive to this panel (30.1 %) highly express gene markers for neuronal and luminal subtypes, but they are negative to gene markers for basal squamous subtype (Figure 4C), suggesting different etiologies between neuronal/luminal and basal squamous subtypes.

2.4. PPT2 Levels Determine Vulnerabilities of Cancer Cells

To know whether expression of the six-gene panel determine vulnerabilities of cancer cells, we analyzed gene dependency scores after Crispr/Cas9 gene knock-out in different cell lines (Depmap; https://depmap.org/portal/) and correlated the effects with expression levels of the defined genes. As shown in Figure 5A, PPT2 was found as the only gene whose expression levels determined cell survivals (p = 0.0263), suggesting the potential addiction to PPT2 expression in bladder cancer cells. Using PPT2 as a target, four potential drugs by AUC (area under the curve) and two by IC₅₀ analyses were screened out from Cancer Cell Line Encyclopedia database (CCLE; https://sites.broadinstitute.org/ccle/) based on their differential therapeutic effects on cells with different PPT2 levels (Figure 5B,C). Anchusin (also known as alkannin or shikonin; a chemokine receptor inhibitor) and patupilone (EPO906 or epothilone B; a microtubule-stabilizing agent) were considered as useful drug candidates for clinical applications due to their wider ranges of differential AUC values or effective drug concentrations. Notably, the expression levels of CCL7 and CCR3 correlate with PPT2 expression levels in MIBC (Figure 5D). Thus, one of possible targeting pathways associated with anchusin treatment would be the CCL7-CCR3 axis, an axis which has been previously reported to drive cancer metastasis [30,31,32,33,34,35,36]. On the other hand, PPT2 expression also correlated with enhancement of several beta-tubulins, including TUBB2A, TUBB3, TUBB4B and TUBB6 (Figure 5E), which could be the targets of patupilone.

2.5. PRPF19 and Its Network Is a Potent Target for Treating Aggressive MIBC

Among those selected genes, PRPF19 (also known as PSO4) was found highly expressed in all bladder cancer cell lines tested, and all the cell lines show high gene dependency on PRPF19 (Figure 5A). The mean gene dependency score was significantly lower (less than -2.0; p < 0.0001) as compared to other genes (Figure 6A). Since the median dependency score of known oncogenes is around −1.0 [27,37], our data highlight the potency of PRPF19 overexpression in MIBC development. PRPF19 is a key gene that can form a variety of PRPF19 complexes in cells. Among the defined binding partners, PRPF19 is best known to form the core component of spliceosome with PRPF8, CDC5L, BCAS2, PLRG1, and CDC40 [38,39]. Interestingly, Crispr/Cas9 knock-out on those complex partners also showed high essentiality for cell survival, especially for PRPF8 and CDC5L (Figure 6B). Clinical study indicated that PRPF19 expression levels were significantly higher in bladder cancer tissues than in adjacent normal counterparts (Figure 6C; p < 0.0001). In addition, PRPF19 expression showed high biological relevance that the DNA content correlated nicely with its gene expression levels in MIBC samples (Figure 6D; p < 0.0001). Those features strongly suggest PRPF19 as a potent cancer driver, which can serve as an ideal target for new drug development.

To address molecular bases of the selective dependency on PRPF19 overexpression, we screened cell-killing patterns of genes that were co-upregulated with PRPF19 in MIBC (Spearman's correlation > 5.0; Supplementary Table S8). As shown in Figure 6E, PRPF19 expression levels negatively correlated with gene effect scores of POLA2 knockout, suggesting these two proteins are functionally relevant [40]. In addition to a strong association with PRPF19 expression (Figure 6F), POLA2 expression level also correlate with its DNA copy number (Figure 6G). Furthermore, patients with overexpressed POLA2 levels showed shorter overall survivals as compared to patients with low POLA2 levels (Figure 6H), suggesting a similar pro-cancer role as PRPF19 does.

Interestingly, we noticed that the essential genes with gene effect scores less than 2.0 in human bladder cancer cells contribute to five major complexes/functions, including ribosome, spliceosome, cell cycle, proteasome, and DNA replication (Figure 7A). Of note, PRPF19 plays an important role to physically link spliceosome with DNA replication complex in the essential gene network via PRA1, where POLA2 is one of key genes in the DNA replication complex. To explore possible drug candidates with the molecular actions of anti-PRPF19 or anti-POLA2, we correlated the gene effect patterns with drug-killing profiles across 578 human cancer cell lines by using PRISM dataset with 4,694 compounds (https://depmap.org/repurposing) [19,20,21]. Through supervised linear regression analyses, our study revealed selumetinib and STA-5326 (also known as apilimod) as possible therapeutic leads against aggressive MIBC overexpressing PRPF19 as well as BAN-ORL-24 and STA-5326 for the ones overexpressing POLA2 (Figure 7B). Among them, STA-5326 treatment was found effective for both of PRPF19 or POLA2-overexpressing cells, suggesting the involvement of STA-5326 effects on spliceosome and DNA replication complexes.

3. Discussion

In this study, we have identified 26 genes whose overexpression can define aggressive MIBC in the TCGA cohort (Cell 2017) by using a statistic platform, HDMAC, for high dimensional regression analyses (Figure 1) [18]. Through clinical association study and cross-validation with a Korean cohort (GSE13507), a six-gene panel (ARID3A, ARMH4, P4HB, PPT2, PRPF19 and SLC1A6) was further defined as a novel signature for MIBC patients with shorter survivals (Figure 2 and Figure 3). Five major pathways were commonly shared by those six genes in MIBC, which show molecular features more similar to neuronal and luminal, but not basal squamous subtypes as defined by transcriptome classification (Figure 4) [13,14,33]. Gene dependency scores based on Crispr/Cas9 gene knock-out studies indicated PPT2 and PRPF19 as ideal targets for therapeutic interventions whose expression can determine vulnerabilities of cancer cells (Figure 5 and Figure 6). Anchusin and patupilone were identified as candidate leads for treating patients with PPT2 overexpression (Figure 5). POLA2 was found co-upregulated with PRPF19 overexpression with functional consequences by gene knockout (Figure 6). In addition, both genes were amplified and overexpressed in MIBC with a clear biological relevance from DNA copy numbers to gene expression levels, which can determine clinical outcomes of the patients (Figure 6). Notably, gene knockout of either genes contributed to strong vulnerabilities of cancer cells, suggesting the involvement of the interplay between splicesome and DNA replication complex in MIBC development of (Figure 7). Through drug mechanism-of-action investigations, Selumetinib and STA-5326 were suggested as good lead compounds for treating patients with PRPF19 overexpression (Figure 7B). STA-5326 was also found feasible for treating patients with POLA2 overexpression. Our findings therefore provide novel biomarkers for patient stratification and possible strategies for clinical interventions.

PPT2 is a lysosomal enzyme with thioesterase activity that removes fatty acyl groups from both long and short chain fatty acyl CoA substrates [41]. With homologous sequence to LTBP (latent transforming growth factor-beta-binding protein) as well as the binding capability toward EGFL8 and FGFR2 [42,43], PPT2 has been proposed to participate in a variety of oncogenic signaling, including TGF-β, β-catenin/WNT and EGF signaling [44]. In addition, PPT2 can also interact with the ILF2/ILF3 complex [45], which is required for mitotic progression and cytokinesis by regulating the expression of a cluster of mitotic genes by competing with staufen-mediated mRNA decay [46]. Notably, this complex was recently reported to play roles in cancer metabolic reprograming through inflammatory signaling cascades, especially enriched in phosphatidylcholine (PC) biosynthesis, fatty acid metabolism, and amino acid metabolism pathways [47]. In our study, we also found that patients with higher PPT2 levels show upregulated FGFR2 signaling, cell cycle activity and PC biosynthesis (Figure 4A). These findings may explain the potential benefits of patupilone (as a microtubule stabilizer) and anchusin (a suppressor of CCL7-CCR3 axis) treatments for cancer patients with PPT2 overexpression, which promotes cell cycle progression and inflammation associated metabolic adaptions.

PRPF19 forms a core complex with CDC5L, known as Prp19C, inside the spliceosome complex crucial for spliceosome assembly and catalytic activation of RNA splicing process [48,49]. Emerging evidence indicate multifaceted roles of this complex in diverse cellular processes, including DNA repair/genome maintenance, DNA replication/cell division, transcription elongation, protein ubiquitination, and lipid droplet biosynthesis [38,39,50,51,52,53]. Upregulation of PRPF19 overexpression is therefore linked with the development of several types of cancers [54,55,56,57,58], and the PRPF19/p38MAPK axis was found as the key pathway to enhance tumor aggressiveness via promoting EMT [57]. In fact, p38MAPK also plays an important role to balance the growth and death signals in response to cellular stresses, including blockage of autophagy-lysosomal system [59,60,61]. This may explain the effectiveness of Selumetinib, a potent MAPK inhibitor, for treating cancer cells with upregulated PRPF19 (Figure 7B).

On the other hand, PRPF19 can function as an E3 ligase to promote ATR-mediated DNA repair and homologous recombination, leading to chemo-/radio-resistance [54,62,63,64,65]. At clinics, co-upregulation of both PRPF19 and POLA2 was found as an important signature for mesothelioma patients with platinum resistance [66]. Our study revealed that PRPF19 overexpression in MIBC was found associated with the enrichment of proteasome activity, protein degradation during mitosis and chromosomal segregation (Figure 4A), supporting the potent roles of PRPF19 as being an E3 ligase during cancer development. STA 5326 (Apilimod) was originally discovered as an IL-12/IL-23 inhibitor [67], whiles this compound was later on identified as a PIKfyve inhibitor that can suppresses the formation and function of lysosomes, leading to impaired energy homeostasis and mitotic autophagy during cytokinesis [68,69,70]. Since PRPF19 controls the turnover of several proteins involved in DNA repair, replication and chromosome segregation through a phosphorylation-and-ubiquitylation circuitry, STA 5326 is therefore a good drug candidate for treating MIBC with PRPF19 and POLA2 (Figure 7B). Interestingly, combination of p38MAPK and PIKfyve inhibitors can synergistically reduce tumor growth of colon cancer in mice by disrupting autophagy, whiles no obvious effects were found on the viability of normal cells [61]. The potential benefits of combined therapy of Selumetinib and STA 5326 for MIBC cells are under investigation.

Significantly, both druggable targets suggested in this study, PPT2 and PRPF19, are involved in protein degradation pathways, highlighting the crucial roles of autophagy-lysosomal and proteasomal systems in malignant transformation and development of MIBC. Increased numbers of lysosomes or enhanced activity of lysosomal/proteasomal enzymes have been previously reported as an important feature to increase the invasiveness and drug resistance in bladder cancer cells [71,72,73,74]. Several proteins involved in protein degradation, such as Nedd4 like E3 ligases, have been proposed as ideal targets for new drug development against bladder cancer [60,73,75]. PRPF19 contains WD40 repeats, enabling it to function as a RNA-binding ubiquitin ligase. Through connecting post-transcriptional regulation with the ubiquitin system, PRPF19 participates in diverse sub-complexes to regulate gene expression and protein translation [38,39]. With more potent roles than other E3 ligases in oncogenesis, PRPF19 may serve as a novel and efficient lethal target for cancer treatment, especially for MIBC. Compounds that block the interface of WD40 repeats in PRPF19 (PDB 4LG8) with its substrates might be potential leads for PRPF19-based therapy [39] which will be investigated in our future study.

In this study, the proposed six-gene panel is highly associated with neuronal subtype with neuroendocrine-like features, even in the absence of neuroendocrine histology (Figure 4C). Previously, similar high-content drug screening was performed by using PRISM Repurposing dataset to define subtype-stratified novel lead compounds against MIBC [76]. Although several novel options were suggested with subtype-specific actions for basal, luminal and mesenchymal-like subtypes, none was suggested for treating MIBC neuronal subtype. Interestingly, Kim et al. indicated that MIBC patients defined as neuronal subtype showed extraordinarily high respond rate toward ICI-based therapy [77]. However, such conclusion was not supported by immune scores of TME with some concerns in the reliability of prediction with a small sample size (n = 11) [78,79,80]; in reverse, MIBC patients with neuronal phenotype likely obtained less benefits from immunotherapy due to immune-excluded features and low PD-L1 expression [13,33,78]. Investigations of ICI-based clinical trials also support the point of view that tumors with neuroendocrine features have been associated with an immune-desert TME, limiting the therapeutic effects of ICI-based therapy [81,82,83]. The findings from our study thus provide a novel strategy for treating patients with aggressive MIBC who might be ICI-non-responders by targeting PPT2 or PRPF19.

While our study has uncovered novel prognostic biomarkers for the most aggressive subtype of MIBC, along with several proposed therapeutic strategies, it is essential to acknowledge certain limitations. Firstly, our research identified racial disparities in the gene expression of selected biomarkers and their associated cancer outcomes. To validate these suggested biomarkers, it is imperative to include more study cohorts with diverse ethnic backgrounds. Furthermore, the determination of the cut-off point for each biomarker may vary according to race and ethnicity. Additionally, the functional impacts of the selected biomarkers, as well as the therapeutic potentials of the suggested lead compounds in our study, require further confirmation through gene knock-in/out cell-based or animal-based studies. Particularly, the underlying mechanisms behind the enhanced activity of protein degradation machineries in determining the aggressiveness of MIBC remain unclear and necessitate further exploration through additional studies.

4. Materials and Methods

4.1. Patients and Gene Expression Data

Clinical data and gene expression data of patients diagnosed as MIBC were collected from TCGA cohort (Cell, 2017) [14] and an Asian cohort (GSE13507) [12]. Patients with incomplete gene expression data, clinical information, and survival statuses were excluded. For TCGA cohort (n = 413), 367 patients were diagnosed as MIBC with complete genetic and clinical information. The RNAseq data normalized with diploid samples as z-scores were downloaded from cBioportal (https://www.cbioportal.org/). To compare the differential expression levels of a defined gene between cancer and normal tissues, normalized RNA-seq data were obtained from TCGA and GTEx projects via GEPIA2 server (http://gepia2.cancer-pku.cn/#index). For the GSE13507 cohort (n =165), 62 patients were diagnosed as MIBC with complete genetic and clinical information. Gene expression data based on Illumina human-6 v2.0 gene beadchip platform were downloaded from GEO databank (https://www.ncbi.nlm.nih.gov/geo/). LIMMA (Linear Models for Microarray Data) R package was utilized to estimate the expression level of each gene by averaging data from multiple probes. The resultant data were converted to z-scores through quantile normalization using data from normal tissues as the references and log₂ transformation.

4.2. Potential Biomarkers for Aggressive MIBC Based on HDMAC Analysis

The HDMAC platform (https://ripsung26.shinyapps.io/rshiny/) was utilized to define new biomarkers for MIBC patients [18]. Alterations in gene expression (z-scores > 2.0 as overexpression or <-2.0 as down-regulation) with frequencies more than 5% among the patients were considered as the key events. In this HDMAC package, the platform can compute the significance between the expression level and patients’ survival, and select the most significant candidate genes through univariate analysis by proportional hazards Cox model with the FDR-adjusted p-values less than 5% [22], followed by multivariate analyses with lasso and adaptive lasso Cox’s proportional hazards models [23]. A generalized additive model was used to investigate the linear/nonlinear effect of gene expression on patients’ survival [24]. This approach enables us to identify target genes highly relevant to patients’ survival.

4.3. Clinical Association Study and Statistical Analyses

Clinical data of MIBC patients, including gender, age, smoke, stage, subtype, nodal status, etastasis status, histological grade and angiolymphatic invasion were retrieved from TCGA dataset. GraphPad Prism (GraphPad software, La Jolla, CA) was utilized for statistical analyses. Statistical differences between two groups were estimated by independent sample t-test and differences among more than two groups were estimated by one-way analysis of variance (ANOVA). Kaplan–Meier method was used to analyze the survival curves, and a log-rank test was performed for statistical comparison. To obtain an easily explainable multivariate survival model for clinical use, we identified the optimal cutoff point for gene expression levels of each target gene by Cox proportional hazard model [25]. The optimal cutoff points were defined as the ones with the most significant (likelihood ratio test or log-rank test) split. Pearson coefficient (r) and Spearman's rank correlation were used to measure a linear correlation of gene expression levels between two genes.

4.4. Gene Set Enrichment Analysis (GSEA)

To map the key pathways associated with overexpression of the defined biomarkers, TCGA cohort was divided into altered (with gene overexpression) and unaltered (without gene overexpression in any gene of the 6-gene panel) groups. Normalized gene expression data were applied for enriched pathway analyses by using GSEA Java (version 4.1.0) from Broad Institute, Cambridge, MA [26]. The significance of enrichment scores (ES) were normalized by permutation test, followed by FDR and FWER corrections. Over-representation analysis was performed to analyze the molecular linkage of target genes with genes enriched in the defined pathways.

4.5. Gene-Dependent Cell Vulnerabilities

The gene expression levels of candidate genes in 38 BLCA CCLE cell lines were downloaded from the Cancer Dependency Portal (DepMap) (https://depmap.org/portal/ccle/). Cell viability data after CRISPR–Cas9 knock-out were assessed based on the genome-scale data from the Project Achilles (DepMap Public 22Q4). Gene effect scores were calculated and converted to CERES score by normalizing the effects with gene copy number using the CERES algorithm [27]. The dependency scores less than -1.0 were identified as essential for cell survival. For genes with dependency scores larger than -1.0, gene-dependent cell vulnerabilities were analyzed based on the correlations between gene expression level of a particular candidate gene and the associated gene effect score in each cell line. A correlation plot with a Pearson coefficient (r) less than -0.3 was considered as positive for gene-dependent cell vulnerability.

4.6. High-Content Drug Screening

The drug response data including AUC (319 anti-cancer drugs) and IC₅₀ values (206 drugs) were accessed by using drug dilution curves of the CCLE Cancer Dependency Portal (DepMap; GDSC1). The AUC values were presented on a scale between 0 and 1, where lower AUC values indicate increased sensitivity to drug treatment. Through fitting the curves to cell viability for each drug and cell line, the gene-drug response relationships were plotted and calculated according to expression levels of the defined candidate gene in different BLCA cell lines. Effective drugs were ranked by their differential AUC or IC₅₀ values and tested for significance by unpaired two-tailed Student’s t-test with Welch correction. Effective drugs with wider ranges of differential AUC values (> 0.6) or effective drug concentrations (>10 digits) were considered as potential leads in this study.

4.7. Drug Screening by Mechanism-of-Action Matching

The drug responses of 578 human cancer cell lines to 4,694 compounds were assessed fromPRISM drug databank (https://depmap.org/repurposing). Supervised linear regression was performed to fit the best association between the cell viability and the gene effect scores by CRISPR knockout (DepMap Public 22Q4). Positive compounds (p < 0.05) were selected and subjected to likelihood ratio test or log-rank test to define the optimal cutoff points with the most significant split. In this study, the selected cutoff points for gene effect score were -2.40 for PRPF19 and -1.80 for POLA2; the cutoff points for cell viability were -1.78 for Selumetinib, -1.37 for STA-5326 and -0.75 for BAN-ORL-24. The matching association was further analyzed by 2 x 2 chi-square by using the defined cutoff points.

5. Conclusions

Utilizing the HDMAC platform, we have identified a six-gene panel capable of predicting worse clinical outcomes in patients with MIBC, which has been subsequently validated in an independent cohort. Patients exhibiting higher expression levels of these biomarkers experienced significantly shorter survival times compared to those with fewer biomarkers. Our functional screening efforts pinpointed PPT2 and PRPF19 as potential druggable targets involved in protein degradation via autophagosome/lysosome or proteasome pathways. Consequently, our study offers a promising diagnostic strategy for identifying aggressive MIBC cases, along with actionable targets for therapeutic intervention.