3.1 CAF scoring and survival analysis in BUC

Using the MCPcounter algorithm, the CAF score for BUC patients was calculated. Subsequently, BUC patients were divided into two groups based on the median CAF score, and survival analysis revealed that patients with lower CAF scores had better prognosis. (Fig. 1A). Meanwhile, similar analysis was conducted in E-MTAB-4321, revealing that patients with lower CAF scores also had better prognosis (Fig. 1B). These results suggest that CAF scoring may be associated with the prognosis of BUC patients. Subsequently, WGCNA analysis was conducted on TCGA data and E-MTAB-4321 data. In TCGA, genes were classified into 15 gene modules based on their expression profiles (Fig. 1C). Since CAFs are typically expressed in the cellular stroma, the correlation between these modules and both CAF scores and StromalScore was analyzed (Fig. 1D). The results showed that the yellow module had high correlation with CAF and StromalScore. Similarly, WGCNA analysis of E-MTAB-4321 revealed that the genes were classified into 15 gene modules (Fig. 1E), and the results showed that the purple module had high correlation with CAF and StromalScore (Fig. 1F).

Fig. 1

Overall survival and WGCNA analysis based on CAF scoring in TCGA and E-MTAB-4321 databases. A Survival analysis in TCGA cohort. B Survival analysis in E-MTAB-4321 cohort. C WGCNA analysis of TCGA cohort. D Correlation analysis between module and CAF scoring. E, F WGCNA analysis of the E-MTAB-432 cohort

3.2 Development and validation of CAF score

Intersect the genes identified in the modules from TCGA and E-MTAB-4321 (Fig. 2A), and 307 genes were finally obtained. Patients in TCGA were randomly divided into training and validation cohorts. LASSO Cox regression analysis was conducted in the training cohort to successfully construct a prognostic model consisting of five genes (Fig. 2B, C). Subsequently, using the coefficients obtained from LASSO Cox regression analysis, risk scores were calculated for each patient. Based on the median risk score, patients in the training cohort were divided into high-risk and low-risk groups. The expression of these five genes shows significant differences between high-risk and low-risk groups in both training and validation cohorts (Fig. 2D, E). The heatmap shows the distribution of 5 genes with the risk score in the training cohort (Fig. 3A). Patients with high risk score have poorer prognosis in the training cohort (Fig. 3D, G, J). The ROC curves showed area under the curve (AUC) values of 0.768, 0.790, and 0.746 for 1 year, 3 years, and 5 years, respectively in the training cohort (Fig. 3M). Similar results were obtained in the validation cohort, indicating that this prognostic model has high predictive accuracy (Fig. 3B, E, H, K, N). To assess the stability of this prognostic model, the same method was applied to calculate risk scores for patients in E-MTAB-4321, and they were then divided into high-risk and low-risk groups. The results similarly indicated that patients with high risk score had poorer prognosis (Fig. 3C, F, I, L). The ROC curves showed the AUC values of 0.763, 0.719, and 0.748 for 1 year, 3 years, and 5 years, respectively (Fig. 3O).

Fig. 2

Construction of the CAF-related gene prognostic signature. A Gene intersections in modules significantly associated with CAFs of TCGA and E-MTAB-4321 cohorts. B The LASSO Cox analysis of CAF-related genes in the training cohort. C The value of the penalty parameter (λ) was determined based on the minimum partial likelihood deviationby tenfold cross-validation. D, E Differential expression of candidate genes between high-risk and low-risk group in training cohort (D) and validation cohort (E)

Fig. 3

Prognostic analysis of 5 genes signature. A-C Heatmap of 5 CAF-related genes in the training (A), test (B) cohort of TCGA cohort and E-MTAB-4321 (C) cohort. D-F The distributions of the risk score in the training (D), test (E) cohort of TCGA cohort and E-MTAB-4321 (F) cohort. G-L The distributions of OS status and OS of patients between high-risk and low-risk groups, patients in the high-risk group had higher score values and mortality in the training (G, J), test (H, K) cohort of TCGA cohort and E-MTAB-4321 (I, L) cohort. M–O ROCs for 1, 3, and 5 year survival time based on the risk score in the training (M), test (N) cohort of TCGA cohort and E-MTAB-4321 (O) cohort

3.3 Independent prognostic analysis.

Univariate and multivariate Cox regression analyses were conducted on clinical characteristics and risk scores within the TCGA cohort, revealing that the risk score is an independent prognostic factor for BUC patients (Fig. 4A, B, p < 0.001). Similar results were also obtained in the E-MTAB-4321 cohort (Fig. 4E, F, p < 0.001). Compared to the use of Age, Gender and Stage to predict the prognosis of BUC, the use of risk score has higher sensitivity and specificity in the TCGA and EMATAB datasets, owing to its’ higher ROC and the C index (Fig. 4C, D, G, H). Furthermore, a nomogram was successfully constructed, and calibration curves demonstrated its high predictive accuracy (Fig. 4I, J).

Fig. 4

Independent prognostic analysis of the signature and construction of nomogram. A, E The univariate Cox regression analysis in the TCGA (A) and E-MTAB-4321 (E) cohort. B, F The multivariate Cox regression analysis in the TCGA (B) and E-MTAB-4321 (F) cohort. C, G ROC curve of the prognosis of BUC patients in the TCGA and E-MTAB-4321 cohort. D, H c-index curve of the prognosis of BUC patients in the TCGA and E-MTAB-4321 cohort. I The nomogram that includes the risk score and clinic characteristics. J The calibration curves for the predicted 1-, 3-, and 5-years OS rates

3.4 Pathway enrichment analysis

Utilizing GO and KEGG analyses to investigate potential key pathways between the two groups. As shown in Fig. 5A, B, GO identified most extracellular structure pathways were enriched, including extracellular matrix, extracellular structure organization, external encapsulating organization, structure organization muscle contraction. There are many enrichment genes between high and low risk groups in the Collagen-Containing Extracellular Matrix gene dataset. And in this genetic concentration, the amount of enrichment genes is relatively higher and the P value is the lowest, indicating that the pathway is very important (Fig. 5B). The KEGG analysis showed that the pathway was significantly enriched in the ECM − receptor interaction, Hypertrophic cardiomyopathy, Protein digestion and absorption, Dilated cardiomyopathy (Fig. 5C, D). At the same time, more genes were enriched in the PI3K-Akt signaling pathway, highlighting the difference in the regulation of this signaling pathway in high and low risk groups (Fig. 5D). In the E-MTAB-4321 cohort, the enriched pathways were also quite similar (Fig. 5E-H).

Fig. 5

Functional enrichment analysis. A, B GO enrichment analysis in the TCGA cohort. C, D KEGG enrichment analysis in the TCGA cohort. E, F GO enrichment analysis in the E-MTAB-4321 cohort. G, H KEGG enrichment analysis in the E-MTAB-4321 cohort

3.5 Immune profiling analysis

One of the hallmarks of cancer progression is the formation of an immunosuppressive and tolerant TME, so it is important to clarify the relationship between the TME and tumor progression. We further investigated the relationship between this prognostic model and TME. The results showed that most immune cells were lower in the low-risk group (Fig. 6A). Subsequently, the expression differences of immune checkpoint-related genes were compared between the two groups. The results showed that most immune checkpoint-related genes had higher expression levels in the high-risk group (Fig. 6B, C). The above results indicate that the low-risk group is closely related to a lower abundance of immune cell infiltration.

Fig. 6

TME and checkpoint related genes analysis. A The comparison of immune cell fractions between groups of the TCGA cohorts via the ssGSEA method. B, C Radar plots show diferences in immune checkpoints between the high- and low-risk group in the TCGA cohort

3.6 TMB analysis

Genomic mutations were compared between the two groups. And ARID1A were the top five mutations genes in both two groups (Fig. 7A, B). TMB analysis indicated that the low-risk group had a higher TMB score (p = 0.04; Fig. 7C). Survival analysis indicates that the patients with lower TMB have poorer prognosis. (Fig. 7D). The combined analysis of TMB and risk score revealed that the patients with low TMB scores and high risk scores had the worst prognosis (Fig. 7E), suggesting that this model combined with TMB score could be used to evaluate the prognosis of BUC patients.

Fig. 7

Tumor mutational burden (TMB) analysis. A, B Top 20 mutated genes shown between the high- and low-risk groups. C Boxplot illustrating TMB differences between the groups. D Survival analysis of overall survival between the TMB subgroups. E TMB score and risk score were combined for survival analysis

3.7 Low-risk group is more sensitive to immunotherapy

To further explore the role of this model in immunotherapy, we calculated the TIDE scores for each patient. The results showed that patients with high risk score had higher TIDE scores, suggesting that these patients may have higher immunotherapy efficacy (Fig. 8A). In addition, low-risk group had lower infiltrating abundance of MDSC and CAF. However, TAM M2 infiltrating abundance was higher (Fig. 8B-D). Further calculations were carried out in E-MTAB-4321 cohort and the results were consistent (Fig. 8E-H).

Fig. 8

Differences in TIDE (A), MDSC (B), CAF (C), and TAM M2 (D) between high- and low-risk groups of TCGA cohort. Differences in TIDE (E), MDSC (F), CAF (G), and TAM M2 (H) between high- and low-risk groups of E-MTAB-4321 cohort

3.8 Drug Susceptibility analysis.

To explore the role of this model in predicting drug therapy, we used the "pRRophetic" package to calculate the IC50 for each drug and compared them between the two groups. The results indicate that patients with high risk scores exhibited lower IC50 values for AKT inhibitor VIII, Pazopanib, Sunitinib, Dasatinib (Fig. 9A–H). However, patients with high risk scores had higher IC50 values for JNK Inhibitor VIII, Mitomycin C (Fig. 9I–L).

Fig. 9

Drug sensitivity analysis. Correlation analysis between risk scores and IC50 levels of A AKT inhibitor VIII, C Pazopanib, E Sunitinib, G Dasatinib I JNK Inhibitor VIII, and K Mitomycin C. Boxplots depict differences in estimated IC50 levels of B AKT inhibitor VIII, D Pazopanib, F Sunitinib, H Dasatinib J JNK Inhibitor VIII, and L Mitomycin C between high- and low-risk groups