The transcript per million (TPM) format RNA-seq data and clinical information relevant to Breast Invasive Carcinoma (TCGA-BRCA) samples were obtained from TCGA datasets (https://portal.gdc.cancer.gov/). The gene expression profiling dataset (GSE47994) was obtained from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/) database (30). We then screened a total of 694 LBC patients from the TCGA-BRCA dataset and 101 LBC patients from the GSE47994 dataset who had prognostic data available, with follow-up times exceeding 365 days.

Four methods were utilized to estimate the abundance of CAFs, including the Estimate the Proportion of Immune and Cancer cells (EPIC) algorithm based on cell-type deconvolution using constrained least square optimization (31), the xCell algorithm based on gene signature enrichment (32), the microenvironment cell populations-counter (MCP-counter) based on marker gene expressions (33), and the Tumor Immune Dysfunction and Exclusion (TIDE) algorithms, which also can predict anti-PD1 and anti-CTLA4 responses in tumor patients (34). The first three methods were executed via the deconvolute() function of the immunedeconv R package (version 2.0.3) (35), while the TIDE method was implemented through http://tide.dfci.harvard.edu/. Additionally, the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm was utilized to calculate the stromal score, which indicates the level of stromal infiltration in each tumor sample, through the estimate R package (version 1.0.13) (36).

The WGCNA R package (version 1.72) was utilized to construct co-expression networks and identify hub genes that target cancer-associated fibroblast (CAF) infiltrations and stromal scores (25). The input genes for network construction were selected based on the median absolute deviation (MAD), with the top 5,000 genes selected for both the TCGA-BRCA and GSE47994 cohorts. The Pearson’s correlation similarity matrix was calculated between each pair of genes (s_ij, where ij represents the gene pairs) and raised to a soft-thresholding power β ( S ij β ), based on the scale-free topology network criterion. The adjacency matrix was then clustered using the topological overlap measure (TOM) and dissimilarity (1-TOM) between genes, and a dynamic tree cut algorithm was applied to the dendrogram to identify gene modules with a minimum of 30 genes in each module. The first principal component of each module’s expression was summarized as a module eigengene (ME), and the Pearson’s correlations between MEs and EPIC-quantified CAF infiltrations, as well as the stromal score, were assessed to identify the most correlated module for further analysis. Hub genes were then identified by overlapping the most correlated module genes between the TCGA-BRCA and GSE47994 cohorts.

The clusterProfiler R package (version 3.14.3) was utilized to analyze the hub genes’ biological functions, which included biological processes (BPs), molecular functions (MFs), and cellular components (CCs), as well as pathways according to GO and KEGG databases (37). p < 0.05 was considered statistically significant.

For CAF risk model construction, 694 LBC cases from TCGA-BRCA were selected based on their large sample size. The validation cohort consisted of 101 LBC cases from GSE47994 cohort. Univariate Cox regression analysis was performed to identify prognostic stromal CAF hub genes for overall survival (OS). Genes with p < 0.05 were selected for LASSO Cox regression analysis with 1,000 iterations using glmnet R package to reduce the number of genes (38). Then, the CAF risk model was constructed as follows: CAF risk score = ∑ (β_i * Exp_i), where β_i is the LASSO coefficient of ith gene and Exp_i is the expression value of ith gene. Using the median CAF risk score of the training cohort as the threshold, LBC patients from both cohorts were classified into high- and low-CAF-risk groups and the OS difference between the two groups was estimated using Kaplan–Meier curves and the log-rank test.

We collected CAF specific and nonspecific markers from published literature (39, 40). To verify the reliability of our CAF model markers in LBC, we examined the Spearman’s correlations between the CAF risk score and the stromal score, as well as various CAF infiltration estimates (EPIC, xCell, MCP-counter, and TIDE). We also analyzed the correlations between CAF model genes and published CAF markers in both TCGA-BRCA and GSE47994 cohorts.

Using the Genomics of Drug Sensitivity in Cancer [GDSC (https://www.cancerrxgene.org/)] database (41), half-maximal inhibitory concentration (IC50) values of common drugs (bleomycin, lapatinib, paclitaxel, camptothecin, cisplatin, docetaxel, methotrexate, and sunitinib) in each LBC sample were estimated based on the transcriptome data by ridge regression with ten-fold cross-validation in pRRophetic R package (version 0.5) (42). The TIDE (http://tide.dfci.harvard.edu/) online algorithm was then adopted for immune checkpoint blockade therapy response predictions (34). The chi-squared test was used to examine differences in response rates between high- and low-CAF-risk groups. The predictive efficacy of the CAF risk signature was evaluated by ROC curves and area under the curve (AUC) values.

To explore the enriched hallmark and KEGG pathway gene sets between high- and low-CAF-risk groups in GSE47994, gene set enrichment analysis (GSEA) was performed using the enrichplot and clusterProfiler R packages. The gene sets used were derived from the Molecular Signatures Database (MSigDB), specifically the “c2.cp.kegg.v7.4.symbols” and “h.all.v7.4.symbols” gene sets (43). Additionally, the enrichment scores of ECM receptor interaction, regulation of actin cytoskeleton, and TGF-β signaling pathway hallmark gene sets were calculated using ssGSEA (44). Finally, Spearman’s correlation analysis was performed to assess the correlation between the CAF risk score and gene set enrichment scores.

To validate the findings at the cellular level, mRNA expressions of the identified markers in 38 fibroblasts and 51 BC cell lines were downloaded from the Cancer Cell Line Encyclopedia [CCLE (https://portals.broadinstitute.org/ccle)] database (45). Expression patterns of these markers in fibroblasts and CRC cell lines were examined using heat maps and Wilcoxon tests. Additionally, immunohistochemical (IHC) staining images of these markers in LBC tissues were downloaded from the Human Protein Atlas [HPA (https://www.proteinatlas.org/)] online database (46).

R software (version 4.2.2; https://www.r-project.org/) was used for all statistical analyses. The median CAF risk score was used as the cutoff value for each cohort to divide LBC patients into high- and low-CAF-risk subgroups. Pairwise comparisons were performed using the Wilcoxon test. Pearson correlation coefficient analysis was performed to evaluate the correlation between genes. Overall survival comparisons were made using the Kaplan-Meier curve with the log-rank test, which were adopted via the survival and survminer R packages. p < 0.05 was considered statistically significant.

The flowchart of this research is displayed in Figure 1 . CAF infiltrations were predicted by multiple methods, including EPIC, xCell, MCP-counter, and TIDE. The stromal score was calculated by the estimate algorithm. Their prognostic values on overall survival (OS) were evaluated via log-rank tests. Kaplan-Meier curves indicated that multiple higher CAF infiltrations and stromal scores were significantly associated with poorer OS in LBC patients. CAF_EPIC, CAF_TIDE, and stromal scores were significantly associated with poorer OS in TCGA-BRCA ( Figures 2A–C ), and CAF_EPIC, CAF_mcp-counter, and stromal scores were significantly associated with poorer OS in GSE47994 ( Figures 2D–F ), which highlights the importance of further studies exploring CAF and stromal-associated genes in LBC. In this study, the EPIC-estimated CAF abundances and stromal scores were summarized as phenotype data for subsequent analysis, and the data from the other three estimated CAF infiltrations were used for external validation of the identified CAF model.

WGCNA analysis was conducted on both TCGA-BRCA and GSE47994 datasets. To build a scale-free topology network, we estimated the soft threshold power (β) of 7 in TCGA-BRCA (scale-free R2 = 0.86) ( Figure 3A ) and 17 in GSE47994 (scale-free R2 = 0.86) ( Figure 3B ). In TCGA-BRCA, the hierarchical clustering tree identified 8 co-expression models ( Figure 3C ), and the magenta module exhibited the strongest positive correlation with the CAF proportion (Cor = 0.54, P = 7e-54) and stromal score (Cor = 0.78, P = 5e-144) ( Figure 3E ). In GSE47994, the dynamic hybrid cutting clustered 6 co-expression models ( Figure 3D ), with the brown module showing the strongest positive correlation with the CAF proportion (Cor = 0.88, P = 5e-25) and stromal score (Cor = 0.92, P = 2e-31) ( Figure 3F ). Therefore, we focused on these two modules for further investigations. A total of 1718 and 158 genes were included in the magenta and brown modules, respectively. In the magenta module, scatter plots indicated strong correlations between MM and GS for CAF (Cor = 0.64, p =1.3e-198) and stromal scores (Cor = 0.88, p < 1e-200) ( Figure 3G ). Similarly, in the brown module, strong correlations were observed between MM and GS for CAF (Cor = 0.64, P = 1.1e-19) and stromal scores (Cor = 0.82, p = 7e-40) ( Figure 3H ). Consequently, we selected 1718 genes from the TCGA-BRCA magenta module and 158 genes from the GSE47994 brown module as highly associated with CAF and stromal scores.

The above CAF-related genes were overlapped and screened to 135 hub genes ( Figure 4A ). These genes were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. The major enriched GO terms were related to extracellular matrix organization, extracellular structure organization, external encapsulating structure organization (BP), collagen-containing extracellular matrix and endoplasmic reticulum lumen (CC), and extracellular matrix structural constituents and collagen binding (MF) ( Figure 4B ). The main enriched KEGG pathways were human papillomavirus infection, the PI3K-Akt signaling pathway, and protein digestion and absorption ( Figure 4C ).

The 694 LBC samples from TCGA- BRCA were used as the training cohort owing to the larger sample size, and 101 GSE47994 samples were used as the validation group. By performing univariate Cox regression analysis of the 135 common hub genes, 20 OS-related genes with p < 0.05 were screened out and subjected to the following LASSO Cox regression analysis ( Figures 4D, E ). Five genes were finally identified for the CAF risk model construction: CAF risk score = expression of RIN2* 0.103 + expression of THBS1* 0.022 + expression of IL1R1* 0.068 + expression of RAB31* 0.055 + expression of COL11A1* 0.082 ( Figure 4F ). The median CAF risk score of the training cohort was 1.85, which was used as the cut-off value to classify LBC patients from each cohort into high- and low-CAF-risk groups. LBC patients in each cohort were divided into high– and low–CAF-risk groups with the median risk score as the cutoff value. Kaplan–Meier curves revealed that LBC patients in the high–CAF-risk group experienced worse OS than those in the low–CAF-risk group in both TCGA- BRCA (HR = 2.394, 95%CI: 1.531–3.745, log-rank p < 0.001) ( Figure 4G ) and GSE47994 (HR = 1.627, 95%CI: 1. 036−2.558, log-rank p = 0.032) ( Figure 4H ). These results indicated CAF and stromal-related signature genes were crucial prognostic markers in LBC.

To evaluate the robustness of the CAF model as an indicator of CAF infiltrations, Spearman’s correlation analyses were performed between the CAF risk score and stromal score as well as CAF abundances predicted by EPIC and three other methods: xCell, MCP-counter, and TIDE. The CAF risk score was strongly and positively correlated with multi-estimated CAF infiltrations and the stromal score in both TCGA-BRCA ( Figure 5A ) and GSE47994 ( Figure 5B ) cohorts. These results confirmed that the CAF risk score was a reliable predictor of CAF infiltrations. To further validate the correlation of the expression levels of the five genes with CAFs, their expression levels were compared with a set of collected CAF markers in both TCGA-BRCA ( Figures 5C, E ) and GSE47994 ( Figures 5D, F ) cohorts. A high and positive correlation was observed between the expression levels of the five genes and most of the CAF markers in both cohorts. These results demonstrated that the five genes were representative of CAFs.

The standard treatment for LBC patients involves radical surgery followed by adjuvant chemotherapy and endocrine therapy (47). IC50 values for multiple anti-tumor drugs, including those used in breast cancer treatment, were estimated using the GDSC database for both the TCGA-BRCA ( Figure 6A ) and GSE47994 ( Figure 6B ) cohorts. Wilcoxon analyses indicated significant differences in IC50 values between high- and low-CAF-risk LBC patients. As for commonly used chemotherapy drugs for breast cancer, although not both datasets showed statistical significance, the results indicated that high-CAF-risk patients were sensitive to Alpelisib and Epirubicin, while low-CAF-risk patients were sensitive to Docetaxel, Fulvestrant, Lapatinib, Palbociclib, Ribociclib and Tamoxifen. However, Cyclophosphamide and Zoledronic acid exhibited different trends between the two datasets. In addition, the results from both datasets showed that high-CAF-risk patients were insensitive to several other drugs, including Axitinib, Dabrafenib, Irinotecan, Sorafenib, Topotecan, and Venetoclax. This suggests that higher CAF risk scores are more likely to induce resistance to these drugs in breast cancer patients. Immunotherapy is among the most important advances in recent oncology, particularly for triple-negative and HER-2-positive breast cancer (48). Trials are also underway to assess the efficacy of immune checkpoint inhibitors such as pembrolizumab for Luminal breast cancer (49–51). To evaluate the CAF risk score as an immunotherapy predictor for LBC patients, the TIDE method was used. In TCGA-BRCA, the non-responder subgroup (n = 492) exhibited significantly higher CAF scores than the responder subgroup (n = 202) (p < 2.2e-16; Figure 6C ). Low-CAF-risk patients (144/347) displayed higher immunotherapy sensitivity and lower TIDE scores than high-CAF-risk patients (58/347) (p < 0.001; Figures 6D, E ). In GSE47994, the non-responder subgroup (n = 63) also had a significantly higher CAF score than the responder subgroup (n = 38) (p = 1.9e-7; Figure 6G ). Low-CAF-risk patients (31/54) exhibited higher immunotherapy sensitivity and lower TIDE scores than high-CAF-risk patients (7/47) (p < 0.001; Figures 6H, I ). The AUC values of 0.711 in TCGA-BRCA ( Figure 6F ) and 0.811 in GSE47994 ( Figure 6J ) indicate the excellent performance of the CAF model for predicting immunotherapy response.

To investigate the functional enrichment of the CAF signature, GSEA was conducted on the TCGA-BRCA dataset to compare the high- and low-CAF-risk groups. The analysis revealed significant enrichment in KEGG signaling pathways associated with ECM receptor interaction, focal adhesion, pathways in cancer, regulation of actin cytoskeleton, and tight junction in the high-CAF-risk group ( Figure 7A ). Moreover, the genes in the high-CAF-risk group were significantly enriched in Hallmarker gene sets related to angiogenesis, apical junction, coagulation, epithelial-mesenchymal transition (EMT), and inflammatory response ( Figure 7B ). Additionally, ssGSEA results indicated that the CAF risk score was positively correlated with enrichment scores for ECM receptor interaction, regulation of actin cytoskeleton, and TGF-β signaling pathway in both TCGA-BRCA ( Figure 7C ) and GSE47994 ( Figure 7D ).

Based on the CCLE database, the mRNA expressions of the five hub genes (RIN2, THBS1, IL1R1, RAB31, COL11A1) were verified to be higher in fibroblast cell lines than those in BC cell lines (Wilcoxon test, all p < 0.001; Figures 8A, B ). Additionally, to determine the protein expression characteristics of these CAF signature genes, the IHC images from the HPA database were analyzed. The data demonstrated that these proteins (RIN2, THBS1, IL1R1 and RAB31) were deeply stained in BC stroma ( Figure 8C ). These verifications suggest that these genes might be CAF-specific markers.

It is well known that the HER2 gene plays an important role in breast cancer, and overexpression or amplification of HER2 can lead to excess HER2 protein on the surface of breast cancer cells, leading to uncontrolled cell growth, tumor development and progression. To evaluate the correlation between the expression of the five hub genes (RIN2, THBS1, IL1R1, RAB31, COL11A1) and HER2 gene expression, we analyzed the gene expression data from TCGA-BRCA ( Figure 8D ) and GSE47994 ( Figure 8E ). Pearson correlation coefficient analysis revealed a significant positive correlation between the expression levels of these genes and HER2 gene expression (p < 0.05), with the exception of THBS1 and IL1R1 in GSE47994.