Identification and validation of a new prognostic signature based on cancer-associated fibroblast-driven genes in breast cancer

Abstract

BACKGROUND

Breast cancer (BC), a leading malignant disease, affects women all over the world. Cancer associated fibroblasts (CAFs) stimulate epithelial-mesenchymal transition, and induce chemoresistance and immunosuppression.

AIM

To establish a CAFs-associated prognostic signature to improve BC patient outcome estimation.

METHODS

We retrieved the transcript profile and clinical data of 1072 BC samples from The Cancer Genome Atlas (TCGA) databases, and 3661 BC samples from the The Gene Expression Omnibus. CAFs and immune cell infiltrations were quantified using CIBERSORT algorithm. CAF-associated gene identification was done by weighted gene co-expression network analysis. A CAF risk signature was established via univariate, least absolute shrinkage and selection operator regression, and multivariate Cox regression analyses. The receiver operating characteristic (ROC) and Kaplan-Meier curves were employed to evaluate the predictability of the model. Subsequently, a nomogram was developed with the risk score and patient clinical signature. Using Spearman's correlations analysis, the relationship between CAF risk score and gene set enrichment scores were examined. Patient samples were collected to validate gene expression by quantitative real-time polymerase chain reaction (qRT-PCR).

RESULTS

Employing an 8-gene (IL18, MYD88, GLIPR1, TNN, BHLHE41, DNAJB5, FKBP14, and XG) signature, we attempted to estimate BC patient prognosis. Based on our analysis, high-risk patients exhibited worse outcomes than low-risk patients. Multivariate analysis revealed the risk score as an independent indicator of BC patient prognosis. ROC analysis exhibited satisfactory nomogram predictability. The area under the curve showed 0.805 at 3 years, and 0.801 at 5 years in the TCGA cohort. We also demonstrated that a reduced CAF risk score was strongly associated with enhanced chemotherapeutic outcomes. CAF risk score was significantly correlated with most hallmark gene sets. Finally, the prognostic signature were further validated by qRT-PCR.

CONCLUSION

We introduced a newly-discovered CAFs-associated gene signature, which can be employed to estimate BC patient outcomes conveniently and accurately.

Key Words: Breast cancer; Prognosis; Gene signature; The Cancer Genome Atlas; The Gene Expression Omnibus

Core Tip: Breast cancer (BC), a leading malignant disease, affects women all over the world. Cancer associated fibroblasts (CAFs) stimulate epithelial-mesenchymal transition, and induce chemoresistance and immunosuppression. We introduced a newly-discovered CAFs-associated gene signature using IL18, MYD88, GLIPR1, TNN, BHLHE41, DNAJB5, FKBP14, and XG, which can be employed to estimate BC patient outcomes conveniently and accurately.

INTRODUCTION

Breast cancer (BC) is responsible for the most cancer-related deaths in women all over the world[1]. The Prediction Analysis of Microarray 50 (PAM50) classifies BC into five categories, namely, luminal A, luminal B, Her2-enriched, normal-like, and basal-like BC[2]. Despite this, patients with the same molecular and clinical signatures can exhibit vastly different outcomes as well as responses to chemotherapy and immunotherapy[3], suggesting the presence of unknown mediators that influence both the cancer prognosis and therapeutic response. Classical assessments such as those involving tumor-node-metastasis (TNM) staging and pathology are subjective at best and are poor at predicting patient outcomes and therapeutic responses. Hence, it is critical to develop novel approaches that are able to accurately estimate BC patient overall survival (OS) and thus improve outcomes.

Tumor pathogenesis and progression are modulated by both the cancer cells themselves and the tumor microenvironment (TME)[4]. The tumor stroma plays a critical role in enhancing tumor aggressiveness by modulating the deposition of the extracellular matrix (ECM) as well as cellular metabolism[5]. Of all TME stromal cells, cancer-associated fibroblasts (CAFs) have been the most widely examined. CAFs make up a majority of the stroma and are essential for tumor invasion, angiogenesis, and ECM remodeling, promoting cell-cell associations and the release of pro-invasive factors[6-8]. STAT3 induces breast cancer growth via ANGPTL4, MMP13 and STC1 secretion by cancer associated fibroblasts[9]. Cancer-associated fibroblasts facilitate premetastatic niche formation through lncRNA SNHG5-mediated angiogenesis and vascular permeability in breast cancer[10]. High expression of TGF-a, PKMYT1 and decreased SFRP1 and SFRP2 in the fibroblasts were associated with disease recurrence, invasive disease or decreased survival[11]. Hence, it is critical to explore the physiological properties of CAFs and their potential as therapeutic targets in BC.

Weighted gene co-expression network analysis (WGCNA) is an extensive bioinformatics algorithm that incorporates highly expressed genes and their associates into several gene modules to examine their association with the phenotype of interest[12]. In the past, WGCNA was successfully applied to identify CAF indicators[13,14]. To date, there have been no reports on WGCNA-based analysis of CAFs and stromal invasion in BC. Here, we reported that a CAF subtype functions as an indicator of BC patient prognosis, and we identified CAF-associated genes with prognostic significance. We also report an eight-gene signature using CAF-associated genes that can accurately estimate both OS and therapeutic response in BC patients. Our findings may highlight a novel and robust prognostic indicator for the personalized therapy of BC patients.

MATERIALS AND METHODS

Data sources and preprocessing

RNA-seq data [log2 (FPKM+1)], mutation annotation format file, and corresponding clinical information of 1072 breast invasive carcinoma patients in the The Cancer Genome Atlas (TCGA) database were obtained using Xena (https://xenabrowser.net). The Ensemble Gene was transformed into gene symbols to retrieve transcripts via the annotation file (hg38, gencode.v22.annotation.gene.probemap) of the GENCODE database[15]. Overall, 1072 primary BC samples and 99 adjoining healthy breast tissues were obtained from the TCGA. Table 1 summarizes in detail the clinical data of these primary BC patients. OS is defined as the time from diagnosis to death from any cause.

Table 1 Clinical and pathological characteristics of patients with breast cancer in The Cancer Genome Atlas cohort.

Variables	n	%
Age, yr
< 60	570	53.17
≥ 60	571	46.74
NA	1	0.09
Pathologic-T
T1	274	25.56
T2	621	57.93
T3	133	12.41
T4	40	3.73
TX	4	0.37
Pathologic-N
N0	502	46.83
N1	355	33.12
N2	118	11.01
N3	76	7.08
NA	21	1.96
Pathologic-M
M0	896	83.58
M1	22	2.05
NA	154	14.37
Stage
Stage I	176	16.41
Stage II	608	56.72
Stage III	246	22.95
Stage IV	20	1.87
NA	22	2.05
PAM50
Basal	136	12.69
Luminal-A	413	38.53
Luminal-B	185	17.25
Her2	66	6.16
Normal	22	2.05
NA	250	23.32

OS information was available for 1050 patients. The GSE96058, GSE18728, and GSE21653 datasets were downloaded from the The Gene Expression Omnibus (GEO) database (www.ncbi.nlm.gov/geo). The GSE96058 and GSE21653 datasets served as the validation cohorts.

The CAF-related gene expression signature (reference genomes) was downloaded from previous literature. This clearly separated CAFs into five subtypes[16]. Only the gene expression values greater than 1 were included. Supplementary Table 1 lists the CAF-related expression signature.

Calculation of CAF and microenvironment cell abundance

The Cell Type Identification by Estimating Relative Subsets of Known RNA Transcripts (CIBERSORT) algorithm was employed for the enrichment estimation of the five CAF and immune cell types for BC and normal samples from the TCGA dataset[17]. The ESTIMATE algorithm computed the BC patient stromal and immune scores using the bulk gene expression data[18]. Kruskal-Wallis tests were used to compare differences in CAF and immune cell proportions, as well as stromal and immune scores in different cancer stages and normal tissues.

Using univariate analysis, we assessed the association between the five CAF phenotypes and OS. The correlations between prognostic CAF subtype, stromal and immune scores were computed via Pearson’s correlation coefficients.

Construction of a weighted gene co-expression network and identification of CAF-associated genes

The WGCNA package was employed to construct the weighted gene co-expression network (WGCN)[12] in R software, based on mRNA expression profiling. Hierarchical clustering was used to obtain modules, with each module harboring a minimum of 30 genes (min Module Size = 30). We computed the eigengenes, clustered the modules based on hierarchy, and combined comparable modules (abline = 0.25). Pearson’s correlation coefficients were calculated to assess the relationship between the modules and patient clinical signatures. The gene significance (GS) was described as the log10 transformation of the matched P value (GS = lg P) of the association between the gene levels and matched clinical data. To further evaluate this module, the module eigengenes (ME) dissimilarity was computed via the module Eigengenes function in the R WGCNA package. This was described as the primary portion of a given module and was deemed as the gene expression signature representative in a module. The module membership (MM) gene values in modules were computed by evaluating associations between gene signatures and MEs. Hub genes were chosen for MM > one-third difference between the maximum and minimum MM, and GS > one-third difference between the maximum and minimum GS.

Functional enrichment analysis

Gene ontology (GO), including cellular component, biological process, and molecular function (MF), and Kyoto Encyclopedia of Genes and Genomes (KEGG) network enrichment analyses were conducted on CAF-related gene hub genes and was visualized using the clusterProfile package[19] in the R software. Adjusted P < 0.05 was set as the significance threshold.

Construction and verification of a CAFrelated gene prognostic signature CAFrelated in PDF format

To build a prognosis-related signature, univariate analysis was employed to screen for OS-associated genes with P values < 0.05. Subsequently, the marked prognostic genes underwent filtration using least absolute shrinkage and selection operator (LASSO)-penalized Cox regression analysis. Only genes with nonzero coefficients in the LASSO model were selected for additional multivariate analysis. Next, the prognostic risk model was generated based on the following: Risk score (RS, CAFPS) = β gene1 × expr gene1 + β gene2 × expr gene2 + ... + β gene n × expr gene n. In this formula, β represented the correlation coefficient in multivariate analysis, and expr represented mRNA expression. The patients were then separated into a high- (HR) or low-risk (LR) cohort, based on the median RS. The same formula was employed for the RSs calculation in the GEO and TCGA datasets.

The Kaplan-Meier log-rank test and the time-dependent receiver operating characteristic (ROC) curve analysis were employed to confirm the eight-gene signature performance. The RS distribution and OS status scatter plots, as well as the expression heatmaps of the eight CAF-related risk genes between the HR and LR cohorts, were constructed to assess the prognostic signature, and the stand-alone prognostic analysis was performed to assess the independent nature of the signature in predicting OS.

Development of nomogram

To translate the prognostic significance of the eight-gene signature into clinical application, we established a nomogram involving the RS and matched clinical data of BC patients, which included age and pathological stage. In addition, the calibration curves were generated to evaluate the agreement between the true and estimated OS.

Estimation of patient therapeutic response, according to the 8-gene signature

To further confirm the predictability of CAFPS in estimating patient chemotherapeutic response, GSE18728 was employed. Differences in CAFPS between responder and non-responder cohorts were assessed via the Wilcoxon test.

Enrichment analyses

Furthermore, the “h.all.v7.4. symbols” gene sets were derived from the Molecular Signatures Database (MSigDB)[20]. The ssGSEA was employed to compute the hallmark gene set enrichment scores[21]. Spearman correlation analysis was used to assess the relationships between the CAFPS and gene set enrichment scores.

Patients and tissue samples

All 80 breast cancer tissues and paired adjacent normal tissues were obtained from patients treated at the first hospital of qinhuangdao between January 2022 and December 2022. The paired adjacent normal tissues were dissected by the surgeons 5-cm away from the tumor edge. Tissue samples were stored at liquid nitrogen until total RNA was extracted. This study was approved by the Ethics Committee of the first hospital of qinhuangdao, and the informed consent forms were obtained from all the patients.

Quantitative real-time polymerase chain reaction

Total RNA was isolated from the BC and matched adjacent tissues using TRIzol (Thermo Fisher Scientific). The Nanodrop 2000 (Thermo Fisher Scientific) was used to detect the purity and concentration of the total RNA. According to the manufacturer's protocol, quantitative real-time polymerase chain reaction (qRT-PCR) was performed by using Taq Pro universal SYBR qPCR Master Mix (Vazyme, Nanjing, China).The primers used were as follows: FKBP14 forward primer 5'-AGCTGATCAACATCGGCAAT-3' and reverse primer 5'- CCAAGAACCCCCATTATTGA -3', βactin forward primer 5'CATGTACGTTGCTATCCAGGC3' and reverse primer 5'CTCCTTAATGTCACGCACGAT-3'. The relative levels of FKBP14 was evaluated by the 2^-ΔΔCT method using βactin as the control.

Statistical analysis

All data analyses were performed in R version 3.6.4 (https://www.r-project.org) and associated packages. The Wilcoxon test was used to examine the CAFPS-clinicopathological properties link. Kaplan-Meier curves and log-rank tests were used for evaluating OS using the survival and survminer R packages.

RESULTS

Clinical significance of the CAF proportion in BC

The proportions of five CAF and immune cell types in BC and normal samples from the TCGA database were obtained using the CIBERSORT algorithm. We observed marked differences in the proportions of the five types of CAFs (Figure 1), CD4 memory resting T cells, CD4 memory activated T cells, regulatory T cells, resting natural killer (NK) cells, and the stromal score between different cancer stages and normal tissues (Supplementary Figure 1). Based on the median value of the CAF proportion, patients were divided into high and low CAF groups. The prognosis of patients in different groups was compared (Figure 2A-F). Patients with lower proportions of the CAF S3 subset were found to have better OS (P = 0.015; Figure 2D). Figure 2G and H illustrates a positive correlation between the CAF-S3 proportion and the immune (r = 0.57, P < 2.2e-16) and stromal scores (r = 0.33, P < 2.2e-16).

Figure 1 Cell fractions of cancer-associated fibroblasts in different cancer stages and normal tissues. A: Cancer associated fibroblasts (CAFs); B: CAFs mesenchymal stem cell (MSC) inflammatory CAF (iCAF)-like s1; C: CAFs MSC iCAF-like s2; D: CAF Transitioning s3; E: CAFs myCAF-like s4; F: CAFs myCAF-like s5.

Figure 2 Kaplan-Meier curves of overall survival in relation to cancer associated fibroblasts, and the relationship between the cancer associated fibroblasts transitioning s3 and stromal and immune scores. A: Cancer associated fibroblasts (CAFs); B: CAFs mesenchymal stem cell (MSC) inflammatory CAF (iCAF)-like s1; C: CAFs MSC iCAF-like s2; D: CAFs Transitioning s3; E: CAFs myCAF like s4; F: CAFs myCAF like s5; G-H: The relationship between the CAFs transitioning s3 and the immune and stromal scores.

Use of WGCNA to develop a principal-gene module related to CAF-S3 cells

To better represent the scale-free topology of the co-expression axis in the WGCNA, we chose the soft-threshold β = 6 (Figure 3A). Ten co-expressed gene modules were clustered (Figure 3B), of which the black, pink, and red modules seen in Figure 3C were most relevant to CAF-S3 cells (Figure 3C). Of these, 921 genes [the black module (427 genes), pink module (255 genes), and red module (239 genes)] were selected as hub CAF-related genes (Figure 3D-F).

Figure 3 Co-expression pathways generated by weighted gene co-expression network analysis. A: A soft-thresholding power (β) of 6 was chosen, following the scale-free topology criterion; B: Clustering dendrograms depicting genes with comparable expression signatures were grouped into co-expression modules; C: Module-trait associations illustrating the relationships between individual gene module eigengenes and matching phenotype; D-F: Scatter plots of the module membership and gene significance of individual genes in the black, pink, and red modules.

Functional enrichment analysis

GO analysis showed that CAFrelated genes were mainly enriched in biological processes related to the ECM and cell adhesion, including focal adhesion, positive regulation of cell adhesion, extracellular structure organization, and ECM organization (Figure 4). KEGG pathway analysis indicated marked enrichment of the genes in pathways associated with cell adhesion (Figure 4).

Figure 4 Gene ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses. The gene ontology analysis of enriched biological process, cellular component, and molecular function terms, and the Kyoto Encyclopedia of Genes and Genomes analysis of hub genes. KEGG: Kyoto Encyclopedia of Genes and Genomes; MF: Molecular function; CC: Cellular component; BP: Biological process.

Screening of prognostic genes and generation of a CAF prognostic gene signature

Only gene expression values greater than 1 were included in this analysis. Univariate analysis showed that 72 CAF-related genes were strongly associated with OS. Further LASSO and multivariate analyses showed that eight CAF-associated genes (IL18, MYD88, GLIPR1, TNN, BHLHE41, DNAJB5, FKBP14, and XG) were independent prognostic factors of BC patient OS (Figure 5A-C).

Figure 5 The mechanism of variable selection. A and B: The performance of least absolute shrinkage and selection operator (LASSO) analysis; C: Forest plot of multivariate analysis for the establishment of the prognostic signature.

A risk score (RS) system was established using the expression levels of these eight genes and their corresponding coefficients obtained from the multivariate Cox regression analysis. We then calculated the RS using the formula: RS [CAF prognostic gene signature (CAFPS)] = (-0.168 × IL18 levels) + (-0.420 × MYD88 levels) + (-0.338 × GLIPR1 levels) + (-0.236 × TNN levels) +(-0.129 × BHLHE41 levels) + (-0.577 × DNAJB5 levels) + (0.695 × FKBP14 levels) + (0.573 × XG levels).

Evaluation of biomarkers

Genetic alterations in the biomarkers were assessed using data from the TCGA database. Most of the eight biomarkers exhibited high CNV frequencies in BC (IL18: 37.26%, MYD88: 22.05%, GLIPR1: 17.11%, TNN: 27.76%, BHLHE41: 16.16%, DNAJB5:14.07%, FKBP14: 15.21%, and XG: 14.83%; Figure 6).

Figure 6 Mutation and copy number alteration analysis of the hub genes. CAN: Copy number alteration.

Predictive accuracy of the eight-gene prognostic signature

Using the median RS, we separated BC patients into either low- (LR) or high-risk (HR) groups to evaluate the predictive accuracy of the eight-gene prognostic signature. This showed that the LR group had significantly better OS, compared with the HR group (HR 2.951, 95%CI 2.220-3.924, P < 0.0001; Figure 7A). The AUCs of the prognostic signature showed that it was effective for predicting the OS of patients with BC (3-year AUC = 0.698, 5-year AUC = 0.719, 8-year AUC = 0.752) (Figure 7B). We generated risk curves and scatter plots to display the RS and survival status of individual BC patients, finding that poorer outcome was strongly related to an elevated RS (Figure 7C and D). Furthermore, a heatmap depicting the expression patterns of the risk-associated genes in both groups showed that FKBP14 and XG levels were significantly increased in HR patients, while the levels of IL18, MYD88, GLIPR1, TNN, BHLHE41, and DNAJB5 were elevated in LR patients (Figure 7E). Multivariate analysis indicated that the CAFPS was an independent prognostic factor of BC patient prognosis after adjustment for other clinical factors (HR 2.466, 95%CI 1.843-3.298, P < 0.001) (Figure 8).

Figure 7 Evaluation of the prognostic signature using the The Cancer Genome Atlas database. A: Kaplan-Meier plots of overall survival between the high- and low-risk cohorts; B: Time-dependent receiver operating characteristic curve analysis; C: Risk score distribution; D: Survival distribution; E: Expression heatmap of the 8 cancer associated fibroblast-associated risk genes.

Figure 8 Assessment of the independent prognostic value of the prognostic signature using the The Cancer Genome Atlas database. A: Univariate analyses of the signature and clinical features; B: Multivariate analyses of the signature and clinical features.

The eight-gene prognostic signature was then validated using the GSE96058 and GSE21653 datasets. Kaplan-Meier survival curves confirmed that the HR cohort experienced worse OS relative to the LR cohort (GSE96085 cohort: P < 0.0001; GSE21653 cohort: P = 0.01) (Figure 9A and B).

Figure 9 Validation of the prognostic signature. A: Kaplan-Meier plots of overall survival of high- and low-risk cohorts in the GSE96058 database; B: Kaplan-Meier plots of overall survival of high- and low-risk cohorts in the GSE21653 database.

We then investigated whether RS could predict outcomes in various subgroups, including age ≥ 60, age < 60, phase I-II, phase III-IV, TNBC, and non-TNBC. This showed that all subgroups were accurately estimated by RS. Stratification analyses in the TCGA indicated worse OS for HR patients in each stratum, including age, stage, and molecular subtype, compared with LR patients. These data indicated that the CAFPS was a robust and independent predictor of OS in different populations (Supplementary Figure 2).

Correlation between the CAFPS and clinicopathological factors

To further assess whether CAFrelated genes participated in the development of BC, we assessed the associations between CAFPS and patient clinicopathological factors. This showed that the CAFPS was significantly correlated with pathological T stage, pathological M stage, and molecular subtype (Figure 10). The higher the T stage, the higher the corresponding RS. Patients with metastasis and Her2-enrichment had higher RS values.

Figure 10  Association between CAFPS and clinicopathological factors. A: Pathologic-T; B: Pathologic-N; C: Pathologic-M; D: Pathologic-stage; E: Immune subtype; F: PAM50.

Generation of the CAFPS-based nomogram

A nomogram was generated using the eight-gene CAFPS and several patient clinical characteristics, namely, age and pathological status, to predict the 3 and 5-year OS of BC patients (Figure 11A). All patients were given a point each for individual prognostic parameters, and a higher total score indicated worse prognosis. We next assessed the predictive efficacy of the nomogram using time-dependent ROC curve analysis. This showed that the AUCs of the nomogram were 0.805 and 0.801 for the 3- and 5-year OS rate predictions, respectively (Figure 11B and C). Moreover, higher AUCs were seen compared with the pathological stage. These results demonstrated the satisfactory predictive efficacy of the nomogram in determining the outcomes of BC patients (Figure 11D and E).

Figure 11  Nomogram generation for the prediction of 3- and 5-year survival. A: Development of a nomogram for estimating 1-, 2-, 3-and 5-year survival probabilities of breast cancer patients using the The Cancer Genome Atlas (TCGA) database; B and C: The 3- and 5-year area under the curves of nomogram, stage, CAFPS using the TCGA database; D and E: The 3- and 5-year calibration curves for predicting overall survival.

Prognostic value of CAFPS in relation to drug response

Patients in the GSE18728 dataset underwent chemotherapy. Of these, 11 responded while 17 did not. The CAFPS scores of the patients who responded were significantly lower compared to those of the non-responding patients (Wilcoxon test, P = 0.048, Figure 12A). Moreover, patients with reduced CAFPS scores responded better than those with elevated scores (Figure 12B). These findings indicate that the CAFPS was effective for estimating patient response to chemotherapy. This also explained the worse outcome and rapid cancer progression of patients with elevated CAFPS scores who underwent chemotherapy.

Figure 12  Prediction of chemotherapy response. A: The cancer associated fibroblast risk score of estimated chemotherapy-responders and non-responders in GSE18728; B: Distributions of responders and non-responders in high- and low-risk groups in GSE18728.

Evaluation of the prognostic gene signature, immune checkpoints, and hallmark gene networks

We next assessed differences in the levels of immune checkpoint molecules between the HR and LR cohorts. Relative to the HR cohort, the levels of PDCD1 (PD1) (P < 2.22e-16), CD274 (PD-L1) (P < 2.22e-16), CTLA4 (P < 2.22e-16), LAG3 (2.7e-13), IDO1 (P < 2.22e-16), and TIGIT (P < 2.22e-16) were significantly raised in the LR cohort (Figure 13). Another interesting finding was that almost all of the hallmark gene networks were significantly associated with CAFPS (Figure 14A). Moreover, the single-sample gene set enrichment analysis (ssGSEA) showed that the CAFPS RS was negatively associated with the TGF-BETA (r = -0.091, P = 0.0032), IL-6-JAK-STAT3 (r = -0.37, P < 2.2e-16), and TNFA-SIGNALING-VIA-NFKB (r =-0.34, P < 2.2e-16) signaling pathways (Figure 14B-D).

Figure 13  Immune checkpoint gene expression in both risk cohorts in The Cancer Genome Atlas cohort. A: PDCD1; B: CD274; C: CTLA4; D: LAG3; E: IDO1; F: TIGIT.

Figure 14  Correlation between CAFPS and signaling pathways. A: Spearman’s correlation analyses revealing that CAFPS is significantly correlated with most hallmark gene sets; B: The ssGSEA data revealing that the cancer associated fibroblast (CAF) risk score is directly associated with TGF-BETA signaling pathways enrichment scores using The Cancer Genome Atlas (TCGA) database; C: The ssGSEA data revealing that the CAF risk score is directly associated with IL-6-JAK-STAT3 signaling pathways enrichment scores using the TCGA database; D: The ssGSEA data revealing that the CAF risk score is directly associated with The TNFA-SIGNALING-VIA-NFKB signaling pathways enrichment scores using the TCGA database.

External experimental validation of prognostic signature

We found that the expression levels of FKBP14 was significantly higher in breast cancer samples than in the paired adjacent normal tissues. Patients were divided into high expression group and low expression group based on the median values of FKBP14 expression. Among patients with high FKBP14 expression, the proportion of patients with tumor size > 2 cm, lymph node metastasis and TNM stage III was higher than that of patients with low FKBP14 expression (P < 0.05) (Table 2).

Table 2 Relationship between expression level of FKBP14 and clinical features in breast cancer.

Variables	n	Low expression of FKBP14 (n = 40)	High expression of FKBP14 (n = 40)	χ2	P value
Age (yr)				0.056	0.813
≤ 60	53	27 (67.5)	26 (65)
> 60	27	13 (32.5)	14 (35)
Tumor size				10.769	0.002
≤ 2 cm	28	21 (52.5)	7 (17.5)
> 2 cm	52	19 (47.5)	33 (82.5)
Lymph node metastasis				25.208	< 0.001
No	48	35 (87.5)	13 (32.5)
Yes	32	5 (12.5)	27 (67.5)
Stage				11.250	0.001
I-II	64	38 (95)	26 (65)
III	16	2 (5)	14 (35)
Subtype				3.660	0.056
Triple negative	17	12 (30)	5 (12.5)
Non-triple negative	63	28 (70)	35 (87.5)

DISCUSSION

Current methods are relatively ineffective for predicting the outcomes of BC patients, likely due to complications related to phenotypes and underlying mechanisms. However, advances in bioinformatics have allowed the identification of new biomarkers to aid BC diagnosis and the estimation of patient outcomes[22,23]. CAFs, non-neoplastic stromal cells found within the TME, aid in tumor formation and metastasis[24,25]. Multiple studies have identified CAFs as indicators of drug resistance and worse outcome in a range of tumors[26-28]. Moreover, it was revealed that factors secreted by CAFs can transmit signals to tumor cells that promote drug resistance[29]. Furthermore, CAFs also release cytokines, exosomes, and growth factors, which influence components of the TME, thereby modulating the therapeutic response. Emerging evidence has revealed a strong association between the number and activity of CAFs and patient prognosis[13,14,28]. Hence, targeting CAFs may be beneficial to anti-cancer therapy.

Prior investigations have identified fibroblasts based on their expression of FAP and αSMA[30]. However, fibroblasts are distinct from cancer cells and have thus not been further stratified. Single-cell RNA sequencing (scRNA-seq) can examine the entire cellular signature of tumor cells[31,32]. Recent CAF examinations recognized two polarized states regulated by ECM formation or inflammatory secretomes[33-35].

Previous reports indicated that cancer-associated fibroblast genes expression can estimate breast cancer patient prognosis and predict the responses to immunotherapy[36-38]. Wang et al[36] reported that signature based on cancer-associated fibroblast genes could divide patients into low- and high-risk groups, accompanied by different OS, clinical features, and immune infiltration characteristics. Huang et al[37] identified various heterogeneous CAF cell populations in breast cancer patients. An CAF-associated gene signature was developed to predict the responses to immunotherapeutics. The five-gene prognostic CAF signature presented by Xu et al[38] was effective in estimating clinical immunotherapy response. In this study, we confirmed that the proportion of CAF-S3 cells was closely associated with patient prognosis, based on the CAF-related gene expression signature obtained from previous reports[16]. Subsequently, we regarded the proportion of CAF-S3 cells as a clinical characteristic and developed a CAF-based gene signature using eight genes (IL18, MYD88, GLIPR1, TNN, BHLHE41, DNAJB5, FKBP14, and XG) using WGCNA, univariate, LASSO, and multivariate analyses to predict patient prognosis. This demonstrated that the LR cohort had markedly better OS compared with the HR cohort. AUC analysis of the prognostic signature demonstrated its accuracy in predicting the OS of patients with breast cancer. Furthermore, multivariate analysis showed that the signature was an independent prognostic factor. The CAFPS OS predictive efficacy was further validated using the GSE96058 and GSE21653 datasets from the GEO and subgroups of the TCGA databases. These findings indicate that the CAFPS was a reliable predictor of BC patient prognosis.

We then generated a nomogram to estimate BC patient outcomes. Based on our results, our model was highly effective in predicting BC patient prognosis. The results of the calibration curves demonstrated that the nomogram was able to accurately predict the prognosis of patients with BC. However, when compared with traditional clinical features, the prognostic signature showed even better potential for clinical application. According to the model, a reduced CAFPS RS was significantly related to enhanced chemotherapeutic response in BC patients. In addition, CAFPS was strongly correlated with most hallmark gene sets, explaining the rapid progression and worse outcomes of patients with elevated CAFPS values. Immune checkpoint molecules can also reflect the immune status of tumor microenvironments. In the present study, the levels of immune checkpoint molecules were found to be significantly elevated in the LR cohort, implying that patients with low CAFPS RS may respond better to immune checkpoint therapy.

For a proper evaluation of the prognostic signature, we then assessed individual biomarkers. IL-18 is a member of the IL-1 family of proinflammatory cytokines. IL-18 modulates Th1 and Th2 immune responses, as well as the activation of NK cells and macrophages[39]. It is reported that the IL-18 protein is strongly expressed in human BC and lymph node tissues, and is directly associated with lymph node metastasis in BC[40]. MYD88 serves a critical function in both innate and adaptive immune responses. It is also involved in crosstalk with the interleukin-1 and Toll-like receptor networks. MYD88 expression is high in metastatic tumors and has been suggested as a potential therapeutic target[41]. GLIPR1 is markedly up-regulated in BC and has been identified as an oncoprotein in BC. In BC cells, GLIPR1 expression is promoted by prolactin and estrogen, promoting cell proliferation[42]. TNN is abundantly expressed in human breast tumor extracts, although it is absent in adjoining normal tissue. Therefore, it is possible that TNN promotes tumor development by inducing BC cell migration[43]. BHLHE41 is a critical modulator of the invasive and metastatic phenotype in triple-negative BC. It is modulated by the p63 metastasis suppressor, and it suppresses TNBC aggressiveness by inhibiting hypoxia-inducible factor 1a (HIF-1a) and HIF-2a (HIFs)[44]. DNAJB5 belongs to the DnaJ/heat shock protein 40 (HSP40) family, and it partners with HSP70 chaperones. Stable DNAJB5 expression in cholangiocarcinoma cells that overexpress microRNA 21 re-sensitizes them to HSP90 inhibitors[45]. FKBP14 is a member of the FK506-binding protein (FKBP) family and is highly expressed in gastric cancer. Patients with elevated FKBP14 Levels typically have poor outcomes. FKBP14 is recognized as an oncogene in gastric cancer and has the potential as a marker for gastric cancer identification, invasion, and prognosis[46]. Relative to normal tissues, FKBP14 is ubiquitously expressed in ovarian cancer tissues where it was shown to be indicative of tumor size and grade in ovarian cancer tissues[47]. In our study, Among patients with high FKBP14 expression, the proportion of patients with tumor size > 2 cm, lymph node metastasis and TNM stage III was higher than that of patients with low FKBP14 expression. XG belongs to the CD99 family and its expression has been correlated with worse OS in Ewing's sarcoma patients with benign tumors, and it is speculated to regulate metastasis. Moreover, we further verified the abnormal expression of the risk genes (FKBP14) through experiments. Nevertheless, there is very limited information regarding the function of CAFs in BC. Hence, it is critical to conduct additional investigations into the eight CAF-related gene markers and their associated underlying mechanisms to better elucidate infiltration, metastasis, and drug resistance in BC patients. Our research has several limitations: (1) The size of our analyzed dataset was relatively small, and a larger patient population must be studied to validate our results; and (2) This paper did not go into detail regarding the mechanism associated with the eight identified genes.

CONCLUSION

Here, we generated a CAF-associated gene prognostic signature using IL18, MYD88, GLIPR1, TNN, BHLHE41, DNAJB5, FKBP14, and XG as prognostic markers of BC. This novel signature can be employed to estimate BC patient outcomes conveniently and accurately. These results provide a basis for further studies on the role of CAFs in BC.

ARTICLE HIGHLIGHTS

Research background

Breast cancer (BC), a leading malignant disease, affects women all over the world. It is still urgent to explore new biomarkers to estimate and enhance prognosis of BC patients.

Research motivation

The present study for the first time investigated the 8 cancer associated fibroblasts (CAFs)-associated genes as the potential biomarkers of the prognosis of patients using bioinformatics.

Research objectives

This study aims to establish a CAFs-associated prognostic signature to improve BC patient outcome estimation.

Research methods

We retrieved the transcript profile and clinical data of 1072 BC samples from The Cancer Genome Atlas (TCGA) databases, and 3661 BC samples from the The Gene Expression Omnibus. CAFs and immune cell infiltrations were quantified using CIBERSORT algorithm. CAF-associated gene identification was done by weighted gene co-expression network analysis. A CAF risk signature was established via univariate, LASSO regression, and multivariate Cox regression analyses. The receiver operating characteristic (ROC) and Kaplan-Meier curves were employed to evaluate the predictability of the model. Subsequently, a nomogram was developed with the risk score and patient clinical signature. Using Spearman's correlations analysis, the relationship between CAF risk score and gene set enrichment scores were examined.

Patient samples were collected to validate gene expression by quantitative real-time polymerase chain reaction (qRT-PCR).

Research results

Employing an 8-gene (IL18, MYD88, GLIPR1, TNN, BHLHE41, DNAJB5, FKBP14, and XG) signature, we attempted to estimate BC patient prognosis. Based on our analysis, high-risk patients exhibited worse outcomes than low-risk patients. Multivariate analysis revealed the risk score as an independent indicator of BC patient prognosis. ROC analysis exhibited satisfactory nomogram predictability. The AUC showed 0.805 at 3 years, and 0.801 at 5 years in the TCGA cohort. We also demonstrated that a reduced CAF risk score was strongly associated with enhanced chemotherapeutic outcomes. CAF risk score was significantly correlated with most hallmark gene sets. Finally, the prognostic signature were further validated by qRT-PCR.

Research conclusions

We introduced a newly-discovered CAFs-associated gene signature, which can be employed to estimate BC patient outcomes conveniently and accurately.

Research perspectives

The mechanisms of 8 CAFs-associated genes in BC is still unclear, which needs further confirmation through molecular biology and clinical experiments.

ACKNOWLEDGEMENTS

The authors thank TCGA and GEO for sharing the BC data.