Developing angiogenesis-related prognostic biomarkers and therapeutic strategies in bladder cancer using deep learning and machine learning

Establishment of ARGS model and regulation mechanism of BLCA angiogenesis

Differentially expressed analysis identified 147 angiogenesis-related genes (98 downregulated, 49 upregulated). A heatmap subsequently visualized ARGS expression patterns in BLCA tissues and normal bladder tissues (Fig. 2A). The differential expression of down and upregulated angiogenesis-related genes is presented in the volcano map (Fig. 2B). GO enrichment analysis revealed that these genes are involved in the regulation of endothelial cell, epithelial cell, cell migration, and tumor angiogenesis. These pathways showed a downregulation trend, suggesting a high regulation of tumor angiogenesis and cell migration (Fig. 2C). KEGG enrichment analysis indicated a decreased trend in pathways related to PI3K-AKT signaling pathway and ECM-receptor interaction (Fig. 2D). Unicox regression analysis identified 41 risk genes associated with BLCA prognosis (Fig. 2E), and LASSO regression further reduced the dimensionality to 12 genes used to construct the ARGS scoring model (Fig. 2F). Among them, CALR, COL14A1, EMP1, ENO1, HMGA1, PDGFRA, SLC3A2, and TCF4 exhibited higher expression levels in the high-risk group, while GNG5, HSPE1, VHL, and MST1R exhibited higher expression levels in the low-risk group (Fig. 2G).

Training of the ARGS model and construction of the BLCA prognostic model

ARGS risk stratification revealed significantly shorter OS in high-risk patients versus low-risk counterparts (P < 0.05) (Fig. 3A). Risk scores demonstrated a significant negative correlation with OS duration. The area under curve (AUC) for 1, 3, and 5-year survival indicated that in both the training and validation cohorts, the ARGS model effectively assessed the impact of angiogenesis on survival in BLCA patients (AUC > 0.6, Risk > 0.6) (Fig. 3B). The clinical relevance analysis of angiogenesis and the analysis of independent risk factors yielded results consistent with the AUC model, and the differences were statistically significant in the training, control, and validation models (P < 0.001) (Figs. 3E and 4D). PCA and t-SNE results demonstrated that the ARGS model effectively stratified BLCA patients, with a well-dispersed pattern (Fig. 3C). KM survival analysis revealed that patients in the high-risk group had significantly shorter OS compared to patients in the low-risk group (P < 0.001) (Fig. 3D). Consistent results were observed in both the training and validation cohorts.

Prognostic significance of the ARGS

Taking the median ARGS risk score as the threshold, we divided the BLCA patients into two groups, the high-risk and low-risk groups (Fig. 3C). The heatmap shows the differential expression of 12 essential genes between two ARGS subgroups (Fig. 3A). The risk scores and surveillance profiles of each BLCA patient are shown in Fig. 3B. KM analysis demonstrated that in both the treatment and validation cohorts, the OS and meditation duration of survival of patients in the high-risk cohort were lower than those in the low-risk cohort (P < 0.05) (Fig. 3D). The AUC values at 1, 3 and 5 years were 0.695, 0.696 and 0.738 in the TCGA cohort, as indicated in Fig. 3C. In the GEO cohort, the AUC values at 1, 3 and 5 years were 0.712, 0.729 and 0.731. The KM analysis showed that the patients had worse OS and median survival times than those in the low-risk cohort (P < 0.05) (Fig. 3D). The correlation between ARGS scores and clinicopathological factors was elucidated in the TCGA cohort (Fig. 4D). These results suggest that the construction of the ARGS score model plays a crucial role in the progression for BLCA and is a good predictor of patient prognosis in BLCA.

Nomogram for survival probability in BLCA patients

To accurately predict the accuracy of OS models for different subgroups of patients, risk scores and other clinicopathologic characteristics of BLCA patients, including age, sex, grade, and stage, were integrated by constructing a column-line graph of prognostic risk-related models (Fig. 4B). We assessed the effectiveness of the KM survival curves and the validity of the ROC phantom based on the prophylactic risk score and the nomogram risk score at 1, 3, and 5 years, respectively (Fig. 4C). The calibration curves showed that the predicted and actual values fit well, and the model exhibited good predictive validity for the survival of actual BLCA patients. The nomogram AUC was 0.773 and 0.851 and the risk AUC was 0.747 and 0.786 for the trial and validation groups (Fig. 4C), respectively, and the results indicated that the column line plot was effective in predicting BLCA survival rates (Fig. 3B). The ARGS scoring model still demonstrates robust predictive ability for the prognosis of BLCA patients across various machine learning algorithms (Fig. 4A).

Enrichment analysis were conducted to explore the functional implications of the ARGS model

In GSVA, the risk score of each sample showed that the low-risk group was mainly enriched in the peroxisome, linoleic acid metabolism, and alpha-linolenic acid metabolism, which are pathways primarily related to metabolism. In contrast, the high-risk group was enriched in ECM-receptor interaction, B cell receptor signaling pathway, T-cell receptor signaling pathway, Toll-like receptor signaling pathway, and other immune-related pathways (Fig. 5A). In GO enrichment, the primary BP of prognostic risk angiogenesis-related genes is involved in peptide antigen assembly with MHC class I protein complex (GO:0002502) and phenylalanine transport (GO:0015823). For CCs, angiogenesis-related genes were mainly present in the collagen type XIV trimer (GO:0005596), the apical pole of neurons (GO:0044225), and the cell pole (GO:0060187). The MF of angiogenesis-related genes were mainly enriched in platelet-derived growth factor alpha-receptor activity (GO:0005018), structural constituent of chromatin (GO:0030527), and macrophage colony-stimulating factor receptor activity (GO:0005011) (Fig. 5B, C). KEGG analysis showed that prognosis-related risk angiogenesis-related genes were associated with the HIF-1 signaling pathway (ko04066), RAS signaling pathway (ko04014), and PI3K-AKT signaling pathway (ko04151) (Fig. 5B, C). The above results suggest that prognostic risk angiogenesis-related genes are mainly associated with immunity, cell proliferation, and vascular growth. Therefore, it can be inferred that in BLCA, tumor angiogenesis may be induced by risk genes regulating the HIF-1/PI3K-AKT/RAS signaling pathway, and the tumor immune microenvironment is closely associated with this process (Fig. 5B, C). T cells and B cells are likely to be the major participants in this process. The 18 transcription factors (KLF9, KLF2, HAND2, HLF, GATA6, FOXF1, EGR3, EGR2, EBF1, ZEB1, ZBTB16, SOX5, SOX17, SOX10, RUNX1T1, NR4A3, MYOCD) associated with ARGS are primarily related to tumor cell proliferation (Fig. 5C).

Expression, prognostic value, and GSEA analysis of ARGS genes in BLCA

The KM survival analysis showed that overexpression of PDGFRA, TCF4, EMP1 and low expression of MST1R, HSPE1, VHL suggest poor prognosis in BLCA patients (P < 0.05) (Figure S1), consistent with Fig. 2E. PDGFRA, TCF4, and EMP1 are risk factors for BLCA, while MST1R, HSPE1, and VHL are protective factors for BLCA. However, no significant differences were observed in the impact of mRNA expression levels of SLC3A2, HMGA1, GNG5, ENO1, and CALR on the prognosis of BLCA patients. The IHC map showed that except for PDGFRA, VHL, EMP1, and GNG5, which did not obtain slices, the expression levels of other genes in tissues were similar to those obtained from KM survival analysis and risk gene analysis. CALR, COL14A1, SLC3A2, ENO1, and TCF4 were highly expressed in BLCA tissues, while HMGA1, MST1R, and HSPE1 were highly expressed in normal bladder and urinary tract epithelial tissues (Figure S2).

GSEA identified significant associations between ARGS genes and critical oncogenic pathways (Figure S3). CALR, PDGFRA, and EMP1 were enriched in cytokine-cytokine receptor interaction and ECM-receptor interaction, suggesting their roles in tumor microenvironment remodeling and immune modulation. Notably, COL14A1, MST1R, SLC3A2, TCF4, and ENO1 clustered at the leading edge of ECM-receptor interaction, implicating their potential involvement in extracellular matrix dysregulation and tumor invasiveness.

Conversely, GNG5 showed inverse enrichment in natural killer cell-mediated cytotoxicity and xenobiotic metabolism, hinting at immune evasion and chemoresistance mechanisms. HMGA1 and VHL were linked to cell cycle, cytokine signaling, and JAK-STAT pathways, aligning with their known roles in proliferation and cell-cycle dysregulation. Additionally, HSPE1 was highly enriched in oxidative phosphorylation, possibly reflecting metabolic reprogramming in BLCA progression. These findings highlight the multifaceted regulatory roles of ARGS genes in BLCA pathogenesis, spanning immune suppression, ECM remodeling, and metabolic adaptation.

Immune cell correlation analysis

Taking the median risk score of BLCA patients as a threshold, risk correlation analysis of immune-related cells and TME scores was performed with various software, and 12 types of immune cells were highly correlated with risk scores(Fig. 6A). The results showed that increased levels of the infiltration of immune cells, such as M0 macrophages, neutrophils, and M2 macrophages, were strongly associated with poor prognosis in patients in the high-risk group. In ARGS, there are mainly cancer associated fibroblasts, M1 macrophages, M2 macrophages, macrophages, monocytes, monocytes, activated myeloid dendritic cells, myeloid dendritic cells, neutrophils, CD4 + Th2 T cells, CD8 + T cells and regulatory T cells (Tregs) (Fig. 6D), and the increased level of infiltration of immune cells described above in the high-risk group suggests that the poorer prognosis of BLCA patients in the high-risk group correlates with the increased level of infiltration of immune cells in the TME (Fig. 6F).

The t-SNE results showed that the grouping model could evaluate the differences between the high-immunity group and the low-immunity group very well (Fig. 6B). Furthermore, both the immunity score and stromal score were lower in the low immunity score cohort than in the high immunity score cohort (Fig. 6C). Significant immunological features were found in the high and low immunity scoring cohorts, with statistically meaningful variations, as demonstrated in Fig. 6A. Immunophenotyping analysis was performed on all samples, and ssGSEA scores based on 29 immune genomes were categorized into four subtypes (referred to as C1-C4) (Fig. 6G), including wound healing (C1), IFN-γ dominance (C2), inflammation (C3), and lymphocyte depletion (C4). Immunological infiltration decreases in the decreasing sequence: C2 > C1 > C3 > C4. For example, C2 is dominant. C1 was mainly clustered in the high-risk group of patients. C3 and C4, however, were mainly clustered in the low-risk group of patients (Fig. 6G). Accordingly, patients in the low-risk category are better equipped to be immunotherapy beneficiaries. While patients in the high-risk bracket are refractory to immunotherapy. As far as the high-immunity subgroup and low-immunity subgroup are involved, the immunotherapy outcomes and prognosis of the high-immunity subgroup are better as compared to the patients in the low-immunity subgroup.

After analyzing the immune checkpoints, we found that TNFSF15, TNFRSF25, TNFRSF14, TMIGD2, and LGALS9 were negatively correlated with the risk scores. This suggests that the costimulatory activation of the TNF receptor and ligand superfamily was suppressed (Fig. 6E), and the downregulation of VHL, HSPE1, and MST1R in the tumor tissues was strongly correlated with suppression of TNF receptor and ligand superfamily. According to TCIA immunotherapy analysis data (https://tcia.at/)²³, when CTLA-4 and PD-1 were blocked, a comparison of immunocyte positive scores (IPS) revealed that patients in the low-risk group had a better immunotherapy outcome and prognosis than patients in the high-risk group. The percentage of tumor-associated immune cells was higher than any intensity of PD-L1 membrane- and cytoplasm-stained tumor-associated immune cells, as determined by the percentage of all tumor-associated immune cells. Patients with BLCA in the high-risk group received both CTLA-4 and PD-1 blockade as the optimal immunotherapy strategy. Moreover, patients in the low-risk group have a better immunotherapy outcome than patients in the high-risk group when PD-1 alone is blocked, and a poor immunotherapy outcome is a crucial factor in the poor prognosis of patients in the high-risk group (Fig. 6H).

Drug sensitivity analysis and evaluation of biomarker selection

To identify therapeutic candidates for BLCA, we conducted drug sensitivity analysis using the prognostic ARGS risk signature. We finally obtained a total of 18 known drugs with the potential to reverse poor prognosis in BLCA patients (R > 0.4), such as TGX221, Saracatinib, and pazopanib. WH-4-023 is the most likely to reverse poor prognosis in BLCA patients (Fig. 7A). STOCK1N-35,696, was identified as a novel potential therapy for BLCA (Fig. 7B). Drug sensitivity analysis showed that chemotherapy is less effective in patients in the high-risk group than those in the low-risk group.

The PCA model displayed the results of batch normalization for five groups of BLCA and normal bladder tissue datasets, indicating good separation of the samples (Fig. 8A). The BATCH algorithm was used to select genes with significant differences between BLCA tissue and normal tissue, using a logFC = 1 and P = 0.05 threshold for multiple learning screenings (Fig. 7C). LASSO regression identified 18 potential biomarkers associated with BLCA (Fig. 7D), while SVM selected 25 potential biomarkers with the smallest selection error (ERROR) for change point screening (Fig. 7E). RF selected 45 potential biomarkers with a statistical count of 1000 and gene importance score > 2 (Fig. 7F). WGCNA exhibited a stable topological fit index reaching 0.9 when the soft threshold was set to 7 (Fig. 7G). Key modules were selected, and it was found that the blue module had the strongest correlation with BLCA, encompassing 33 key module genes (Fig. 7D-G). The intersection of the potential biomarkers obtained from the four machine learning screenings yielded 5 common genes (MYH11, FGF9, BIN1, SRPX, JAM3) (Fig. 8B). The AUC curves for the control and validation groups were 0.834 and 1.000, respectively (Fig. 8C), with high predictive validity. The expression changes of these genes in the samples were represented by line graphs, showing significant differences in the expression of potential biomarkers between the control group and the experimental group (Fig. 8E). Figure 8D shows that, in the training group queue, the AUC for MYH11 was 0.782 (95% CI = 0.729–0.827), for FGF9 it was 0.782 (95% CI = 0.732–0.830), for BIN1 it was 0.782 (95% CI = 0.733–0.830), for SRPX it was 0.763 (95% CI = 0.708–0.813), and for JAM3 it was 0.767 (95% CI = 0.717–0.812). Among them, MYH11 had the highest cutoff value, and FGF9 had the highest model sensitivity.

Functional evaluation of biomarkers for anti-angiogenic therapy in BLCA

Consequently, MYH11 was found to be primarily involved in pathways such as vascular smooth muscle contraction, ECM-receptor interaction, and focal adhesion in BLCA. These pathways showed a negative correlation with MYH11 expression, suggesting that the downregulation of MYH11 is associated with tumor proliferation, migration, and angiogenesis in BLCA (Fig. 8F). Similarly, SRPX was found to be involved in signaling pathways such as complement and coagulation cascades, regulation of actin cytoskeleton, and focal adhesion, sharing similar functions with MYH11 and primarily associated with tumor proliferation and migration in BLCA. The downregulation of BIN1 was related to signaling pathways such as complement and coagulation cascades and cell adhesion molecules, which are mostly involved in cell proliferation and migration. FGF9 and JAM3 were also implicated in tumor proliferation and migration in BLCA. These biomarkers can serve as evaluation indicators for the therapeutic efficacy of anti-angiogenic treatment in BLCA, with MYH11 showing the highest sensitivity.

Exploring novel ADC-based treatments for BLCA using homology modeling and pharmacophore modeling

Using the Degree value as a screening criterion, the PPI network revealed the top 20 key genes that induce BLCA angiogenesis. Among them, FN1 is the most likely protein for BLCA angiogenesis and has high potential as a target for BLCA anti angiogenic therapy (Table 1). The protein crystal model of FN1 was derived from the Bos taurus (Bovine).

To identify the critical chemical features of highly cytotoxic ADC payloads targeting angiogenesis in BLCA, we constructed pharmacophore models incorporating three key interaction features: hydrogen bond acceptors, hydrophobic regions, and aromatic ring centers. Among the 10 generated models, Model 1 demonstrated superior alignment with known ADC payload mechanisms, achieving the highest rank (79.289) and FitValue (3.99953). The thermodynamic stability of Model 1 was further supported by favorable energy values (absolute energy = 596.419 kcal/mol; clean energy = 707.921 kcal/mol). These pharmacophore features correspond to empirically validated characteristics of ADC payloads, where hydrogen bond networks mediate target engagement while hydrophobic domains influence payload retention and bystander effects (Table 2). Model 1 included features of hydrogen bond acceptor interaction force and hydrophobic feature interaction force (Fig. 9A).

Using SWISS-MODEL, we generated 10 homology models of human FN1 based on the bovine FN1 template (PDB:1FNF). Model 9 showed superior structural quality, with 97.23% of residues in Ramachandran-favored regions (vs. 70.11–95.97% for other models) and excellent validation scores (MolProbity = 0.77, QMEANDisCo = 0.81 ± 0.05). Despite moderate GMQE (0.13) due to FN1’s complex domain architecture, Model 9 preserved key functional motifs including the integrin-binding RGD domain (Table 3). The overall quality factor analysis of FN1 yielded a coefficient of 80, indicating good predictive ability and a stable protein crystal structure (Fig. 9D). The Ramachandran plot occupancy was 97.23% (> 90%), indicating reasonable dihedral angles for the target protein’s Cα (Fig. 9C). The energy landscape maps for the full residue and Chi1-Chi2 displayed active red residues that matched the amino acid residues involved in the molecular docking, further confirming the potential of FN1 as a cancer target (Fig. 9B).

ADME property analysis revealed that Ouabain, Rutin, SSD, Panaxoside A and Baicalin have poor synthetic accessibility and do not comply with Lipinski’s rule (Fig. 9E; Table 5). Except for Wogonin, the remaining compounds exhibit good solubility and lipophilicity. All compounds have low polarity, indicating stability (Table 5). Among them, Ceastrol has high drug-likeness and low levels of toxic side effects. Its low synthetic accessibility suggests ease of artificial synthesis or extraction.

All compounds targeting FN1 will be retained for the subsequent pharmacophore overlay of ADC companion agents. Comprehensive toxicity evaluation of potential ADC payloads (Table 4) revealed significant variation in safety profiles. Gypenoside III showed the highest toxicity risk (Acute Tox. = 4, STOT RE = 4), followed by Panaxoside A (scores = 4/4) and Genistein (scores = 4/2). In contrast, SSD exhibited favorable toxicity characteristics (scores = 3/3), supporting its selection as the lead candidate. These findings were consistent with molecular docking results, where higher toxicity compounds demonstrated stronger binding affinity but poorer safety margins. The toxicity hierarchy (Gypenoside III > Panaxoside A > Genistein > SSD) informed subsequent prioritization of compounds for ADC development.

Virtual screening of cytotoxic ADC payload candidates using pharmacophore models and molecular docking

Based on Model 1 pharmacophore model, a screening was conducted to identify potential natural compounds with high cytotoxicity for ADC. SSD, Gypenoside III, Panaxoside A, and Celastrol were identified as potential natural compounds targeting FN1 that could potentially serve as high cytotoxicity ADCs. SSD has a FitValue of 2.46255, Gypenoside III has a FitValue of 2.40565, Panaxoside A has a FitValue of 2.23743, and Celastrol has a FitValue of 1.52411, indicating that SSD has the highest potential as an ADC companion agent. After molecular docking, CDOCKER energy analysis revealed that the small molecule ligand SSD had the lowest binding energy (−15.186 kcal/mol) when bound to the receptor protein FN1, indicating the most stable ligand-protein complex (Table 6). In conclusion, SSD has the highest potential as an ADC companion agent and can form stable complexes with FN1. Therefore, SSD is considered to have the ability to act as an ADC companion agent against BLCA angiogenesis.

TOPKA analysis of SSD

The prediction results show that SSD exhibits low-risk characteristics in most toxicity endpoints (Table S1). SSD was clearly predicted as a “non-mutagenic substance” (Ames-test, Bayesian score = − 31.60, P < 0.05) and a “non-carcinogenic substance” (all rodent models, P < 0.05), with Bayesian scores for carcinogenicity in rats and mice ranging from − 8.05 to −0.89, showing significant statistical significance. Acute and chronic toxicity analysis showed that the LD50 of SSD in rats was 4.83 mg/(g·BW·day) (P < 0.05), and the chronic LOAEL was 2.62 mg/kg·BW (P < 0.05), indicating its high safety under acute and long-term exposure. In addition, the biodegradability of SSD was predicted to be “degradable” (Bayesian score = 5.69, P < 0.05), indicating that it may be easily metabolized in the environment and has a low ecological accumulation risk. The MD of all endpoints was significantly higher than the training set heart value (MD > 10), and the P-value was much lower than 0.05, further supporting the high confidence of the prediction results.