Identification of a cancer-associated fibroblast signature for predicting prognosis and immunotherapeutic responses in bladder urothelial carcinoma

Abstract

Background

Cancer-associated fibroblasts (CAFs) are the most important cellular components in bladder urothelial carcinoma (BLCA) and are involved in the development and immunosuppression of BLCA. Therefore, we aimed to construct a CAF-associated signature for predicting the prognosis and immunotherapy response in patients with BLCA.

Methods

CAF infiltration and stromal score were quantified using two algorithms. Weighted gene co-expression network analysis (WGCNA) was performed to identify the CAF-associated modules and hub genes. Univariate Cox and Least Absolute Shrinkage and Selection Operator regression analyses were used for constructing CAF signatures and calculating CAF scores. The ability of the CAF signature to predict prognosis and response to immunotherapy was validated using the data from three cohorts.

Results

WGCNA identified two CAF-associated modules and constructed a CAF signature containing 27 genes. In all three cohorts, patients with high CAF scores had markedly worse prognoses than those with low CAF scores, and CAF scores were independent risk factors. In addition, patients with high CAF scores did not respond to immunotherapy, whereas those with lower CAF scores responded to immunotherapy.

Conclusion

CAF signature can be used to predict prognosis and immunotherapy response to guide individualized treatment planning in patients with BLCA.

Keywords:

• Bladder urothelial carcinoma
• cancer-associated fibroblasts
• WGCNA
• prognosis
• immunotherapy

Introduction

Bladder urothelial carcinoma (BLCA) is the most common tumor of the bladder and a common tumor of the urinary system [1]. Approximately 780,000 bladder cancer survivors are living in the United States [2]. The vast majority of patients with BLCA are men, and the incidence rate in them is approximately three times that of women [2,3]. Miller et al. reported that most of the patients with confirmed BLCA were elderly, and their median age at the time of diagnosis was approximately 73 years [2]. In recent years, the incidence of BLCA is rising in China, and statistical analysis of data from cancer registration centers in 31 provinces across the country revealed 78,100 new cases of BLCA and 32,100 deaths nationwide in 2014. The clinical treatment methods for BLCA include surgery, adjuvant or neoadjuvant chemotherapy, chemotherapy, immunotherapy, and targeted treatment; however, most patients relapse after the treatment [4,5].

In recent years, advances in gene sequencing technology have provided new clues for the clinical application of BLCA [6]. For example, DNA methylation-related biomarkers are used for non-invasive diagnosis of early bladder cancer [7]. Cancer-associated fibroblasts (CAFs) are a major component of the tumor microenvironment (TME). Tumor cells initiate and maintain the activation of CAFs, which promote the proliferation and metastasis of tumor cells [8]. Chen et al. analyzed single-cell transcriptome RNA sequencing data from eight patients with BLCA and found that inflammatory tumor-associated fibroblasts promoted tumor cell proliferation; bladder cancer cells, co-cultured with these fibroblasts, showed a higher proliferation ability [9]. CAF-associated transcription factor RUNX2 is overexpressed in bladder cancer cells; it promotes the migration and invasion of cancer cells by promoting epithelial-to-mesenchymal transition and stimulating the expression of matrix metalloproteinases, which destroy the basement membrane and type IV collagen to remodel the extracellular matrix (ECM) [10].

CAFs are associated with the progression and poor prognosis of various cancers, such as colorectal [11] and prostate cancers [12]. A multicenter study found that a high expression of platelet-derived growth factor in CAFs is an independent risk factor that affects the prognosis of breast cancer [13]. TGF-β, secreted by CAFs, participates in tumor progression by inhibiting immune cell function and creating an immunosuppressive environment [14]. Moreover, the abundance of CAFs in the TME correlates with the efficacy of immunotherapy [15]. CAFs can recruit myeloid cells and induce their differentiation into myeloid-derived suppressor cells by secreting IL-6 and CXCL12, which further modulate the effect of immunotherapy [16]. CAFs can also upregulate the expression of CTLA-4 on CD8 + T cells and affect their intratumoral infiltration, thereby affecting the immunotherapy response in patients [17]. These findings suggest that CAFs may affect the prognosis of BLCA and immunotherapy response in these patients; therefore, we explored how CAF determines the prognosis and immunotherapy response in patients with BLCA.

Weighted gene co-expression network analysis (WGCNA) is used to analyze gene expression patterns in multiple groups of samples. The method classifies genes with highly similar expression patterns into modules and further analyzes the internal relationship between the key genes in modules (or modules) and the clinical characteristics of patients [18]. It is most commonly applied to high-dimensional datasets in genomics. Further, WGCNA is widely used to screen candidate biomarkers or therapeutic targets. Zhen et al. reported CAF-associated biomarkers of gastric cancer using WGCNA and constructed gene signatures associated with prognosis and immunotherapy response [19]. However, few related studies have been conducted on BLCA. Here, we aimed to identify CAF-associated biomarkers and construct gene signatures for evaluating the prognosis of patients with BLCA and their response to immunotherapy.

Materials and methods

Data processing

We downloaded the BLCA-related sequencing data (fragments per kilobase of exon model per million mapped fragments) and clinical information from the TCGA database. We excluded samples with no follow-up days or incomplete clinical information and included the data of 403 patients with BLCA. Besides, we downloaded the raw data of the GSE32894 [20] and GSE13507 [21] datasets and the clinical information of patients from the GEO database. The raw data were normalized using the limma package [22]. We excluded the samples with incomplete prognostic follow-up data and finally included 224 samples in this study. We used the IMvigor210 dataset [23] for external validation, which contains sequencing data and clinical follow-up information of patients with BLCA treated with atilizumab.

CAF infiltration estimation and stromal score calculation

We used the Estimating the Proportions of Immune and Cancer cells (EPIC) algorithm (as a deconvolution algorithm based on the cell type) to calculate the relative content of CAF in the BLCA samples [24]. The EPIC algorithm analyzed the infiltration ratio of eight types of immune cells, including CAFs, based on the expression data [25]. The Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm [26] calculates immune and stromal scores by analyzing specific gene expression characteristics of immune and stromal cells to predict the infiltration of non-tumor cells. This algorithm was used to calculate the score of stromal cells. In addition, we used Tumor Immune Dysfunction and Exclusion (TIDE) score [27] to reflect the response of patients to immunotherapy in the TCGA-BLCA and GSE32894 datasets. A TIDE score <0 indicated a high likelihood of responding to immunotherapy.

Construction of CAF and stromal-related co-expression networks

Genes with the top 5000 median absolute deviations were screened for WGCNA. Hierarchical clustering was performed before constructing the co-expression network, and an appropriate threshold was set to remove outliers and check for missing values. A suitable soft thresholding power (β) was defined to conform the relationship between genes in the gene co-expression network to the scale-free network distribution. We obtained the expression matrix and β-values to calculate the adjacency of expression levels among all genes to draw the adjacency matrix. Next, the adjacency matrix was converted into a topological overlap matrix (TOM), and the corresponding degree of difference (1-TOM) was calculated. Modules were divided according to the standard of a dynamic shear tree, and the minimum capacity of modules was set to 30. Finally, we calculated the module eigengene (ME), module membership (MM), gene significance (GS), and module significance (MS) of the gene module. GS indicated the degree to which a gene was associated with the clinical information, and its absolute value was directly proportional to the significance of the gene. MS was the average GS of all genes in the module, and a higher value indicated a stronger association between the module and clinical information; this value was used to select a significant module for subsequent analysis. Finally, we determined the CAF-associated genes by intersecting the GSE32894 and TCGA-BLCA datasets.

Functional enrichment analysis

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used for the functional enrichment analysis of the CAF-associated genes. The GO database describes the characteristics of genes and gene products of a species [28]. GO is divided into three non-overlapping domains: molecular function (MF), cellular component (CC), and biological process (BP). KEGG is a comprehensive database that helps to understand biological systems and high-level functions, and it processes and systematically describes biological processes [29].

Construction of a CAF-associated signature

The data of the TCGA-BLCA cohort were used to construct CAF-associated signatures as the training sets, and GSE32894, GSE13507 and IMvigor210 were used as the test sets. Univariate Cox regression analysis was used to determine the effect of prognostic CAF-associated genes on overall survival (OS). Subsequently, genes with p < 0.05 were included in the Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis, and 1000 iterative operations were performed [30]. The CAF score was calculated using the following equation: CAF score = Ʃ (β_i * Exp_i), where β_i denotes the LASSO coefficient of the ith gene and Exp_i denotes the expression value of the ith gene. Patients were divided into high- and low-CAF score groups according to the median CAF score, and Kaplan–Meier survival curves and log-rank tests were used to estimate the difference in OS between the two groups. In addition, the results of univariate and multivariate Cox regression analyses were combined with the clinical information of the patients to verify whether the CAF score was an independent prognostic risk factor.

Analysis of somatic alteration data and immunotherapy response prediction

We downloaded the somatic mutation data of the patients with BLCA and calculated the tumor mutation burden (TMB) using the tmb() function of the maftools package [31]. TMB values were defined as the total number of somatic mutations divided by the size of the sequenced region. This package was also used to visualize somatic mutations in patients with high or low CAF scores. Gene set enrichment analysis (GSEA) [32] determined KEGG functional pathways enriched for differentially expressed genes between high- and low-CAF score groups. The difference in response rate to immunotherapy between high- and low-CAF score groups was compared using the chi-square test, and the predictive power of CAF score was evaluated using the area under the curve (AUC) value of the receiver operating characteristic (ROC) curve.

Statistical analysis

All statistical analyses were performed using the R software (version 4.2.2; https://www.r-project.org/). The Wilcoxon test was used to compare the differences between the two groups of data. The chi-square test was used to analyze the difference in somatic mutation frequency between high- and low-CAF score groups. Spearman analysis was used to determine the correlation between TMB and CAF scores. A p-value of less than 0.05 indicated statistical significance.

Results

WGCNA revealed CAF-associated modules

There were no outliers in the TCGA-BLCA and GSE32894 cohorts. When the threshold was 4, R² was greater than 0.9, which ensured that the network was close to scale-free (Figure 1A and D). This value was also the minimum threshold to smooth the curve (Figure 1B and E); therefore, the β-value was 4. The dynamic cutting method was used to identify different gene modules, and the modules with higher similarity were merged, resulting in seven different gene modules (Figure 1C and F). The blue module had the highest correlation with CAF (Cor = 0.82; p = 4e-92) and stromal scores (Cor = 0.94; p = 1e-174) in the TCGA-BLCA cohort (Figure 1G and H), and scatter plots showed a strong correlation between MM and GS and CAF (Cor = 0.9; p = 1e-200). The brown module had the highest correlation with CAF (Cor = 0.81; p = 2e-53) and stromal scores (Cor = 0.94; p = 1e-108) in the GSE32894 cohort (Figure 1I and J), and scatter plots showed a strong correlation between MM and GS and CAF (Cor = 0.87; p = 1e-200). Therefore, these two modules were selected for subsequent analysis.

Figure 1. WGCNA results. Relationship between R² and different soft thresholds (a). TCGA-BLCA (D). GSE32894. Mean value of adjacency coefficients of genes corresponding to soft thresholds (B). TCGA-BLCA (E). GSE32894. The genes with similar expression patterns were clustered into co-expression modules using a cluster dendrogram (C). TCGA-BLCA (F). GSE32894. Heat map of association between gene modules and CAF traits (G). TCGA-BLCA (I). GSE32894. Scatter plot of the correlation between MM and GS in blue (H) and brown (J) modules.

Functional analysis of the hub genes

We found 1726 genes in the blue module, 915 genes in the brown module, and 723 genes at the intersection of the two modules (Figure 2A). GO analysis revealed that these genes were mostly enriched in positive regulation of cell adhesion for BP, collagen-containing extracellular matrix for CC, and extracellular matrix structural constituent for MF (Figure 2B). The KEGG analysis (Figure 2C) revealed that these genes were mostly enriched in focal adhesion, phagosome, and cell adhesion molecules.

Figure 2. (A). Venn plot showing the intersection of TCGA-STAD blue module and GSE32894 brown module genes. GO (B) and KEGG (C) results of the hub genes. (D-E). LASSO regression analysis with adjustment parameters (λ) calculated according to partial likelihood bias and cross-validation.

Construction of the CAF-associated signature

TCGA had the largest sample size in these three datasets and was used as the training set. Univariate Cox analysis identified 237 genes from the hub genes that affected OS. LASSO analysis (Figure 2D and E) was then performed on these 237 genes to identify CAF signatures containing 27 genes (Supplementary Table S1). The median CAF score was used as the cutoff value to classify patients with BLCA in the three cohorts into high-CAF and low-CAF groups. Kaplan–Meier survival curves revealed that patients with high CAF scores had a worse prognosis than those with low CAF scores in all three cohorts (Figure 3A, D, and G). In addition, the combination of univariate (Figure 3B, E, and H) and multivariate (Figure 3C, F, and I) Cox regression analyses and the clinical information of the three cohorts revealed that the CAF score was an independent prognostic factor after adjusting for some confounding factors. In the TCGA cohort, patients with high CAF scores had significantly shorter progression-free survival than patients with low CAF scores (Supplementary Figure S1A), and CAF scores were independent risk prognostic factors (Supplementary Figure S1B and C). After validation using the GSE13507 cohort, it was found that patients with high CAF scores had lower OS (Supplementary Figure S1D and E), but CAF scores were not independent prognostic risk factors (Supplementary Figure S1F). The TCGA database had the most complete clinical information among the three cohorts; therefore, we combined this information with the prognostic factors identified by multivariate Cox regression analysis to construct a nomogram (Figure 3J). The time-dependent ROC curve verification showed that the AUCs of 1-, 3-, and 5-year OS were above 0.8 (Figure 3K), and the calibration curve (Figure 3L) also supported the good predictive ability of the nomogram.

Figure 3. Kaplan–meier analyses of patients having high and low CAF scores (a). TCGA-BLCA, (D). GSE32894, (G) IMvigor210. Univariate Cox regression analysis (B). TCGA-BLCA, (E). GSE32894, (H) IMvigor210. Multivariate Cox regression analysis (C). TCGA-BLCA, (F). GSE32894, (I) IMvigor210. (J). nomogram of the TCGA-BLCA cohort. (K). time-dependent ROC curve of the nomogram. (L). calibration curve.

GSEA and correlation of CAF scores and somatic variants

ECM receptor interaction, focal adhesion, and dilated cardiomyopathy were the main KEGG pathways in the high-CAF score group of the TCGA-BLCA cohort (Figure 4A). Allograft rejection, metabolism of xenobiotics by cytochrome, and pentose and glucuronate interconversions were the main KEGG pathways in the low-CAF score group (Figure 4B). In addition, we compared the expression levels of some genes associated with the immunotherapy response [33,34] between high- and low-CAF score groups. HAVCR2 and SNCA were highly expressed in the high-CAF score group, whereas TBX2 was highly expressed in the low-CAF score group (Figure 4C). Figure 4D and E shows the distribution of the top 20 genes with the highest mutation frequency in different CAF score groups. The most frequently mutated genes included TP53, TTN, KMT2D, and MUC16, and the mutation frequency of TP53 in the high-CAF group was significantly higher than that in the low-CAF group (Table 1). Spearman analysis showed no significant association between TMB and CAF scores (Figure 4F); however, a statistically significant difference was observed in TMB of patients with high and low CAF scores (Figure 4G).

Figure 4. GSEA of KEGG Pathway in high (A) and low (B) CAF score groups (C). Comparison of some immunotherapy-related genes between patients with high and low CAF scores. Waterfall plot of somatic mutations in patients with high (D) and low (E) CAF scores. (F). correlation analysis of TMB and CAF scores. (G). Comparison of TMB between high- and low-CAF score groups.

Table 1. Association of CAF scores with somatic variants.

Gene	H-wild	H-mutation	L-wild	L-mutation	P value
FBN2	195(97.5%)	5(2.5%)	177(89.39%)	21(10.61%)	0.002
NAV3	195(97.5%)	5(2.5%)	178(89.9%)	20(10.1%)	0.004
ADGRV1	188(94%)	12(6%)	168(84.85%)	30(15.15%)	0.005
NFE2L2	181(90.5%)	19(9.5%)	193(97.47%)	5(2.53%)	0.007
LRP1B	188(94%)	12(6%)	169(85.35%)	29(14.65%)	0.008
DYNC1H1	193(96.5%)	7(3.5%)	177(89.39%)	21(10.61%)	0.010
HUWE1	195(97.5%)	5(2.5%)	181(91.41%)	17(8.59%)	0.015
USH2A	191(95.5%)	9(4.5%)	175(88.38%)	23(11.62%)	0.015
AHNAK	188(94%)	12(6%)	171(86.36%)	27(13.64%)	0.017
SZT2	196(98%)	4(2%)	183(92.42%)	15(7.58%)	0.018
TP53	92(46%)	108(54%)	115(58.08%)	83(41.92%)	0.021
CCDC88A	193(96.5%)	7(3.5%)	179(90.4%)	19(9.6%)	0.024
MED13	190(95%)	10(5%)	175(88.38%)	23(11.62%)	0.027
SPTAN1	187(93.5%)	13(6.5%)	171(86.36%)	27(13.64%)	0.028
PRKDC	193(96.5%)	7(3.5%)	180(90.91%)	18(9.09%)	0.036
CHD4	195(97.5%)	5(2.5%)	183(92.42%)	15(7.58%)	0.037
RELN	190(95%)	10(5%)	176(88.89%)	22(11.11%)	0.040
DNAH3	190(95%)	10(5%)	176(88.89%)	22(11.11%)	0.040
DIDO1	190(95%)	10(5%)	176(88.89%)	22(11.11%)	0.040
ELF3	183(91.5%)	17(8.5%)	167(84.34%)	31(15.66%)	0.042
AHR	194(97%)	6(3%)	182(91.92%)	16(8.08%)	0.046

Association of CAF scores and response to immunotherapy

A TIDE score greater than 0 indicated a non-response to immunotherapy, whereas a score less than 0 indicated a response to immunotherapy in both cohorts. The IMvigor210 dataset contains different categories of real-world patient responses to immunotherapy. Complete response and partial response represent a response to immunotherapy, whereas stable disease and progressive disease represent no response to immunotherapy. In the TCGA-BLCA, GSE32894, and IMvigor210 cohorts, patients with immunotherapy response (complete/partial) had significantly lower CAF scores than those without immunotherapy response (Figure 5A, C, and F). Moreover, the proportion of patients with immunotherapy responses in the low-CAF group was significantly higher than that in the high-CAF group (Figure 5B, D, and G). The AUC for CAF scores to predict immunotherapy response was 0.715 for the TCGA-BLCA cohort, 0.715 for the GSE32894 cohort, and 0.586 for the IMvigor210 cohort (Figure 5C, E, and H).

Figure 5. Comparison of CAF scores between immunotherapy responders and nonresponders (A). TCGA-BLCA, (C). GSE32894, (F) IMvigor210. Distribution of responders and nonresponders in the high- and low-CAF score groups (B). TCGA-BLCA, (D). GSE32894, (G) IMvigor210. AUC of ROC curve (C). TCGA-BLCA, (E). GSE32894, (H) IMvigor210.

Discussion

At present, there is still controversy about the occurrence and development of BLCA, but some studies have identified its potential mechanisms [35,36]. The occurrence and development of tumors are not only modulated by their genes and epigenetic defects but also by the TME [37,38]. Fibroblasts in the TME play a key role in tumor development [39]. The population of CAFs is complex and heterogeneous, which may be attributed to their multiple origins [40]. The TME includes cellular and non-cellular components (ECM) [41]. The structure of ECM affects important cell functions, such as cell proliferation, cell differentiation, and tumor formation [42,43]. CAFs play an important role in promoting cancer; therefore, they may be important targets for cancer treatment. For example, losartan reduced the collagen and hyaluronic acid produced by CAFs in breast and pancreatic cancers, thereby affecting vascular perfusion and drug delivery in these malignant tumors [44].

Here, we used EPIC and ESTIMATE algorithms to quantify the infiltration of CAFs in the BLCA samples. Then, CAF-associated modules and hub genes were identified through WGCNA of the data of TCGA-BLCA and GSE32894 cohorts. The biological functions and signaling pathways of these genes were mostly related to ECM. Next, we constructed a signature of 27 genes and calculated the CAF score through univariate Cox and LASSO regression analyses. Finally, we found that the prognosis of patients with high CAF scores in the three datasets was significantly poor than that of patients with low CAF scores, and CAF scores were independent prognostic risk factors. GSEA suggested that genes in high-risk groups were enriched in the ECM receiver interaction. CAF score had no significant correlation with TMB; however, it can be used to judge the response of patients to immunotherapy, thereby suggesting that CAF may be a new biomarker for predicting prognosis and response to immunotherapy in patients.

Some gene signatures constructed in this study are related to the occurrence, development, or prognosis of BLCA. For example, the function of CD8 + T cells is inhibited in patients with bladder cancer with high expression of HLA-E, which, in turn, leads to poor immunotherapy response [45]. Liu et al. found that the loss of EMP1 promotes the migration, resistance to ferroptosis, and oxidative stress of bladder cancer cells, which, in turn, favors their metastasis [46]. The overexpression of CERCAM is related to the poor prognosis of patients with BLCA. The overexpression of CERCAM markedly enhances the viability and invasive ability of bladder cancer cells in vitro, whereas silencing of CERCAM inhibits the growth of subcutaneously implanted tumors in vivo [47]. The high expression of TM4SF1 correlates with the high TNM stage and poor prognosis of patients with muscle-invasive bladder cancer, and knockdown of TM4SF1 can significantly inhibit the proliferation of cancer cells [48].

Some authors have reported the use of gene models to predict the prognosis of BLCA. For example, Deng et al. constructed an IFN-γ-related gene signature based on the TCGA database and combined it with age and stage to construct a nomogram. The authors found that the AUC values were 0.72, 0.71, and 0.74 for predicting 1-, 3-, and 5-year survival rates, respectively [49]. Wu et al. constructed a pyroptosis-related gene signature, which predicted 1-, 3-, and 5-year survival rates with AUC of 0.70, 0.72, and 0.71, respectively, in the TCGA database and 0.74, 0.67, and 0.65, respectively, in the GSE13507 validation [50]. Liu et al. constructed a gene signature related to cuproptosis, and the AUC value in the TCGA database was 0.708–0.757. The authors used four data sets for testing, but the AUC values of these test sets were between 0.500 and 0.768 [51]. Yang et al. constructed a ferroptosis-related gene signature, which predicted 1-, 3-, and 5-year survival rates with AUC of 0.694, 0.723, and 0.757 for 1-, 3-, and 5-year survival rates, respectively, in the TCGA database and 0.644, 0.637, and 0.628 for 1-, 3-, and 5-year survival rates, respectively, in the GSE13507 cohort [52]. The AUC value of most gene signatures was 0.7; the AUC of the CAF score in three cohorts reached 0.7 and that of the nomogram reached 0.8, indicating that our signatures had strong prediction ability for prognosis.

PD1 and PD-L1 have been widely used in clinical applications to predict the response of patients to immunotherapy; however, variable results have been obtained in different ICI studies [53–55]. A significant correlation was observed between the PD-L1 expression and drug efficacy in the phase 2 study on duvacumab, but the results still need to be verified through phase 3 controlled trials involving a larger sample size [56]. In addition to the above two biomarkers, TMB and related microsatellite instability and DNA mismatch repair deficiency may be new biomarkers for BLCA immunotherapy [57]. Presently, some gene signatures are also used to predict the effect of immunotherapy. For example, some authors have used the IMvigor210 cohort to test the predictive power of signatures. Zhang et al. constructed a signature associated with immune infiltration with a response rate of 26.7% in patients with high scores and 12.8% in those with low scores [58]. Feng et al. constructed a signature related to ferroptosis with a response rate of 37% in patients with low scores and 18% in those with high scores (AUC was 0.568) [59]. Zeng et al. constructed a 5-methylcytosine-related signature with a response rate of 24% in patients with low scores and 19% in those with high scores [60]. However, the predictive power of biomarkers in these studies was low, and a combination of multiple factors may provide a better prediction.

TP53,TTN and KMT2D are the three most common mutations in BLCA patients in this study and the frequency of mutations at these three loci has also been found to be high in other studies involving multiple cohorts [61,62]. Some scholars have found that the use of TP53/PIK3CA/ATM mutation classifier can accurately predict the benefit of immunotherapy in patients with bladder cancer, and have used the data of two cohorts to prove it [63]. Some mutations with significant differences among different CAF score groups were also found in this study. NFE2L2 was found to be a ferroptosis-related gene and associated with the prognosis of BLCA patients [64]. LRP1B is associated with cancer stem cell characteristics of BLCA and is also a poor prognostic biomarker [65]. AHNAK was found to be associated with cisplatin resistance in BLCA patients [66]. Timely and accurate treatment of BLCA patients is very important for the prognosis [67,68]. The biomarkers found in this study may bring new ideas for the treatment of BLCA.

Our study has some limitations. First, we used data from three cohorts for external validation. Our study involves clinical regression analyses based on public databases, and the results have not been validated in prospective clinical trials. In the future, we aim to collect relevant data and follow-up information, improve the external validation of the data, and obtain a more efficient and accurate prediction model. Second, the biological functions of the CAF-associated hub genes in our constructed gene signature could not be verified in further experiments. Therefore, we will explore the molecular mechanism and regulatory role of these hub genes in the occurrence and development of BLCA.

Conclusion

We identified the CAF-associated hub genes using WGCNA, constructed a relevant gene signature, and calculated the CAF score of each patient with BLCA. The signature accurately predicted the prognosis and immunotherapy response of patients with BLCA, which may provide new insights into the CAF-targeting therapeutic strategies for BLCA.

Medical ethics statement

Ethical approval was waived by the Institutional Ethics Committee because data were obtained from public databases, and all patients were de-identified to maintain confidentiality.

Acknowledgments

We thank Bullet Edits Limited for the linguistic editing and proofreading of the manuscript.

Data availability statement

The datasets generated and analyzed in this study can be found in the Gene Expression Omnibus repository (accession number: GSE32894 and GSE13507), The Cancer Genome Atlas-bladder urothelial carcinoma (TCGA-BLCA), and IMvigor210 databases.

Disclosure statement

The authors declare that they have no competing interests.

Additional information

Funding

The author(s) reported there is no funding associated with the work featured in this article.