Scientific Papers

Multi-cohort validation of Ascore: an anoikis-based prognostic signature for predicting disease progression and immunotherapy response in bladder cancer | Molecular Cancer

Description of Image

Mutation patterns and expression of ARGs in BLCA

In this study, we identified a total of 332 ARGs and analyzed their mutation prevalence across 398 BLCA patients. Notably, mutations in ARGs were detected in 385 of these individuals. Fig. S1A presents the top 10 mutated ARGs, signifying their pervasive role in BLCA pathogenesis. Particularly, TP53, PIK3CA, and RB1 emerged as the most frequently mutated ARGs. To delve into the implications of these mutations, we partitioned the BLCA patients into two distinct groups: the Wild and Mutant groups. GSVA was utilized to calculate the enrichment score of KEGG pathways in individual BLCA patients, and the ‘limma’ tool was employed to perform differential analysis. Results show that pathways upregulated in the Mutant group, including “Cell Cycle”, “DNA Replication”, and “Mismatch Repair” (Fig. S1B). Our data indicate potential defects in DNA repair, cell cycle regulation, and genomic stability associated with ARG mutations. Patients with mutations in ARGs may exhibit more aggressive tumor behavior.

To gain further insights into the involvement of ARGs in BLCA, we compared the 332 ARGs expression between normal and tumor bladder tissues using TCGA cohort. Our analysis revealed 112 genes upregulated in normal tissues and 86 in tumor tissues (Table S1). Univariate Cox regression yielded 54 prognosis-related ARGs (Fig. S1C), categorized as either ‘risk’ or ‘protective’ genes. The relationship between these 54 ARGs is visualized in Fig. S1D, showing the correspondent expression between 54 prognosis-related ARGs.

Consensus clustering on ARGs in BLCA

To underscore the clinical relevance of ARGs, consensus clustering was conducted using the “ConsensusClusterPlus” R package, revealing two distinct patient clusters when k = 2 (Fig. 2A, Fig. S2A, B). Kaplan–Meier survival analysis confirmed significant survival differences between the clusters (P = 0.005; Fig. 2B). Clinical variables like age, pathological stage, and lymph node metastasis also exhibited pronounced variance between the two clusters (Table 1). Patients in Cluster1 were older, had higher stage grading, higher T grading, and were more likely to have lymph node metastasis. Importantly, the two clusters displayed divergent expression patterns of the 54 prognosis-related ARGs (Fig. 2C), with “risk” genes highly expressed in Cluster1, whereas “protective” genes in Cluster2.

Fig. 2
figure 2

Consensus Clustering based on Prognostic Anoikis-Related Genes (ARGs) in BLCA. A Consensus matrix depicting the clustering results when k (cluster number) is set to 2. B Kaplan–Meier curves illustrating the overall survival differences between the two clusters (P = 0.005). C Heatmap displaying the expression levels of 54 prognostic ARGs, along with clinical characteristic annotations for each cluster. D Volcano plot of differentially expressed genes (DEGs) between clusters, with Cluster 2 as control (| logFC |> 1, P < 0.05). E Gene Ontology (GO) analysis highlighting the biological processes (BP), cellular components (CC), and molecular functions (MF) enriched between Cluster 1 and Cluster 2

Table 1 Clinical characteristics of two clusters

Next, in order to gain insights into the molecular characteristics behind the distinction, 1520 differentially expressed genes (DEGs) between the clusters were identified (Fig. 2D). Functional analysis through GO and KEGG analysis (Fig. 2E, Fig. S2C) showed that Cluster1 was closely related to pro-invasive functions such as “positive regulation of cell activation” and “cytokine-cytokine receptor interaction”. On the other hand, Cluster2 may be associated with metabolic alterations. Additionally, GSVA revealed that Cluster1 was closely related to epithelial‐to‐mesenchymal transition (EMT) and nuclear factor‐κB (NF‐κB) signaling (Fig. S2D), both of which were linked to resistance to anoikis or apoptosis [36, 37]. These findings indicated significant differences in biological functions between the two clusters categorized by ARGs and demonstrated the rationality and implications of such categorization in bladder cancer.

Exploring immune infiltration and responsiveness to immunotherapy between two clusters

The Tumor immune microenvironment (TIME) plays a critical role in cancer development and progression [38]. Variability in TIME among patients is closely linked to responsiveness to immunotherapy [39]. Despite limited research exploring the relationship between anoikis and TIME, understanding the TIME landscape in relation to different ARGs expression patterns in BLCA is crucial. Utilizing three algorithms: “ESTIMATE”, “CIBERSORT”, and “ssGSEA”, we depicted variations in TIME components between clusters (Fig. 3A). Cluster 1 displayed elevated Stromal, Immune, and ESTIMATE scores, yet had reduced Tumor Purity, presenting a more complex immune microenvironment characterized by increased tumor heterogeneity. Moreover, Cluster1 showed a higher abundance of immunosuppressive cells such as regulatory T cells and CD4 T cells, while Cluster2 exhibited higher expression of CD8 T cells, suggesting an active and potentially effective anti-tumor immune response. These findings indicate significant differences in the TIME between the clusters, with Cluster 1 being more complex.

Fig. 3
figure 3

Immune Infiltration and Responsiveness to Immunotherapy across Clusters. A ESTIMATE scores and immune cell populations compared between Cluster 1 and Cluster 2. (B-E) TIDE analysis including TIDE score B Exclusion score Dysfunction score D and potential immunotherapy responders E (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001)

TIDE is an algorithm designed to predict the response of immunotherapy, with a higher TIDE score meaning a reduced benefit from immunotherapy and an increased risk of immune escape. In our analysis, Cluster1 exhibited higher “TIDE scores”, “Exclusion scores”, and “Dysfunction scores” (Fig. 3B-D). Furthermore, Cluster2 had a higher proportion of patients predicted to respond to immunotherapy (P < 0.0001; Fig. 3E). These results support the notion that patients in Cluster2 may derive more benefit from immunotherapy compared to those in Cluster1.

Establishment and validation of an anoikis-based signature

Recognizing the impact of ARGs on BLCA patient outcomes and immunotherapeutic responsiveness, a prognostic signature was pursued to understand bladder cancer’s underlying complexities. Through univariate Cox regression analysis, 654 of the previously mentioned 1520 DEGs were identified as preliminary prognostic, which was narrowed down using the LASSO algorithm, leading to 16 genes (Fig. S3A, B). After subsequent multivariate Cox regression analysis (Fig. S3C), a prognostic signature called “Ascore” was constructed from four genes: CERCAM, EMP1, GNLY, and PTPRR (Fig. 4A). The formula is as follows: Ascore = (0.218 × CERCAM expression) + (0.291 × EMP1 expression) + (-0.259 × GNLY expression) + (-0.192 × PTPRR expression).

Fig. 4
figure 4

Establishment and Validation of the Ascore Prognostic Signature. A Multivariate Cox coefficients for four ARGs (CERCAM, EMP1, GNLY, PTPRR) in the prognostic signature. B Ascore distribution among BLCA patients, sorted from lowest to highest. C Survival status categorized by Ascore for each BLCA patient. D Heatmap displaying expression levels of four genes in different Ascore groups. E Sankey diagram correlating clusters, Ascore groups, and BLCA survival status. F Kaplan–Meier analysis comparing overall survival between high and low Ascore groups in BLCA (P < 0.0001). G Receiver Operating Characteristic (ROC) curves depicting Ascore signature’s predictive performance for 1, 3, and 5-year overall survival in BLCA, with the Area Under the Curve (AUC) values of 0.709, 0.724, and 0.745, respectively. (H–K) Kaplan–Meier analysis and time-dependent ROC curves in two external validation sets: GSE32548 and GSE32894

Next, we calculated Ascore for each individual patient in the BLCA cohort and sorted them from lowest to highest (Fig. 4B). Based on the median Ascore value, patients were divided into high and low Ascore groups. The survival status of the patients (Fig. 4C) and the expression levels of the 4 genes in the signature (Fig. 4D) were also displayed. Notably, patients with higher Ascore values correspond to worse survival outcomes, accompanied by higher expression levels of the “risk” genes, CERCAM and EMP1. To visually depict the relationship between clusters, Ascore subgroups, and survival status in BLCA, we utilized a Sankey diagram (Fig. 4E). The diagram revealed that patients with higher Ascore were more likely to be categorized in the Cluster1 and had a poorer prognosis. This observation was further supported by Kaplan–Meier analysis (Fig. 4F), with patients in the low Ascore group having significantly superior overall survival (OS, P < 0.0001). Similar trends were observed in other indicators such as Disease-Free Survival (DSS, P < 0.0001), Progression-Free Interval (PFI, P < 0.0001), and Disease-Free Interval (DFI, P = 0.0453) in BLCA (Fig. S3D-F).

To evaluate the predictive ability of our prognostic signature, we generated Receiver Operating Characteristic (ROC) curves for OS at 1, 3, and 5 years (Fig. 4G). The area under the curve (AUC) values were 0.709, 0.724, and 0.745, respectively, indicating the good predictive performance of our model, especially in long-term survival. Furthermore, we validated the accuracy of our signature in two independent external validation sets (GSE32548, GSE32894) and yielded satisfying results, with 5-years AUC values of 0.726 and 0.806, respectively (Fig. 4H-K).

To improve clinical utility, we also created a nomogram using Ascore and clinical parameters. Univariate and multivariate Cox regression analyses assessed the effects of clinical features and Ascore on BLCA survival. Both age and Ascore emerged as significant survival predictors (Fig. S4A). Using these predictors, we devised a nomogram estimating 1, 3, and 5-year survival probabilities for BLCA patients (Fig. S4B). Calibration curves attested to the nomogram’s predictive accuracy (Fig. S4C-E). Decision curve analysis (DCA) affirmed its potential patient benefits (Fig. S4F). Overall, these suggest that the Ascore-based nomogram, rooted in anoikis concepts, has substantial clinical predictive value for bladder cancer.

High Ascore indicates advanced disease and anoikis resistance

We further analyzed the relationship between Ascore and clinical features and found older patients as well as those with advanced disease often had higher Ascore values. (Fig. S5A). To determine the biological significance of Ascore, we screened out DEGs between Ascore groups (Fig. S5B) and performed KEGG analysis. The results showed that PI3K-AKT signaling pathway, which was reported to promote proliferation and inhibit anoikis in bladder cancer, was enriched in the high Ascore group (Fig. S5C). While the low Ascore group exhibited opposite with an abundance in the PPAR signaling pathway [40, 41]. Furthermore, anoikis-related gene sets from the GO database further validate our findings (Fig. S5D), with Ascore positively correlated with the “negative regulation of anoikis” while negatively correlated with the “positive regulation of anoikis”. Our findings illustrated that a high Ascore represents an anoikis-resistant status.

We further investigated using single-cell RNA sequencing data from HRA000212, which integrated 8 samples of bladder cancer for subsequent analysis (Fig. S6A). After manually eliminating doublet, a total of 41,387 cells were annotated into seven groups, as shown in the t-SNE and UMAP plot (Fig. 5A, B, Fig. S6B, C): Epithelial (EPCAM); Endothelial (PECAM1); iCAFs (COL1A1, PDGFRA); mCAFs (COL1A1, RGS5); Myeloid (LYZ); B cells (CD79A); T cells (CD3D).

Fig. 5
figure 5

Single-Cell RNA Sequencing Analysis of Ascore Distribution and Biological Significance in Bladder Cancer. A t-SNE plot showing seven main cell types distribution in the integrated dataset, with doublets manually annotated. B Dot plot of marker genes’ expression levels in each cell type. C Ascore and four genes (CERCAM, EMP1, GNLY, PTPRR) expression and distribution across cell types. D t-SNE plot showing Ascore expression levels and patterns in each cell type. E Left plot: Six main epithelial cell subgroups visualized by t-SNE dimensionality reduction. Right plot: Ascore distribution and expression in epithelial cells, highlighting Subgroup0 and Subgroup2. F GO analysis of biological function differences between Subgroup0 and Subgroup2

We proceeded to explore the distribution of four genes comprising the Ascore in the dataset (Fig. 5C, Fig. S6D). Our findings revealed that CERCAM was predominantly expressed in inflammatory cancer-associated fibroblasts (iCAFs), a type of cancer-associated fibroblasts known for their strong pro-proliferation properties. On the other hand, EMP1 exhibited a wide distribution in bladder cancers but showed a particular concentration in endothelial cells. GNLY and PTPRR were primarily found in T cells and epithelial cells, respectively. Subsequently, we calculated the Ascore in each individual cell and observed endothelial cells and iCAFs exhibited higher levels. However, the distribution of Ascore varied significantly in epithelial cells (Fig. 5C, D, Fig. S6D), which sparked our curiosity.

We further sub-clustered epithelial cells into 6 groups (Fig. 5E), and observed high Ascore concentration in the Subgroup0, in contrast to the Subgroup2. Further biological function analysis elucidated the underlying differences. As shown in Fig. 5F, Subgroup0 epithelial cells exhibited a higher propensity for proliferation and migration, and displayed resistance to apoptosis. Moreover, these cells appeared to be involved in suppressing immune responses, particularly T cells activation. On the other hand, in Subgroup2, cells appeared to be in a hypoxic tumor microenvironment.

Collectively, our findings based on both RNA-seq and scRNA-seq proved that higher Ascore correlated with more aggressive and unfavorable phenotypes, and reflected anoikis resistance in bladder cancer.

Immune landscape variations between different Ascore groups in TCGA-BLCA

Since Subgroup0 has shown an association with immunosuppression. We further explored the relationship between Ascore and TIME in the TCGA-BLCA cohort, finding immunosuppressive cells like CD4 + T cells dominated in the high Ascore group (Fig. S7A). On the contrary, CD8 + T cells and NK cells expression was abundant in the low Ascore group, indicating an immune active microenvironment.

TMB is a metric reflecting the number of mutations in cancer. Many researches have shown that high TMB predicts a more likely T-cell response and a greater likelihood of benefiting from treatment with ICIs [42]. In our study, high Ascore patients displayed lower TMB (P = 0.0034; Fig. S7B), indicating less potential benefit from immunotherapy. Additionally, the TIDE results substantiate this finding, with patients in the high Ascore group presenting higher TIDE scores (P < 0.0001; Fig. S7C) and a lower proportion of response (P < 0.0001, Fig. S7D).

Evaluating Ascore as a prognostic tool and predictor of immunotherapy response

To comprehensively assess the relationship between Ascore and immunotherapy, we turned to the IMvigor210 cohort, a well-known large phase II clinical trial investigating treatment response to anti-PD-L1 (atezolizumab) immunotherapy in patients with advanced urothelial cancer. Here, we noticed that a higher Ascore was associated with adverse survival outcomes post-immunotherapy (P = 0.0189; Fig. 6A). We incorporated additional risk factors, such as gender, Eastern Cooperative Oncology Group (ECOG) performance status, liver metastases presence, and smoking history, into our prognostic evaluation. Both univariate and multivariate analyses positioned Ascore as an independent prognostic indicator (Fig. S8A) with robust predictive performance for survival (AUC = 0.673; Fig. 6B), outperforming other known prognostic factors. Intriguingly, this trend remained consistent when we broadened our analysis to include all patients with urothelial cancer receiving immunotherapy within the cohort (Fig. S8B-D).

Fig. 6
figure 6

Ascore’s Prognostic and Predictive Role in Immunotherapy and Chemotherapy Cohorts. A Survival outcomes within post-immunotherapy bladder cancer patients relative to different Ascore groups. B ROC analysis of Ascore’s predictive performance for survival against other prognostic factors like ECOG and liver metastasis in bladder cancer. C Response rates to immunotherapy in bladder cancer patients based on Ascore groups. D ROC analysis illustrating Ascore’s predictive accuracy for immunotherapy response in bladder cancer. E Ascore’s prognostic significance in chemotherapy-treated MIBC patients. F Assessment of Ascore in predicting pathological response to chemotherapy. (ECOG: Eastern Cooperative Oncology Group; ICIs: Immune Checkpoint Inhibitors; PD: progressive disease; SD: stable disease; PR: partial response; CR: complete response; TMB: Tumor Mutational Burden; IC: Immune Cell)

We also noticed that divergent Ascore expressions were associated with varying responses to immunotherapy, with a significant proportion of low Ascore bladder cancer patients responding favorably (P = 0.002; Fig. 6C). We then evaluated the ability of Ascore to predict immunotherapy response in bladder cancer, and the ROC analysis (Fig. 6D) illustrated Ascore’s superior predictive capacity (AUC = 0.717) compared to TMB (AUC = 0.700), PD-L1 expression levels on immune cells (AUC = 0.628), and PD-L1 expression in tumor tissues (AUC = 0.582). Similar trends were found in the entire urothelial cancer cohort (Fig. S8E, F). These findings suggested Ascore’s utility as both a prognostic factor and a potential predictor of immunotherapy response.

Furthermore, to discern whether Ascore’s predictive ability was exclusive to immunotherapy, we included a group of 149 MIBC patients who were administered preoperative platinum-based neoadjuvant chemotherapy. We found that Ascore, while prognostically significant (P = 0.0213; Fig. 6E), did not anticipate the pathological response to chemotherapy (P = 0.221; Fig. 6F). This finding provides unique insights into Ascore’s specificity, implying a more pronounced predictive capability in the realm of immunotherapy response.

Validation of Ascore as a non-invasive prognostic biomarker in circulating tumor cells

To evaluate the real-world utility of Ascore in prognosticating survival and immunotherapy response in patients with bladder cancer, we conducted a retrospective investigation at Drum Tower Hospital, affiliated with Nanjing University in China. Patients who were diagnosed with bladder cancer and underwent RC were included in our analysis. These patients were stratified into two distinct cohorts, termed Gulou-Cohort1 and Gulou-Cohort2, as depicted in Fig. 7A.

Fig. 7
figure 7

Study Design and Ascore Prognostic Validity Evaluation in Gulou-Cohort1. A Retrospective study stratification flowchart showing Gulou-Cohort1 (134 patients) and Gulou-Cohort2 (40 patients) compositions. Circulating tumor cells were successfully isolated and quantified from 62 Gulou-Cohort1 patients’ blood samples. B Kaplan–Meier curves comparing prognostic outcomes between the high and low Ascore groups in Gulou-Cohort1. C ROC analysis showing Ascore’s predictive value versus CTC counts and ECOG score in Gulou-Cohort1. (CTC: Circulating tumor cell; ECOG: Eastern Cooperative Oncology Group)

Circulating tumor cells (CTCs) are neoplastic cells that have disseminated from the primary tumor site into the circulatory system and are implicated as a prognostic indicator in cancer progression and metastasis [43]. Within Gulou-Cohort1, we enrolled 134 patients who had not received any form of systemic therapy prior to RC. Blood samples were prospectively obtained, leading to the successful isolation and characterization of CTCs in 62 cases (Fig. 7A). Relative RNA expression of four specific genes (CERCAM, EMP1, GNLY, PTPRR) was quantified to compute the Ascore for each patient, subsequently categorizing them into Ascore-High and Ascore-Low groups (Table S2). Importantly, the Ascore-High group was predominantly comprised of older individuals, exhibited more advanced pathological stages, and had a greater likelihood of lymph node metastases (Table 2), but there was no significant association with ECOG score. These observations are congruent with our prior research. Moreover, higher Ascore values in CTCs were associated with worse prognostic outcomes compared to the Ascore-Low group (P = 0.0364; Fig. 7B).

Table 2 Clinical characteristics of high and low Ascore groups in Gulou-Cohort1

To further substantiate the prognostic significance of Ascore in CTCs, we included clinical parameters from Table 2 in a univariate Cox regression analysis and also integrated CTC count for a comparative assessment (Fig. S9). The CTC count, previously established as a significant prognostic marker in various cancers such as prostate cancer and small-cell lung cancer [44, 45], predicts disease outcomes effectively. The univariate Cox regression analysis identified lymph node metastasis, ECOG score, CTC count, and Ascore as risk factors for Gulou-Cohort1 (P < 0.05). Subsequent multivariate regression analysis revealed that both Ascore and ECOG score were independent prognostic factors (Fig. S9). We then evaluated the predictive capabilities of Ascore, CTC count, and ECOG score on the prognosis of Gulou-Cohort1 patients. The ROC curve demonstrated that Ascore, with an AUC of 0.803, outperformed the CTC technique (AUC = 0.735) and ECOG score (AUC = 0.683) in prognostic prediction (Fig. 7C). This validation within our cohort underscores Ascore’s potential as a non-invasive prognostic marker and confirms its clinical applicability.

Ascore predicts anti-PD-1 immunotherapy response in Gulou-Cohort2

Our preceding study proposed that Ascore could predict the responsiveness of anti-PD-L1 immunotherapy in bladder cancer patients. To extend this finding to anti-PD-1 immunotherapy, we analyzed 40 patients who underwent standard anti-PD-1 treatment prior to RC. Pre-treatment tumor tissue was evaluated for the expression of the four aforementioned genes through immunohistochemistry (IHC). Each sample received an H-score, with the scoring criteria for the expression of these genes in tumor cells detailed in Fig. S10A. Individual Ascore values were derived using a previously established formula, as exhibited in Table S3. Figure 8A showed representative IHC images from the Ascore-Low (Patient4) and Ascore-High (Patient7) groups, illustrating a responder and a non-responder, respectively. We found significant differences in pathological response rates between the Ascore-High and Ascore-Low groups (P = 0.004; Table 3). Specifically, lower Ascore values in pre-treatment tumor samples were linked to a more favorable immunotherapy response (45% vs. 12.5%, P < 0.001). In addition, a higher rate of complete response (CR) was observed in the Ascore-Low group compared to the Ascore-High group (25% vs. 7.5%, P = 0.041). Consistently, patients who responded to immunotherapy (CR/PR) had noteworthy lower Ascores (P < 0.0001). Our results revealed a strong relation between Ascore and anti-PD-1 immunotherapy response.

Fig. 8
figure 8

Ascore Predictive Capability for Anti-PD-1 Immunotherapy Response in Gulou-Cohort2. A Representative immunohistochemistry (IHC) images illustrating the expression of four key genes (CERCAM, EMP1, GNLY, PTPRR) in two patients from Gulou-Cohort2 (Scale bars = 100 μm). Patient 4, who responded to anti-PD-1 therapy, had a low Ascore, in contrast to non-responder Patient 7, who had a high Ascore. B Distribution of Ascores among different response groups (CR: complete response; PR: partial response; SD: stable disease; PD: progressive disease; ***P < 0.001). C ROC curves comparing the predictive accuracy of Ascore (AUC = 0.913) versus PD-L1 expression in tumor-infiltrating immune cells (ICs) (AUC = 0.662). D Decision curve analysis (DCA) indicating the net benefit of using Ascore compared to evaluating ICs’ PD-L1 expression. E Kaplan–Meier curves c showing a correlation between higher Ascore values in tissue samples and reduced survival rates (P = 0.0194)

Table 3 Pathologic response of high and low Ascore groups in Gulou-Cohort2

To further refine the predictive accuracy of the Ascore signature, we expanded our analysis to include PD-L1 expression in tumor-infiltrating immune cells (ICs) as shown in Fig. S10B. These ICs—comprising macrophages, dendritic cells, and lymphocytes—are considered instrumental in shaping immunotherapy response in bladder cancer. Fig. S10C presents representative IHC images for Patient4 and Patient7, with their PD-L1 expression on ICs being IC2 and IC3, respectively. This challenges the common belief that higher PD-L1 levels correlate with a successful immunotherapy response. ROC analysis further indicated the predictive capability of Ascore (AUC = 0.913) was superior to that of ICs (AUC = 0.662) in forecasting immunotherapy pathologic response (Fig. 8C). DCA curves suggested that the net benefit of the Ascore was obviously higher than PD-L1 expression in ICs (Fig. 8D).

In alignment with data from Gulou-Cohort1, elevated Ascore values in tissue samples correlated with poorer survival rates, reinforcing the general applicability of Ascore as a predictive biomarker for bladder cancer progression (P = 0.0194, Fig. 8E).

Description of Image

Source link