Risk signature based on cancer-associated fibroblasts in bladder cancer using single-cell and bulk RNA data

In this study, scRNA-seq files from two muscle-invasive bladder cancer specimens were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo↗) from accession number GSE130001↗ (7). Clinical information (Table S1), transcriptomic data, and somatic mutation data related to BLCA were obtained from The Cancer Genome Atlas (TCGA) database (www.Cancer.gov/↗), microarray expression data and corresponding clinical data were obtained from the GPL6102↗ platform for GSE13507↗ and GSE32894↗ from the GPL6102↗ platform GPL6947↗ in the GEO database (8). TCGA, GSE13507↗ and GSE32894↗ were used as the training set and external validation set respectively, after removing batch effects via the “SVA” package.

The infiltration of CAF in the TCGA and GSE13507↗ groups was calculated using bulk RNA-seq-based EPIC (9). Similarly, xCell (10), MCP-counter (11), and ESTIMATE (12) were also used for both groups to predict the infiltration of CAF, as well as the StromalScore. “Surv_categorize” function, was used to calculate the optimal cutpoint value to distinguish between the high and low CAF and StromalScore groups in the TCGA and GSE13507↗ samples. The “survival” package served to analyze and compare the survival rates of low and high CAF (or StromalScore) by the Kaplan-Meier method to determine whether CAF and/or StromalScore levels correlate with survival in bladder cancer.

The “WGCNA” package was utilized to construct a weighted gene co-expression network analysis (WGCNA) in the TCGA cohort (13). After filtering the bulk RNA-seq data from TCGA and removing outliers, we constructed Pearson correlation matrices and generated weighted neighbor-joining matrices to emphasize strong correlations and penalize weak correlations. A scale-free network was constructed using the function powerEstimate soft threshold to select the best soft threshold power β = 4. A topological overlap matrix (TOM) was generated. The TOM-based correlation dissimilarity measure was set at a minimum number of genes/modules of 30, resulting in the generation of 11 modules. Next, we performed correlation analysis between modules and EPIC-based CAF infiltration with StromalScore, and modules with the highest correlation coefficients were regarded as candidates for correlation with differentially infiltrated immune cells. After selecting candidate modules, we defined |MM| (|module membership|) > 0.6 and |GS| ( |gene significance|) > 0.6 as the screening criteria for screening key genes in candidate modules.

The scRNA-seq data was then processed with the “Seurat” package. After combining two muscle-invasive bladder cancer specimens using the “merge” function, cells with excluded genes discovered in less than three cells and detected genes detected in fewer than 200 genes were eliminated, with the fraction of mitochondria confined to less than 5%. The scRNA-seq data from high-quality cells were normalized and the “FindVariableFeatures” function looked for highly variable genes for downstream analysis. Principal component analysis (PCA) was then performed on the highly variable genes to identify significant principal components (PCs). After the initial 12 PCs were selected. Cell clustering was then visualized using the t-SNE algorithm using the t-distribution random neighborhood embedding, and the “plotly” package was plotted in 3D. The “FindAllMarkers” function was used to detect marker genes for each cell cluster with a log2-fold change (FC) filter value > 1 and p_val_adj< 0.05. We made statements for each cell cluster based on the following marker genes: Epithelial cells (KRT19, CDH1, EPCAM); Fibroblasts (MMP2, EMILIN1, SFRP2); Myofibroblasts (ACTA2, PDGRB); Endothelial cells (KDR, VCAM1, AQP1, SEMA3G, CLDN5, PLVAP). To improve accuracy, we compared key genes of the WGCNA module, which is highly associated with CAF, with key genes with CAF characteristics.

To determine the function of key genes and reveal their underlying biological functions and potential mechanisms, we performed Gene Ontology (GO) and Kyoto Gene Encyclopedia and Genomic Pathway Enrichment Analysis (KEGG) using the clusterProfiler R package, p< 0.05 and FDR< 0.05.

To obtain CAF-related genes for which prognostic markers could be constructed, univariate Cox regression was performed to examine the correlation between CAF genes and overall survival (OS) in the TCGA dataset. CAF genes with a p-value< 0.05 were considered significant prognostic genes. To minimize the risk of overfitting, we then applied a least absolute shrinkage and selection operator (LASSO) Cox regression model via the “glmnet” R package. Based on the LASSO regression coefficients and gene expression, a risk score was calculated for each bladder cancer in the training set under the following equation:

This was subsequently categorized into high and low CAF risk groups based on the median TCGA CAF-riskScore and generalized to the validation group. Heat maps were generated to visualize the association between CAF-riskScore and candidate genes. The prognostic performance of CAF-riskScore was assessed by time-dependent subject operating characteristic (ROC) curves and survival analysis. Univariate and multivariate Cox regression analyses were performed to determine whether the CAF-riskScore independently served as a significant prognostic indicator. Nomograms were established for the GSE13507↗ and TCGA cohorts. Calibration curves were used to evaluate the predictive performance of these nomograms.

The proportion of immune cell infiltration in TCGA samples was calculated using CIBERSORT (14). We examined the differences in the infiltration of immune cells under risk grouping and the correlation of CAF-riskScore with certain cells in TME.

Spearman correlation estimated the relationship between the CAF riskScore and the CAF infiltration predicted by the above algorithm. Heatmaps were plotted by the “GGally” package. Spearman correlations similarly measured the CAF-riskScore and the relationship between candidate genes and known CAF-associated genes, and the “ggplot2” package did the plotting.

GSEA analysis was performed on the high CAF group from the TCGA cohort based on the Molecular Signature Database (MSigDB) (c2.cp.KEGG gene sets, hallmark gene sets), using default settings, and the top 5 and P< 0.05 pathways were plotted for each human collection gene sets.

Mutated genes were calculated using the Tumor Mutation Burden (TMB). The “maftools” R package was used to visualize the top 20 genes with the highest mutation frequency in both groups.

TIDE is a technique for assessing the possibility for immune evasion based on transcriptional profiling (15). TIDE allows the assessment of the effectiveness of immune checkpoint blockade for ICI therapy including anti-PD1 and anti-CTLA4 therapies. After using all tumor sample means as normalized controls, we downloaded TIDE estimates for each sample high and low CAF-risk groups of TCGA and GSE13507↗ cohorts and from the TIDE website (http://tide.dfci.harvard.edu/↗).“Ggpubr” and “reshape2” package were used for visualization.

Based on expression and sensitivity data from the GDSC2 database, the “oncoPredict” package was used to predict the drug sensitivity of different drugs in the high and low CAF-risk groups (16).

All cell lines were purchased from American-Type Culture Collection (ATCC). Human immortalized uroepithelial (SV-HUC-1) SV-HUC-1 cell line was cultured with Ham’s F-12K (HyClone, China)/10% fetal bovine serum (Gibco, Australia) media while BLCA cell lines (5637, T24) were cultured with RPMI 1640 (HyClone, China)/10% fetal bovine serum media. All cells were cultured in an incubator with 5% CO2 at 37°C.

TRIzol reagent (Thermo Fisher Scientific, USA) was used to extract RNA from the cells and the PrimeScriptTM RT Reagent Kit (TaKaRa, Japan) was used for reverse transcription into cDNA, respectively. After all target RNAs have been reverse transcribed to cDNA, relative quantitation by real‐time PCR was performed using TB Green PCT Master Mix (akara, Japan). And then qRT-PCR analysis was performed using a CFX96 real-time PCR system. GAPDH was used for experimental reference. All primers, the sequences of which are listed in, were synthesized by Sangon Biotech (Shanghai). 1

Statistical analysis was performed using R 4.2.0 and GraphPad Prism 8. The data package used for statistical analysis within R was as described above. Correlation matrices were analyzed using Pearson or Spearman for correlation. Survival analysis was performed using the Kaplan-Meier method and Log-rank tests. Comparisons between the two groups were made using the Wilcoxon test. A chi-square test was used for categorical variables. The differences in p-values< 0.05 were considered to be statistically significant.

Since CAF is of import in tumorigenesis, it is urgent to further elucidate the relationship between CAF and BLCA prognosis. The EPIC, xCell, and MCPcounter algorithms were used to calculate the amount of CAF cells in TCGA samples. BLCA patients were then divided into high and low fibroblast content groups. Kaplan-Meier analysis showed that the 3 algorithms had a longer survival rate in the low CAF content group than in the high CAF infiltration group (Figures 2A–C), thus suggesting that CAF play a major role in BLCA. Similar results were found for the ESTIMATE-based StromalScore (Figure 2D), with fibroblasts being an important component of the stromal cells.

Based on this observation, WGCNA was used to identify CAF-associated genes in bladder cancer. Supported by the depth of CAF infiltration in the EPIC package with StromalScore, WGCNA was utilized to identify CAF-associated genes in bladder cancer. First, samples above 120 were removed from the TCGA cohort (Figure 3A), selection 4 was chosen as the optimal soft threshold power (no scale R2 = 0.9254) (Figure 3B), and WGCNA identified 11 modules as shown in Figures 3C, D: of these, the brown module was significantly associated with high levels of CAF cells, (correlation = 0.9, p<0.0001), and the brown module The correlation between GS and MM is a key measure of the quality of gene module construction. After correlating the proportion of brown modules with CAF, the correlation between GS and MM reached 0.96 (Figure 3E), suggesting that all genes in the pink modules are specifically expressed by CAF in BLCA and that the expression levels of these genes are not easily influenced by other cells. Therefore, we set more stringent screening conditions of GS correlation > 0.6 and MM correlation > 0.6 and selected 981 key genes for downstream analysis (Table S3).

After pre-processing scRNA-seq data from GSE130001↗ based on the stringent quality control metrics described, 1623 high-quality cell samples were separated from the 2 found bladder cancer tissues. The number of genes detected (nFeature) and the depth of sequencing (total UMI, nCount), and the percentage of mitochondrial genes (percent.mt) were plotted (Figure 4A), with a strong positive correlation between nFeature and nCount and a Pearson correlation coefficient of 0.96 (Figure 4B). Subsequently, we used the t-SNE technique on the first 12 major components (Figures 4C, S1) to visualize the high-dimensional scRNA-seq data and successfully classified the cells into 12 clusters, followed by the following marker genes declared for each cell cluster (Figures 4D, E). In addition, the threshold for significantly expressed marker genes in each cluster was logFC > 1, adjPval< 0.05, and the top 10 markedly different genes for each cluster were shown by heatmap (Figure 4F). Of these, 220 genes were selected as CAF-related genes (Table S3).

After de-intersecting the key genes obtained from WGCNA and the CAF marker genes obtained from the single-cell analysis, 124 candidate CAF-related genes were obtained for subsequent analysis (Figure 5A; Table S3). As shown in Figure 6A, enriched GO terms were significantly enriched mainly in extracellular matrix organization, extracellular structure organization, external encapsulating structure organization, transmembrane receptor protein serine/threonine kinase signaling pathway and connective tissue development. Figure 6B shows the top 12 enriched KEGG pathways, which mainly include Focal adhesion, ECM-receptor interaction, Proteoglycans in cancer, TGF−beta signaling pathway and PI3K Akt signaling pathway. These enrichment terms enhance the reliability of the marker gene screen.

In the TCGA cohort, a total of 80 genes exhibited P< 0.05 by entering the 124 CAF marker genes mentioned above into univariate Cox regression analysis (Table S4). The final LASSO regression analysis identified seven prognostic marker genes: LDL Receptor Related Protein 1 (LRP1), Annexin A5 (ANXA5), Serpin Family E Member 2 (SERPINE2), Extracellular Matrix Protein 1 (ECM1), Retinol Binding Protein 1 (RBP1), Gap Junction Protein Alpha 1 (GJA1), FKBP Prolyl Isomerase 10 (FKBP10) (Figures 5B, C). CAF-riskScores were computed for each sample in TCGA based on the coefficients (Table S3) and the expression of prognostic genes. BCLA patients in the training set were then divided into low and high-risk groups based on the median riskScore, (Figure 5D), and extended to the validation group (Figure 5E). To assess the performance of the risk model, ROC curves were plotted and the area under the ROC curve (AUC) was 0.641, 0.654, and 0.666 for 1, 3, and 5 years in the training group, respectively (Figure 5F), and 0.576, 0.670 and 0.656 for 1, 3 and 5 years in the validation group (Figure 5G), respectively. In the GSE32894↗ database, CAF-riskScores also has good predictive performance. AUC for 1, 3 and 5 years was 0.687, 0.669 and 0.697 (Figure S2).

The above results suggest that the CAF-riskScore could be a promising prognostic biomarker for patients with bladder cancer. Therefore, to determine whether CAF-riskScore could be used independently as a prognostic indicator, we performed univariate and multivariate Cox regression analyses in the TCGA and GSE13507↗ cohorts. The results showed that the riskScore (HR=3.513, 95% CI:1.997-6.182, P<0.001), TNM Stage (HR=1.630, 95% CI:1.335-1.989, P<0.001), and age (HR=1.028, 95% CI:1.012-1.044, P<0.001) were independently associated with OS (Figure 7A), and in the GSE13507↗ cohort (Figure 8A), riskScore (HR=3.760, 95% CI:1.349-10.484, P=0.011) was consistently validated as an independent prognostic indicator. To predict OS rates for patients with bladder cancer in the TCGA (Figure 7B) and GSE62254↗ (Figure 8B) cohorts, we created a predictive nomogram for clinicians including the CAF-riskScore. Calibration curves showed that the model predictions for 1-, 3- and 5-year OS probabilities were consistent with the ideal predictions (grey line) in all datasets (Figures 7C, 8C). These results suggest that the nomogram model can be invoked as a reliable tool for predicting OS in bladder cancer patients. Two-by-two comparisons of OS in different risk groups were investigated by log-rank tests. Kaplan-Meier curves showed that the group with high CAF-risk tended to have significantly detrimental survival outcomes compared to the low CAF-risk group (TCGA cohort, hazard ratio (HR) = 2.134, 95% CI:1.565 -2.91, log-rank P<0.001, Figure 7D; GSE13507↗ cohort, HR=1.626, 95% CI:1.607-2.470, log-rank P=0.022, Figure 8D; GSE32894↗ cohort, HR=2.927, 95% CI:1.26-6.801, log-rank P=0.009, Figure S2D). Notably, the CAF-risk grouping for GSE13507↗ was provided by the median CAF-riskScore of the TCGA cohort.

Based on the Cibersort algorithm, we investigated the correlation between CAF-riskScore and TME components at the bulk RNA-seq level. As shown in Figure 9A, according to risk grouping, in terms of immune cells, we found that the high CAF group had significantly higher resting CD4 memory T cells, M0 Macrophages, M2 Macrophages, and Neutrophils, than the low CAF group. In contrast, Plasma cells, CD8 T cells, activated CD4 memory T cells, follicular helper T cells, resting NK cells, Monocytes, and activated Dendritic cells, were significantly higher in the high CAF group. Using correlation graphs to identify the correlation between CAF-riskScore and TME components. Notably, as the riskScore increased, the infiltration of CD8 T cells, CD4 memory T cells, and follicular helper T cells gradually decreased, while M0 Macrophages and M2 Macrophages gradually replaced some of the above cells (Figures 9B-G).

By running the EPIC, xCell, MCPcounter, ESTIMATE, and TIDE algorithms, we investigated the correlation between CAF-riskScore and the fraction of CAF in the bulk RNA-seq level. As shown in Figure 10A, the riskScore obtained from the CAF model was highly positively correlated with the degree of CAF infiltration obtained from EPIC, MCPcounter, ESTIMATE, and TIDE. We also discussed the relationship between CAF-related gene expression and the expression of the CAF-riskScore and its component genes. In the TCGA cohort, CAF-riskScore was positively correlated with all CAF-related genes except S1000A4, and p<0.05 (Figure 10B).

We then investigated the functional pathways in this model by GSEA analysis. Patients in the TCGA cohort were divided into high-CAF and low-CAF groups according to the above approach. We found immune and extracellular communication features in the engaged CAF group, including cell adhesion molecules, cytokine-cytokine receptor interaction, ECM-receptor interaction, and focal adhesion (Figure 10C). In addition, several hallmark gene sets were also significantly enriched in the high CAF-risk group, including allograft rejection, complement, epithelial mesenchymal transition, inflammatory response, and kras signaling up (Figure 10D).

The 20 most commonly mutated genes in the high (Figure 11A) and low (Figure 11B) CAF-risk groups are shown as waterfall plots, with some genes present in both groups, including TP53, TTN, KMT2D, ARID1A, MUC16, KMT2C, PIK3CA, SYNE1, RB1, HMCN1, KDM6A, RYR2, and FLG. MACF1, EP300, FAT4, XIRP2, CSMD3, KMT2A, and CUBN mutations were uniquely present in the high CAF-risk group, while OBSCN, ELF3, SPTAN1, BIRC6, NEB, and STAG2 mutations were uniquely seen in the low CAF-risk group.

Subsequently, using the TIDE online algorithm, we predicted the probability of response to immune checkpoint inhibitors in both datasets. In the TCGA cohort, as shown in Figure 10C, only 16% of patients in the high CAF-risk group had a response to immunotherapy, in the low CAF-risk group this was 60%, (Figure 11C, chi-square test, p<0.001). We then used the ROC curve to assess the relationship between the CAF-riskScore and the two groups of responders versus non-responders, with an AUC of 0.815, 95% CI:0.771-0.857. In the GSE13507↗ cohort (Figure 11D), again, also in the low-CAF risk group, more patients responded to immunotherapy (Figure 11E, 75% vs 28%), and AUC was 0.758 (Figure 11F). The drug sensitivity of each patient was then calculated by the oncoPredict algorithm to seek different chemotherapeutic agents for the CAF-related high and low CAF-risk groups. Based on the GDSC2 cancer cell line database, we found that higher CAF-riskScore increased TCGA in 6 anticancer drugs (acetalax, dihydrorotenone, gallibiscoquinazole, leflunomide, navitoclax, sinularin) sensitivity (Figures 12A-F).

Subsequently, the qRT-PCR measured the seven genes level (Figure 13). Thus, all the expression of CAF genes was higher based on the normal cell line in BLCA cell line.

These results suggest that while a high CAF-riskScore is associated with increased sensitivity to some chemotherapies, immunotherapy may be more effective in bladder cancer patients with low CAF-riskScore.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://figshare.com/articles/dataset/row-data_zip/22227037↗.