Cancer-associated fibroblast-derived gene signature discriminates distinct prognoses by integrated single-cell and bulk RNA-seq analyses in breast cancer

Single-cell transcriptome analysis

In our study, we initially obtained scRNA-seq data from 26 breast cancer samples from the GEO database (GSE176078) [13]. We performed unsupervised clustering of the individual cells using the read count matrix as input, employing the Seurat package (v4.1.1) in R (v4.1.3) [14]. Stringent quality control criteria were implemented, primarily focusing on the number of detected genes and the proportion of mitochondrial gene counts per cell.

Initially, we filtered out cells with fewer than 200 detected genes and cells with more than 15% mitochondrial gene counts. We used the Harmony algorithm to integrate the multi-sample data and correct batch effects [15]. Subsequently, dimension reduction clustering and differential expression analysis were performed following the Seurat-guided tutorial. Principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) dimension reduction were carried out using the top 20 principal components. Cell cluster annotations were based on canonical gene markers.

Weighted gene co-expression network analysis

We used the WGCNA package to obtain genes most related to CAFs content [16]. Samples were clustered to ascertain the overall relevance of all samples in the dataset, and outliers were excluded. The soft thresholding power β was chosen based on the lowest power for which the scale-free topology fit index reached a high value. The minimum gene number/module was set to 50 and, finally, 11 modules were generated. Next, we undertook correlation analyses between modules and traits to find the most relevant modules for CAFs content.

Collection of public datasets

RNA-sequencing expression matrix and clinical information of breast cancer samples and para cancerous tissues from The Cancer Genome Atlas (TCGA) database were downloaded from UCSC Xena (https://xena.ucsc.edu/). Three additional independent datasets (GSE20685, GSE37751 and GSE58812) were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) [17–19]. We downloaded somatic mutation data from Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/). Somatic mutation data sorted in the form of Mutation Annotation Format (MAF) were analyzed and then used to calculate Tumor mutation burden (TMB) using the R package maftools [20].

Construction and validation of a CAF-related prognostic signature

First, 1205 genes related to CAFs were collected from single-cell transcriptome analysis and 487 yellow module genes from WGCNA. Taking the intersection of the two yielded 341 candidate CAF-related genes for breast cancer. To obtain CAF-related genes that could construct a prognostic signature, univariate Cox regression and least absolute shrinkage and selection operator (LASSO) regression analyses were carried out. We eventually obtained 8 genes, including CERCAM, EMP1, SDC1, PRKG1, XG, TNN, WLS, and PDLIM4, and constructed a CAF-related prognostic model based on these genes. To group the breast cancer patients, the risk score of each breast cancer patient in the training set was calculated according to the following formula:

$$\text{Risk score}={\displaystyle \sum \text{ni}}={\displaystyle \sum (\text{Coefi}\times \text{xi})}$$

The breast cancer patients were then categorized into the high-risk and low-risk groups according to the median of risk score. The predictive sensitivity of the risk score was painted via the R package survival ROC for estimation [21]. The model effectiveness was evaluated in the validation set using the same coefficient and cutoff values that were used in the training set.

Biological functional analysis between high/low-risk group patients

The DESeq2 R package was used to perform differentially expressed genes (DEGs) analysis. DEGs were determined with a cutoff of an adjust p-value of less than 0.05 and |Log2 fold change| greater than 1 [22]. The clusterProfiler R package was used to perform gene set enrichment analysis (GSEA) [23]. With the use of Fisher’s exact test, those with false discovery rate FDR-corrected p-values of less than 0.05 were regarded as marked indicators. Single sample gene set enrichment analysis was performed via the R package GSVA [24]. Gene signatures of recurrent cancer cell states were collected from the previous study. The ITH score was calculated using the DEPTH R package [25].

Tumor immune microenvironment in breast cancer patients

To study the infiltration of immune cells, we used TIMER2.0, an efficient algorithm for predicting immune cell infiltration of bulk tumor gene expression data (http://timer.cistrome.org/). In addition, we collected series of tumor immunomodulators from the literatures and calculated the correlation of risk score with them.

Predicting drug responses and immunotherapy sensitivity

We used the R package oncoPredict to assess the predictive ability of risk score chemotherapeutic agents by calculating patients IC50 for various common chemotherapeutic agents. The Wilcoxon rank test was then used to compare the difference in IC50 between the high/low-risk groups.

Univariate and multivariable Cox regression

We performed univariate Cox regression on breast cancer patients with gene expression and overall survival. Multivariate Cox regression was used to evaluate independent risk factors in the same cohort. Genes and factors with a false discovery rate (FDR) <0.05 were considered statistically associated with patient survival. The results of univariate and multivariate Cox regression were acquired and visualized by using the R package forestplot.

Establishment of the nomogram

This study used the Cox regression model along with the R package rms to build an OS prediction nomogram that set 1-, 2-, 3-, and 5-year OS as the endpoints.

Statistical analysis

All statistical analyses were performed using R version 4.1.3 (https://www.r-project.org/) and its adequate packages. Statistical significance was set at p ≤ 0.05.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Interrogating the cellular constitution of breast cancer at single-cell resolution

In order to meticulously investigate the cellular constitution of breast cancer at single-cell resolution and identify cell markers of CAFs, we re-analyzed the scRNA-seq data of tumors from 26 breast cancer samples. Firstly, we integrated these data and corrected the potential batch effects through Harmony algorithm. Following rigorous quality control and data filtering, data of 29733 genes within 85408 cells was obtained. Three major compartments in the TME of breast cancer, including the immune subset, the stromal subset, and the epithelial subset were identified (Figure 1A). UMAP visualization showed that scRNA-seq data from different samples were integrated and mixed uniformly (Figure 1B). Using canonical lineage markers, we annotated each cell subpopulation among three main cellular subsets as epithelial cells, T cells, B cells, plasma cells, myeloid cells, endothelial cells, pericytes, cycling cells, and CAFs (Figure 1C, Supplementary Table 1). For instance, the immune subset consisted of T cells which were identified by expressions of CD2, CD3D and CD3E, B cells with high MS4A1, CD79A, and CD79B expressions, plasma cells with high IGHG1, IGKC, and JCHAIN expressions, and myeloid cells which were identified by significant expressions of LYZ, C1QA and C1QB (Figure 1C). Additionally, we found several stromal subpopulations including endothelial cells, pericytes, and CAFs. We annotated the endothelial cells due to the unique expression of ACKR1, PLVAP and PECAM1, as well as pericytes with expressions of RGS5, ACTA2, TAGLN (Figure 1C, 1D). Moreover, with a special focus on the CAFs, we observed CAFs showed significant expressions of DCN, COL1A1 and LUM, which were believed to be uniquely expressed proteins on the surface of CAFs and played a role in the development and function of CAFs. As expected, the Gene Ontology (GO) enrichment analysis exhibited that CAFs gene markers were enriched in pathways including extracellular matrix organization, extracellular structure organization and collagen fibril organization (Figure 1E). Collectively, we interrogated the cellular constitution of breast cancer in detail at single-cell resolution and identified cellular marker genes of each subpopulation, especially for CAFs.

Figure 1. Interrogating the cellular constitution of breast cancer at single-cell resolution. (A, B) UMAP plot showing the major cell subpopulations in breast cancer. (C) Bubble heatmap showing expression levels of selected signature genes in breast cancer. Dot size indicates fraction of expressing cells, colored based on normalized expression levels. (D) Feature plots to further identify various CAFs, based on the expression levels of marker genes. (E) GO enrichment of CAFs signature genes.

Screening for CAF-related genes by WGCNA in breast cancer

To comprehensively anatomize CAF-related biomarkers, we implemented WGCNA in the bulk RNA-seq dataset of breast cancer. Initially, following the removal of outliers from the TCGA samples, five were selected as the optimal soft-threshold power and 11 modules were identified by WGCNA algorithm (Figure 2A and Supplementary Figure 1A–1C). We compared data from three immune infiltration algorithms, including MCPCOUNTER, EPIC, and XCELL, as we expected, there was excellent congruity on the percentage of CAFs among these data (Figure 2B and Supplementary Figure 1D, 1E). In order to pinpoint the key modules of WGCNA associated with CAFs, we separately calculated the correlation of the three types of immune infiltration data with the modules. As shown in Figure 2C, the yellow module exhibited the highest correlation to the deconvolution result according to the MCPCOUNTER and EPIC, and equal importance in the XCELL result (Supplementary Table 2). Collectively, we designated the gene in the yellow module as a potential biomarker for CAFs in breast cancer (Figure 2D and Supplementary Figure 1F, 1G). In addition, astonishingly, the results of the GO enrichment analysis of the yellow module genes were highly consistent with the enrichment analysis of the CAFs marker genes identified by scRNA-seq data (Figure 2E). In conclusion, our analysis revealed a significant correlation between the yellow module and CAFs in breast cancer through the application of WGCNA.

Figure 2. Screening for CAF-related genes by WGCNA in breast cancer. (A) The cluster dendrogram constructing the gene modules and module merging. (B) Correlation plot of infiltration of CAFs by EPIC and MCPCOUNTER. (C) Correlation analysis of modules with traits yielded 11 modules, with the yellow module considered to be the most relevant module for CAFs. (D) Scatter plot between the yellow module and MCPCOUNTER. (E) GO enrichment of yellow module genes.

Construction of a CAF-related prognostic signature in breast cancer

By intersecting the 1205 CAFs marker genes with the 487 yellow module genes for fetching, 341 candidate CAF-related genes were obtained to be included in the downstream analysis (Figure 3A). Initially, univariate Cox regression analysis revealed that 12 genes were associated with the prognosis of breast cancer, followed by LASSO analysis to derive a prognostic signature composed of 8 genes (Figure 3B, 3C), including CERCAM, EMP1, SDC1, PRKG1, XG, TNN, WLS, and PDLIM4 (Supplementary Table 3). The TCGA breast cancer cohort was used as the training dataset, and risk score for each sample were computed using the coefficients and expression levels of the prognostic signature genes. We then categorized breast cancer patients into high/low-risk groups in the TCGA training cohort based on median risk score and discovered that patients in the low-risk group had markedly superior outcomes (Figure 3D–3F). Following this, the area under the receiver operating characteristic (ROC) curve (AUC) values for 1, 2, 3, and 5 years were 0.76, 0.69, 0.67, and 0.69 respectively (Figure 3G), which fully demonstrated the excellent performance of our risk model. In order to confirm the reliability and stability of the CAF-related prognostic signature, we validated it in several other additional independent validation datasets. Similarly, consistent with the training set, the breast cancer patients in the validation dataset were stratified based on the median risk score, revealing significantly worse outcomes for patients in the high-risk group (Figure 3H–3J). The AUC values for risk score in the GSE20685 dataset were 0.70 for 1 year, 0.73 for 2 years, 0.76 for 3 years and 0.74 for 5 years respectively (Figure 3K). In addition, breast cancer patients with higher risk score also had significantly shorter survival in the GSE37751 and GSE58812 datasets (Supplementary Figure 2). Overall, we have established and validated an innovative and robust CAF-related prognostic signature for predicting breast cancer prognosis.

Figure 3. Screening of CAF-related genes and construction a CAF-related prognostic signature in breast cancer. (A) The Venn graph of the CAF signature genes and yellow module genes. (B) Coefficient profiles in the LASSO regression model. (C) Cross-validation for tuning parameter selection in the LASSO regression. (D) Kaplan-Meier survival analysis was performed on the relationship between the risk score and OS using the TCGA training cohort. (E) The rank of risk score in the TCGA training cohort. (F) Survival status in the TCGA training cohort. (G) Time-dependent ROC curve analysis of the prognostic model (1, 2, 3, and 5 years) in the TCGA training cohort. (H) Kaplan-Meier survival analysis was performed on the relationship between the risk score and OS using the GSE20685 validation cohort. (I) The rank of risk score in the GSE20685 validation cohort. (J) Survival status in the GSE20685 validation cohort. (K) Time-dependent ROC curve analysis of the prognostic model (1, 2, 3, and 5 years) in the GSE20685 validation cohort.

Analyses of clinicopathological characteristics based on the CAF-associated prognostic signature

Besides the distinct survival outcomes observed between high/low-risk populations, notable variances were also observed in their clinicopathologic features. An incremental rise in the risk score was noted with advancing clinical stages, implying a potential correlation between the risk score and the progression of tumors (Figure 4A). Subsequently, we investigated the prognostic role of risk score in patients within tumor stage I-IV. Patients with lower risk score exhibited significantly better prognoses in stages II and III, while outcomes in stages I and IV did not show a similar improvement (Supplementary Figure 3A–3D). This phenomenon might be attributed to the limited number of patients with stage I and IV breast cancer and the growing emphasis on early breast cancer prevention. Among the different PAM50 molecular subtypes, the risk score exhibited a notable elevation in the LumB subtype and minimal in the normal-like subtype, which accounts for only a small fraction of breast cancers (Figure 4B). Similarly, we investigated the prognostic implications of risk score in patients within different PAM50 subtypes. Analysis revealed that patients with elevated risk score experienced poorer prognose in the Lum A/B and Basal-like subtypes, but not significantly in the Her2 and normal-like subtypes (Supplementary Figure 3E–3I). Furthermore, elderly patients had higher risk score, which may indicate that it played a role in the aging process (Figure 4C). Intra-tumor heterogeneity (ITH) is one of the crucial ingredients contributing to the failure of cancer therapies and patient mortality, higher ITH levels are associated with the development of therapeutic resistance [26]. Conspicuously, patients in high-risk group displayed elevated levels of ITH (Figure 4D). Risk score was positively correlated with ITH score, implying that high-risk patients were more resistant to single treatment (Figure 4E). Summarily, above results showed us high/low-risk breast cancer patients presented with unique clinicopathological profiles.

Figure 4. Functional and genomic features of CAF-related risk score-based classification. (A) Levels of risk score in different tumor stages of breast cancer. (B) Levels of risk score in different molecular subtypes of breast cancer. (C) Levels of risk score in different age groups of breast cancer. (D) Boxplot showing the levels of ITH in high/low-risk groups. Paired two-sided Wilcoxon test. (E) Scatter plot showing the correlation between the risk score and ITH score. (F) Bar plot showing different pathways enriched in high/low-risk groups of breast cancer calculated by GSEA. (G) Boxplots showing the signature score of 16 cancer cell states in high/low-risk groups of breast cancer scored by GSVA. Paired two-sided Wilcoxon test. The asterisks represent the statistical P-value (^*p < 0.05; ^**p < 0.01; ^***p < 0.001; ^****p < 0.0001; ^nsp > 0.05).

Functional analyses of the CAF-related prognostic signature

In order to investigate the underlying mechanisms explaining differences in outcomes and clinical characteristics between risk subgroups, we dissected the differences between subgroups in several ways. Initially, Gene set enrichment analysis (GSEA) analysis (Figure 4F) highlighted that breast cancer patients in the high-risk group showed a considerably enrichment in E2F targets, epithelial-mesenchymal transition (EMT) and angiogenesis (Supplementary Figure 4A–4C). Conversely, allograft rejection, pancreas beta cells and KRAS signaling DN were observed in the low-risk group patients (Supplementary Figure 4D–4F). Barkley D and colleagues innovatively developed multiple gene sets containing 16 recurrent cancer cell states which interacted with the TME to form organizational systems capable of promoting tumor progression, metastasis, and treatment failure [27]. To comprehensively gain insight into transcriptional heterogeneity between high and low risk groups, we calculated 16 recurrent cancer cell states using single sample gene set enrichment analysis (ssGSEA). As depicted in Figure 4G, low-risk patients demonstrated elevated score in astrocyte (AC)-like, alveolar, basic, ciliated, and glandular. Nevertheless, the cycle and hypoxia modules were enriched in high-risk patients. These data may illuminate the potential rationales for the poorer prognosis of high-risk patients.

Mutational features of the CAF-related prognostic signature

One of the most fundamental characteristics of cancer is preternatural and uncontrolled cell growth caused by genome-driven mutations, which in turn affects the homeostatic development of a range of essential cellular functions [28]. We described and compared the mutational features of patients in the high/low-risk groups. Somatic mutations were observed in 88.15% of patients, with a higher frequency of TP53 mutation, but lower frequency of PIK3CA and CDH1 mutations in the high-risk group (Figure 5A). TMB refers to the aggregate number of substitutions, insertion or deletion mutations per megabase in the exon coding region of a gene in the genome of a tumor cell [29]. Additionally, as an indicator of tumor mutations, TMB is frequently used to predict the efficacy of immunotherapy [30]. Apparently, patients in the high-risk group exhibited higher TMB levels (Figure 5B), with TMB showing a positively correlated with risk score (Figure 5C). Surprisingly, breast cancer patients stratified by the median TMB level did not show significant differences in overall survival expectancies (Figure 5D). Therefore, we investigated how the combination of TMB and CAF-related risk score could jointly categorize breast cancer patients into groups with significantly different prognoses (Figure 5E). To sum up, frequency of somatic mutation was higher in the high-risk group, and integration of CAF-related risk score and TMB could further refine the prediction of prognosis in breast cancer patients.

Figure 5. Mutational feature of the CAF-related prognostic signature. (A) Waterfall plot represents the mutation distribution of the most frequently mutated genes in high/low-risk groups. (B) Boxplot showing the levels of TMB in high/low-risk groups. Paired two-sided Wilcoxon test. (C) Scatter plot showing the correlation between the risk score and TMB in TCGA cohort. (D) Kaplan-Meier survival analysis was performed on the relationship between TMB and OS. (E) Kaplan-Meier survival analysis was performed on the relationship between combination of TMB and the risk score and OS. The asterisks represent the statistical P-value (^*p < 0.05; ^**p < 0.01; ^***p < 0.001; ^****p < 0.0001).

Tumor immune microenvironment of the CAF-related prognostic signature

In this section, we aimed to investigate the disparities in the tumor immune microenvironment (TIME) between the two risk groups. First, we assessed the distribution of immunocyte percentage in each patient by the CIBERSORT algorithm (Figure 6A). Apparently, M2 macrophages were more prevalent in high-risk patients, and correlation analysis showed a significant positive correlation between the percentage of M2 macrophages and the risk score (Figure 6B). In addition, CD8⁺ T cells, M1 macrophages, and activated NK cells were more enriched in patients in the low-risk group, with their infiltrations were negatively correlated with risk score (Figure 6C–6E). We also calculated immune cell infiltration utilizing additional algorithms, including CIBERSORT. ABS and MCPCOUNTER, and rightfully so these results are consistent with CIBERSORT (Supplementary Figure 5). Subsequently, we collected 65 immunomodulators from various sources in the literature, and as depicted in Figure 6F we calculated the correlation between these immunomodulators and the CAF-associated risk score. In particular, risk score was significantly positively correlated with the immunosuppressive molecules CD276, TNFSF9, and HAVCR2, and significantly negatively correlated with the stimulatory immune checkpoint markers SLEP and TNFRSF14.

Figure 6. Dissection of tumor immune microenvironment features between high/low-risk group. (A) Boxplots showing the proportion of 22 immune cells in high/low-risk groups of breast cancer estimated by CIBERSORT. Paired two-sided Wilcoxon test. (B–E) Scatter plots showing the correlation between the risk score and the proportion of M2 macrophages, CD8⁺ T cells, M1 macrophages, and activated NK cells. (F) Bar plot of the correlation between immunomodulators and the risk score in TCGA cohort. The asterisks represent the statistical P-value (^*p < 0.05; ^**p < 0.01; ^***p < 0.001; ^****p < 0.0001; ^nsp > 0.05).

In addition, we capitalized on the previously mentioned scRNA-seq data which contains 24 paired breast cancer bulk RNA-seq data. According to the median risk score calculated by bulk RNA-seq data, these 24 samples were stratified into high/low-risk group (Figure 7A). Patients in High/low-risk group displayed contrasting cellular subpopulation compositions (Figure 7B). Notably, myeloid cells, known for typically exerted an immunosuppressive effect, were more abundant in high-risk samples (Figure 7B). This finding corresponded with our CIBERSORT classification of immune cell infiltration analysis described above. To comprehensively gain an understanding of myeloid cells, we further re-clustered myeloid cells and distinguished macrophages, monocytes, dendritic cells (DCs) and cycling cells based on multiple cellular markers (Figure 7C, 7D). The results showed us that patients in the high-risk group had significantly higher composition of macrophages, which predominantly contribute to suppressing anti-cancer immunity (Figure 7E). Overall, these data illustrated distinct differences in the immune features of the two groups, with patients in the high-risk group demonstrating an immunosuppressive phenotype.

Figure 7. scRNA-seq analysis of the tumor immune microenvironment features between high/low-risk group. (A) UMAP plot showing the major cell subpopulations of high- and low-risk breast tumors. (B) Relative proportions of diverse cell types across high/low-risk tumors. (C) UMAP plot showing the diverse subsets of myeloid cells in breast cancers. (D) Bubble heatmap showing expression levels of selected signature genes for myeloid cells in breast cancers. Dot size indicates fraction of expressing cells, colored based on normalized expression levels. (E) Relative proportions of diverse subpopulations of myeloid cells across high/low-risk tumors.

Correlation between anti-cancer drug sensitivities and the CAF-related prognostic signature score

In order to delve deeper into the clinical implications of risk score and prognostic genes, we calculated the half-maximal inhibitory concentration (IC50) values for various drugs in breast cancer patients in TCGA using the GDSC database. Figure 8A showed the correlation between drug sensitivities and the risk score as well as prognostic genes. A notable positive correlation was observed between the risk score and IC50 of Oxaliplatin, Sabutoclax, and Irinotecan. However, the IC50 of BI.2536, a potent highly selective PLK1 inhibitor, negatively correlated with the CAF-associated risk score. Interestingly, the prognostic gene PRKG1 was highly correlated with the IC50 of Leflunomide (Supplementary Figure 6A), and TNN was negatively correlated with the IC50 of Nutlin.3a .... _1047, a potent MDM2 inhibitor (Supplementary Figure 6B). Additionally, we compared the IC50 levels of a series of Food and Drug Administration (FDA)-approved breast cancer treatments, and the IC50 levels were higher in the high-risk group for Palbociclib, Oxaliplatin, Camptothecin, Irinotecan, Mitoxantrone, and Sorafenib. These results suggested that breast cancer patients with lower risk score might be sensitive to these FDA-approved drugs (Figure 8B–8G). The I-SPY2 trial platform is a continuous, multicenter, Phase II neoadjuvant trial platform for high-risk, early-stage breast cancer [31]. We calculated the CAF-related risk score for each patient and found higher risk score for pathological complete response (pCR) patients in the Pertuzumab (Paclitaxel + Pertuzumab + Trastuzumab) and TMD1/P (T-DM1 + Pertuzumab) subgroups of the 10 cancer treatment arms for which they were advertised (Supplementary Figure 6C). In the GSE123845 dataset, which was breast cancer data after standard neoadjuvant therapy, patients with higher risk score demonstrated a favorable prognosis (Supplementary Figure 6D). In summary, the above findings defined that the CAF-related risk score could serve as a reputable tool for predicting drug sensitivity and response to cancer therapy for breast cancer patients.

Figure 8. High- and low-risk group patients differ in drug sensitivity. (A) Bubble plot showing the relationship between IC50 of drugs, risk score, and model genes. (B) Boxplot showing the comparison of IC50 of Palbociclib between high/low-risk groups, and scatter plot showing the correlation between the IC50 of drug and risk score. (C) Boxplot showing the comparison of IC50 of Oxaliplatin between high/low-risk groups, and scatter plot showing the correlation between the IC50 of drug and risk score. (D) Boxplot showing the comparison of IC50 of Camptothecin between high/low-risk groups, and scatter plot showing the correlation between the IC50 of drug and risk score. (E) Boxplot showing the comparison of IC50 of Irinotecan between high/low-risk groups, and scatter plot showing the correlation between the IC50 of drug and risk score. (F) Boxplot showing the comparison of IC50 of Mitoxantrone between high/low-risk groups, and scatter plot showing the correlation between the IC50 of drug and risk score. (G) Boxplot showing the comparison of IC50 of Sorafenib between high/low-risk groups, and scatter plot showing the correlation between the IC50 of drug and risk score. The asterisks represent the statistical P-value (^****p < 0.0001).

Constructing a nomogram of the CAF-related prognostic signature

Both univariate and multifactorial Cox analyses revealed that age, tumor stage, and CAF-related risk score were independent prognostic factors in breast cancer patients (Figure 9A, 9B). Therefore, a nomogram was developed incorporating the CAF-related risk score with additional clinical information on other independent predictors proposed through multivariate Cox analysis. In the nomogram, CAF-related risk score contributed significantly to the prediction of survival probability, serving as a quantitative and visual tool for predicting 1-, 3-, and 5-year OS (Figure 9C). The AUCs for 1-, 2-, 3-, and 5-year OS in the nomogram were 0.847, 0.832, 0.763, and 0.758, respectively, suggesting a significantly better prognostic ability than the CAF-related risk score alone (Figure 9D). In addition, calibration curves were presented to assess the performance of the nomograms, which showed that the predicted curves of the model are close to the ideal ones (Figure 9E). These results highlighted that significant predictive impact of the nomogram model on breast cancer patients.

Figure 9. Establishment and assessment of the nomogram survival model. (A) Univariate analysis for the clinicopathologic characteristics and the risk score in TCGA cohort. (B) Multivariate analysis for the clinicopathologic characteristics and the risk score in TCGA cohort. (C) A nomogram was established to predict the prognosis of breast cancer patients. (D) Time-dependent ROC curve analysis of the nomogram (1, 2, 3, and 5 years) in TCGA cohort. (E) Calibration plots showing the probability in TCGA cohort. The asterisks represent the statistical P-value (^*p < 0.05; ^**p < 0.01; ^***p < 0.001; ^****p < 0.0001).