The overall analytical workflow for this study is depicted in Fig. 1. A total of 516 ARGs were retrieved from the Genecards and Harmonizome databases (Supplementary Table 1). Subsequent differential expression analysis identified 136 ARGs that were significantly dysregulated between BC and normal tissues (Supplementary Table 2). A heatmap illustrates the expression patterns of the top 50 up- and down-regulated genes (Fig. 2A), while a volcano plot visualizes the complete set of differentially expressed ARGs (Fig. 2B). Genes that were significantly upregulated included IGF1, a potent survival factor, and FASN (Fatty Acid Synthase), which is crucial for metabolic reprogramming. This is consistent with their known roles in promoting cancer cell survival and resistance to anoikis [18, 19]. Conversely, genes like GKN1, a putative tumor suppressor, were downregulated.

Fig. 1.

Schematic representation of the overall design and step-by-step workflow of this investigation.

Fig. 2.

Identification and characterization of prognostic anoikis-related genes (ARGs) in bladder cancer. (A) Heatmap displaying the expression patterns of the top 50 differentially expressed ARGs (25 up-regulated, 25 down-regulated) between BC and normal tissues. The selection criteria were FDR 0.05 and $|$log2(Fold Change)$|$ $>$1. (B) Volcano plot visualizing the distribution of all differentially expressed ARGs. (C) Forest plot of the results from univariate Cox regression analysis, identifying 54 ARGs that were significantly associated with patient prognosis in BC (p 0.05). (D) To explore the interplay among key prognostic genes, a network plot was constructed, illustrating the co-expression relationships between the 54 identified ARGs. ARGs, anoikis-related genes; BC, bladder cancer.

Univariate Cox regression analysis was performed on the 136 ARGs to evaluate their prognostic significance. This identified 54 genes that were significantly correlated with survival outcomes in the BC cohort (p 0.05), as shown in a forest plot (Fig. 2C) and listed in Supplementary Table 3. Furthermore, we constructed a network diagram to visualize co-expression relationships among the most prominent prognosis-related ARGs (Fig. 2D). We also investigated the genomic alterations of these 54 prognostic ARGs, specifically examining their copy number variation (CNV) frequencies and chromosomal locations (Fig. 3A,B).

Fig. 3.

Construction and validation of two distinct molecular subtypes based on prognostic ARG expression. (A,B) Genomic landscape of the 54 prognostic ARGs, showing their copy number variation (CNV) frequency (A) and chromosomal locations (B). (C) Consensus clustering matrix for the prognostic ARGs, indicating that k = 2 is the optimal number for classifying the cohort into two distinct subtypes. (D) Kaplan-Meier survival curves comparing the overall survival between the two ARG-defined subtypes (ARGcluster A and B), demonstrating a significant difference in patient survival (p 0.001). (E–G) Validation of the subtype classification using PCA (E), t-SNE (F), and UMAP (G), all of which confirmed clear separation of the two clusters. PCA, Principal Component Analysis; t-SNE, t-distributed Stochastic Neighbor Embedding; UMAP, Uniform Manifold Approximation and Projection.

To stratify patients based on the expression profiles of the 54 prognostic ARGs, consensus clustering was performed using the “ConsensusClusterPlus” R package. This analysis optimally partitioned the BC cohort into two distinct molecular subtypes (k = 2), hereafter referred to as ARGcluster A and B (Fig. 3C; Supplementary Table 4). Kaplan-Meier analysis revealed a significant difference in OS between these groups (p 0.001), with ARGcluster B exhibiting a markedly unfavorable prognosis (Fig. 3D). The robustness and distinctness of this classification were further validated by Principal Component Analysis (PCA), t-SNE, and UMAP, all of which demonstrated clear separation between the two subtypes (Fig. 3E–G).

We next explored the clinicopathological and molecular characteristics of these ARG-defined subtypes. A heatmap illustrates the differential expression of the 54 ARGs across the two clusters, alongside associated clinical features (Fig. 4A). Boxplots further visualize these expression differences (Fig. 4B). Moreover, analysis of the tumor microenvironment indicated different immune infiltration patterns between the subtypes (Fig. 4C).

Fig. 4.

Functional and immunological characterization of the ARG-defined molecular subtypes. (A) A heatmap was generated to visualize the differential expression patterns of the 54 ARGs across the two identified molecular subtypes, annotated with corresponding clinicopathological characteristics. (B) Boxplots comparing the expression levels of the prognostic ARGs between the two molecular subtypes. (C) Comparison of the relative abundance of infiltrating immune cell populations in the two subtype clusters. (D) Results of GSVA showing significantly enriched KEGG pathways in ARGcluster B, including “focal adhesion”, “regulation of actin cytoskeleton”, and “ECM receptor interaction”. (E) GSEA plot depicting the top five most significantly enriched pathways associated with ARGcluster B (*p 0.05; ***p 0.001). GSVA, Gene Set Variation Analysis; GSEA, Gene Set Enrichment Analysis; ECM, extracellular matrix.

To elucidate the underlying biological functions driving these differences, we conducted Gene Set Variation Analysis (GSVA) for KEGG pathways. This revealed that ARGcluster B, the group with worse outcome, was significantly enriched in pathways critical for tumor invasion and metastasis, such as “focal adhesion”, “regulation of actin cytoskeleton”, and “ECM receptor interaction” (Fig. 4D). These pathways are central to cell-matrix adhesion and migration, and their upregulation strongly suggests that cells in this subtype have enhanced capabilities for resisting anoikis and invading surrounding tissues. GSEA further corroborated these findings, identifying “CHEMOKINE_SIGNALING_PATHWAY” and “CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION” among the top five enriched pathways in ARGcluster B (Fig. 4E; Supplementary Table 5). This molecular profile provides a strong mechanistic rationale for the poor prognosis observed in this patient subgroup.

To construct a robust prognostic model, the entire patient cohort was randomly partitioned into a training set and a test set. LASSO Cox regression analysis of the training dataset identified the most effective prognostic genes. This process yielded a 10-gene signature, which we named the ARGmodel, consisting of IGF1, CALR, E2F1, MYC, PLK1, SATB1, FASN, ID2, RAC3, and GKN1 (Fig. 5A,B). For each patient, a risk score was computed using the expression levels of these signature genes weighted by their corresponding coefficients (Supplementary Table 6).

Fig. 5.

Development and validation of a prognostic signature based on anoikis-related genes (ARGmodel). (A,B) LASSO regression analysis for the selection of optimal prognostic genes. (A) Tuning parameter (lambda) selection via cross-validation. (B) LASSO coefficient profiles of the selected ARGs. (C–E) Kaplan-Meier survival curves comparing overall survival between two groups stratified by the ARGmodel in the training set (C), test set (D), and entire cohort (E). (F–H) Time-dependent Receiver Operating Characteristic (ROC) curves assessing the predictive accuracy of the ARGmodel for 1-, 3-, and 5-year overall survival in the training set (F), test set (G), and entire cohort (H). AUC values are indicated in the plots. (I) Forest plot of results from multivariate Cox regression analysis evaluating the independent prognostic value of the ARGmodel risk score. *p 0.05; **p 0.01; ***p 0.001.

We then rigorously validated the prognostic performance of the ARGmodel. Kaplan-Meier analysis revealed the model successfully stratified patients into high- and low-risk categories with markedly different OS. This significant stratification was observed in the training set (p 0.001), and subsequently confirmed in both the test set (p 0.001) and the entire cohort (p 0.001) (Fig. 5C–E). As expected, individuals in the high-risk group demonstrated significantly worse survival. Furthermore, ROC curve analysis yielded high Area Under the Curve (AUC) values for 1-, 3-, and 5-year OS predictions, confirming the superior predictive accuracy of the model (Fig. 5F–H). Crucially, multivariate Cox regression analysis established that the ARGmodel risk score was an independent prognostic factor for BC (p = 0.0032), distinct from other clinical variables (Fig. 5I).

We then visualized the characteristics of this signature. A heatmap illustrates the expression profiles of the 10 hub ARGs across patients sorted by risk score (Fig. 6A). The risk scores were also significantly different between the previously identified molecular subtypes (Fig. 6B). Finally, an alluvial diagram was generated to visualize the dynamic relationships between the molecular clusters, the ARGmodel risk stratification, and patient survival status (Fig. 6C).

Fig. 6.

Evaluation of the ARGmodel and construction of a clinically applicable prognostic nomogram. (A) Heatmap illustrating the expression profiles of the 10 model ARGs in patients stratified by risk score. (B) Comparison of risk scores between the two previously identified molecular subtypes (ARGcluster A and B). (C) Alluvial diagram visualizing the relationship between molecular subtypes, ARGmodel risk groups, and patient survival status. (D) A prognostic nomogram developed by integrating the ARGmodel risk score with clinicopathological features. (E) Cumulative hazard curves for patients stratified by the nomogram-derived risk (nomoRisk). (F) Calibration plot assessing concordance between the nomogram’s predicted survival probabilities and actual observations. (G–I) Decision Curve Analysis (DCA) evaluating the clinical net benefit of the nomogram for predicting 1- (G), 3- (H), and 5-year (I) survival. *p 0.05; **p 0.01; ***p 0.001.

We constructed a nomogram that incorporates the ARGmodel risk score, along with significant clinicopathological factors, to improve the clinical applicability of our prognostic signature (Fig. 6D). Patients stratified by the nomogram-derived risk score (nomoRisk) showed distinct outcomes, with the high-risk group exhibiting a significantly higher cumulative hazard (Fig. 6E). The performance of the nomogram was rigorously evaluated. Calibration plots for 1-, 3-, and 5-year survival demonstrated excellent agreement between the nomogram-predicted probabilities and the actual observed outcomes, indicating high predictive accuracy (Fig. 6F). Furthermore, we performed Decision Curve Analysis (DCA) to assess its clinical value. For the prediction of 1-, 3-, and 5-year survival, the nomogram provided a greater net benefit across a wide range of threshold probabilities compared to treating all or no patients, confirming its potential as a valuable tool for clinical decision-making (Fig. 6G–I).

To elucidate the biological mechanisms underlying the prognostic significance of our ARG-based signature, we investigated the tumor immune microenvironment (TME). Deconvolution of immune cell fractions revealed significant differences between the high- and low-risk groups. Notably, the high-risk group exhibited significantly lower infiltration of CD8+ T cells (p = 0.043), alongside higher infiltration of M0 macrophages (p = 0.011) and neutrophils (p = 0.039) (Supplementary Fig. 1A). The depletion of cytotoxic CD8+ T cells in high-risk patients suggests a suppressed anti-tumor immune status, potentially leading to poorer response to immunotherapy. Direct correlation analysis further confirmed that the riskScore was positively associated with abundance of M0 macrophages (p = 0.0017) and neutrophils (p = 0.036) (Supplementary Fig. 1D,E). Notably, the high-risk group, which corresponds to high MYC expression, exhibited a distinct immune profile. Significant positive correlations were observed with pro-inflammatory cells such as neutrophils, activated memory CD4+ T cells, and activated mast cells (Supplementary Fig. 1B,C). This suggests the presence of an active, yet potentially ineffective, inflammatory response. Intriguingly, the high-risk group showed significant negative correlations with regulatory T cells (Tregs) and naive B cells (Supplementary Fig. 1B,C). The reduction in Tregs is counterintuitive to typical immunosuppressive models, but may indicate a state of immune chaos or exhaustion where standard regulatory feedback loops are disrupted. Differences in ESTIMATE scores also supported these findings (Supplementary Fig. 1F). A correlation analysis between MYC and CD274 (the gene encoding PD-L1) expression in the TCGA-BLCA cohort revealed a significant positive correlation (Supplementary Fig. 1G) (Pearson correlation, R = 0.3, p = 8.5 $\times$ 10^-10).

To pinpoint the cellular origins of the signature ARGs, we leveraged the public single-cell RNA-sequencing dataset BLCA_GSE145281 from the TISCH database. This dataset comprises 17 distinct cell clusters within the BC TME (Supplementary Fig. 2A). At single-cell resolution, the expression levels of 7 of the 10 signature genes were successfully mapped. This analysis revealed highly specific expression patterns. For instance, CALR was predominantly expressed in B cells and NK cells, while MYC expression was concentrated in B cells and CD4+ T cells (Supplementary Fig. 2B,C). These findings provide granular detail on which cell types contribute to the ARG signature. The expression of MYC within CD4+ T cells is particularly intriguing, as it may drive their differentiation towards pro-tumorigenic subsets like Tregs or Th17 cells, further contributing to the immunosuppressive microenvironment [20].

Finally, we explored the potential of the ARGmodel to predict therapeutic response in BC. By estimating the sensitivity to a panel of clinical drugs, distinct response patterns were observed between the high- and low-risk groups (Fig. 7). This analysis indicates the risk stratification system based on our 10-ARG signature (IGF1, CALR, E2F1, MYC, PLK1, SATB1, FASN, ID2, RAC3, and GKN1) could potentially guide treatment selection. The pathways regulated by these genes appear to be linked to differential drug sensitivity, marking them as promising targets for the development of personalized therapeutic strategies in BC. We posit that a subset of ARGs may act as valuable biomarkers for OS in BC, and as actionable targets for therapeutic intervention.

Fig. 7.

Differential drug sensitivity analysis between high- and low-risk groups. (A–Z) Response patterns to different clinical drugs between high-risk and low-risk groups.

We next employed a multi-step approach to identify a key driver in our prognostic signature. First, a PPI network was constructed using the STRING database to investigate functional connections among the 54 differentially expressed ARGs with prognostic significance. Subsequent analysis with Cytoscape identified MYC as the principal hub gene within this network based on its high degree of connectivity (Fig. 8A). This network centrality suggests a pivotal regulatory role for MYC. The top 10 hub genes ranked by score are detailed in Supplementary Table 7.

Fig. 8.

Identification of MYC as a key hub gene and validation of its prognostic and clinical significance. (A) PPI network analysis identified MYC as a central hub gene. (B) MYC expression was significantly elevated in tumor tissues compared to normal tissues in both the TCGA-BLCA and GSE13507 cohorts. (C,D) Kaplan-Meier analysis demonstrated that high MYC expression was consistently correlated with poorer overall survival in public cohorts (C). This was confirmed in our independent GXMU-BC cohort (D). (E,F,J,K) MYC expression was significantly higher in MIBC than in non-MIBC, a finding observed at both the transcriptomic (E,F) and protein (J,K) levels. (G) In the IMvigor210 cohort, high MYC expression was associated with a diminished response to immunotherapy (Response: CR/PR; Non-response: SD/PD). (H,I) Representative immunohistochemistry images confirm upregulation of MYC protein expression in GXMU-BC tumor tissues. Scale bars: 50 µm (20$\times$); 20 µm (40$\times$). PPI, protein-protein interaction; NMIBC, non-muscle invasive bladder cancer; MIBC, muscle-invasive bladder cancer.

We proceeded to validate the clinical and prognostic relevance of MYC using the TCGA-BLCA and GSE13507 cohorts, confirming its unique position among the signature genes. MYC emerged as a pivotal factor among the identified hub genes. In the TCGA cohort, MYC was significantly upregulated in BC tissue compared to normal controls (p = 0.016, Fig. 8B), and its elevated expression was strongly predictive of unfavorable OS (p = 1.76 $\times$ 10^-4, Fig. 8C). The prognostic value of MYC was further validated in our independent GXMU-BC cohort, which confirmed its negative impact on OS (p = 0.032, Fig. 8D).

Crucially, MYC was linked to key features of clinical aggressiveness, with high MYC expression being associated with advanced tumor stage and MIBC (Fig. 8E,F). Furthermore, analysis of the IMvigor210 immunotherapy cohort revealed that high MYC expression correlated with a significantly poorer response to immune checkpoint blockade (p = 0.03), suggesting a role in immunotherapy resistance (Fig. 8G). This combination of network centrality, robust prognostic power, and strong association with both tumor invasion and immunotherapy failure provided a compelling rationale for focusing on MYC as a key regulator.

At the protein level, IHC analysis of our 81 GXMU-BC patient samples confirmed aberrant MYC expression in 74% (60/81) of BC tissues (Fig. 8H–K). To further assess MYC as a predictive biomarker, we investigated its association with outcomes from ICI therapy. In two separate cohorts of 43 and 38 ICI-treated BC patients, a significant link was observed between MYC expression and clinical response. Specifically, patients with high MYC expression were significantly less likely to achieve a clinical response (Complete/Partial Response) compared to those with low MYC expression (p = 0.044 and p = 0.036 for the two cohorts, Fisher’s exact test). This finding indicates that MYC may be a potential biomarker for identifying patients who are resistant to ICI therapy (Tables 2,3).