Identification of Molecular Subtypes and Prognostic Features of Breast Cancer Based on TGF-β Signaling-related Genes

Abstract

Objectives:

The TGF-β signaling pathway is widely acknowledged for its role in various aspects of cancer progression, including cellular invasion, epithelial-mesenchymal transition, and immunosuppression. Immune checkpoint inhibitors (ICIs) and pharmacological agents that target TGF-β offer significant potential as therapeutic options for cancer. However, the specific role of TGF-β in prognostic assessment and treatment strategies for breast cancer (BC) remains unclear.

Methods:

The Cancer Genome Atlas (TCGA) database was utilized to develop a predictive model incorporating five TGF-β signaling-related genes (TSRGs). The GSE161529 dataset from the Gene Expression Omnibus was employed to conduct single-cell analyses aimed at further elucidating the characteristics of these TSRGs. Additionally, an unsupervised clustering algorithm was applied to categorize BC patients into two distinct groups based on the five TSRGs, with a focus on immune response and overall survival (OS). Further investigations were conducted to explore variations in pharmacotherapy and the tumor microenvironment across different patient cohorts and clusters.

Results:

The predictive model for BC identified five TSRGs: FUT8, IFNG, ID3, KLF10, and PARD6A. Single-cell analysis revealed that IFNG is predominantly expressed in CD8+ T cells. Consensus clustering effectively categorized BC patients into two distinct clusters, with cluster B demonstrating a longer OS and a more favorable prognosis. Immunological assessments indicated a higher presence of immune checkpoints and immune cells in cluster B, suggesting a greater likelihood of responsiveness to ICIs.

Conclusion:

The findings of this study highlight the potential of the TGF-β signaling pathway for prognostic classification and the development of personalized treatment strategies for BC patients, thereby enhancing our understanding of its significance in BC prognosis.

Introduction

Breast cancer (BC) represents the most prevalent malignancy among women, constituting approximately 30% of all cancer cases. ¹ The prognosis for individuals diagnosed with BC remains unfavorable, primarily due to the disease’s significant heterogeneity and propensity for early metastasis. Current therapeutic modalities, including surgical intervention, chemotherapy, radiotherapy, neoadjuvant therapy, targeted therapy, and immunotherapy, while effective in reducing mortality rates among BC patients, are insufficient in addressing the challenges posed by the high heterogeneity and early metastatic behavior of the disease.^2,3 Furthermore, existing classification systems, such as the TNM staging and PAM50 classification, fall short in accurately predicting patient outcomes. ⁴ Consequently, the identification and development of novel biomarkers to facilitate the reclassification of BC subtypes and to inform treatment strategies is essential for enhancing patient prognosis.

The tumor microenvironment (TME) comprises the extracellular matrix (ECM), stromal cells, and tumor-infiltrating immune cells (TIICs). ⁵ The type and functional mechanisms of TIICs are critical, as they not only affect the efficacy of anti-cancer immunotherapies such as immune checkpoint inhibitors (ICIs), but also serve as key determinants in predicting the clinical Research has established a robust correlation between the prognosis of triple-negative BC and the extent of T cell infiltration within the TME.^{6 -8} Additionally, there is increasing evidence that various signaling pathways modulate the composition and behavior of TIICs in the TME, with the transforming growth factor beta (TGF-β) signaling pathway emerging as particularly significant.^9,10

The TGF-β signaling pathway exhibits a complex dual role in cancer progression. ¹¹ In the early stages of cancer, TGF-β acts as a tumor suppressor by promoting cellular differentiation and inhibiting proliferation. Conversely, in advanced stages, TGF-β can function as a tumor promoter by facilitating epithelial-to-mesenchymal transition (EMT), thereby enhancing cancer initiation and enabling immune evasion. ¹² TGF-β also promotes ferroptosis, a form of regulated cell death that inhibits cancer cell proliferation. Specifically, TGF-β1 activates Smad3, which in turn suppresses the expression of the cystine/glutamate antiporter system Xc- (xCT), leading to decreased intracellular glutathione levels and increased lipid peroxidation, ultimately triggering cellular ferroptosis. ¹³ Moreover, excessive TGF-β activity within the TME can impair immune responses by inhibiting various immune cell types, including cytotoxic T cells and natural killer (NK) cells.^12,14 The interplay between TGF-β and the programmed cell death-1 receptor/programmed cell death ligand 1 (PD-1/PD-L1) pathways has been shown to be both independent and complementary, suggesting that TGF-β may serve as a predictive marker for the therapeutic efficacy of ICIs targeting PD-1/PD-L1.^{15 -17} Although TGF-β has been proposed as a potential therapeutic target for BC treatment in recent years,^18,19 further investigation into the intrinsic relationship between TGF-β and the TME is warranted to better assess the most effective treatment modalities for BC patients and to improve their overall prognosis.

This research established a predictive model for breast cancer that focuses on the TGF-β signaling pathway. A comprehensive analysis revealed a significant association between TGF-β and immune cells as well as cytokines within the tumor microenvironment (TME). The prognostic characteristics associated with TGF-β can be utilized to forecast and delineate prognostic outcomes and immune response status in breast cancer patients, as demonstrated by survival and drug susceptibility assessments. Additionally, cluster analysis based on five TGF-β signaling-related genes (TSRGs) enables the classification of breast cancer patients into two distinct molecular subtypes, which may enhance future clinical patient stratification and therapeutic strategies.

Materials and Methods

Data sets collection

The transcriptome expression matrix and clinical data from the TCGA-BRCA cohort, sourced from the Cancer Genome Atlas (TCGA) database, served as the training sets for this study. The TCGA-BRCA cohort comprised 1,113 BC samples alongside 113 normal breast tissue samples. For validation purposes, data sets GSE20685, GSE42568, and GSE58812 were obtained from the Gene Expression Omnibus (GEO) database. All three validation sets included overall survival (OS) data, with GSE42568 additionally providing recurrence-free survival (RFS) data and GSE58812 offering metastasis-free survival (MFS) data. Ultimately, a training set was established comprising 976 patients from the TCGA-BRCA cohort, while the validation sets included 327 patients from the GSE20685 cohort, 104 patients from the GSE42568 cohort, and 106 patients from the GSE58812 cohort. Furthermore, 225 TSRGs were extracted from a previously published article. ²⁰

Mutation analysis and immune cell infiltration analysis

Somatic mutation data for the TCGA-BRCA cohort were retrieved from UCSC Xena, and the somatic mutation frequencies of TSRGs were visualized using the R package “maftools” in the form of a waterfall plot. ²¹ The single-sample Gene Set Enrichment Analysis (ssGSEA) method was employed to investigate variations in immune cell infiltration between normal breast tissue and BC tissue, with specific details regarding 28 immune cell types obtained from supplementary materials in a published article. ³ Additionally, differences in the TME between normal and BC groups were evaluated using the ESTIMATE algorithm to derive stromal and estimate scores.

Construction of TGF-β risk score model

Initially, prognostic genes associated with OS were identified through univariate Cox regression analysis, applying a significance threshold of P < .05 within the TCGA-BRCA cohort. These OS-related genes were subsequently intersected with the 225 TSRGs, and the least absolute shrinkage and selection operator (LASSO) Cox regression analysis was utilized to minimize redundant genes and mitigate the risk of model overfitting. ²² Ultimately, multivariate Cox regression analysis was conducted to select the most effective prognostic model and to calculate coefficients for each candidate gene. The TGF-β risk score was determined as the product of the expression levels of each candidate gene and their corresponding coefficients. Patients were categorized into high- or low-risk groups based on their median TGF-β risk score. Additionally, boxplots were employed to illustrate the differences in candidate gene expression between the two groups, while Kaplan-Meier survival analysis was utilized to assess the impact of candidate gene expression on the prognosis of BC patients.

Validation of TGF-β risk score model

The three previously mentioned datasets were utilized for validation purposes to further evaluate the accuracy of the TGF-risk score model. Patients were categorized into high- and low-risk groups based on the median risk scores within each validation set, with risk scores derived from the expression levels of selected candidate genes. Subsequently, the survival outcomes of patients in these risk categories were analyzed using Kaplan-Meier survival analysis. Given the constraints posed by variations in existing data, receiver operating characteristic (ROC) curves were employed in the GSE58812 validation set to assess the model's predictive accuracy for patient prognosis at 1, 3, and 5 years, while ROC curves from other training and validation sets were utilized to evaluate the model’s prognostic accuracy at 3, 5, and 10 years.

Construction and validation of the TGF-β risk score related-clinical nomogram

In the construction and validation of a TGF-β risk score-related clinical nomogram, independent prognostic factors were identified by examining the clinicopathological characteristics and TGF-β risk scores of the TCGA-BRCA cohort through univariate and multivariate Cox regression analyses. The “rms” and “regplot” R packages were subsequently employed to develop the nomogram based on these findings. ²³ To validate the predictive capability of the nomogram, calibration curves for 1, 3, and 5 years, ROC curves for 1, 2, and 5 years, and decision curve analysis (DCA) were conducted.

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) was performed using the R package “clusterProfiler,” which included gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for patients categorized into high- and low-risk groups. The results were visually represented through bubble plots and mountain plots, respectively.

Immune landscape estimation in BRCA TME

Given the close association of emerging immunotherapies with immune checkpoints, the relationship between candidate genes and six immune checkpoints (PDCD1, CTLA4, IDO1, TIGIT, LAG3, and PVR) was assessed. Segmented violin plots illustrated the expression differences of these immune checkpoints across various risk groups and datasets. Additionally, to explore the role of TSRGs in the TME of BC patients, the Tumor Immune Estimation Resource (TIMER) database was utilized to analyze the correlation between the expression of these genes and the aforementioned immune checkpoints. Furthermore, transcriptome data from 976 BC patients were uploaded to the tumor immune dysfunction and exclusion (TIDE) website to evaluate responses to immune checkpoint blockade (ICB) in high- and low-risk patients.

Predicting the clinical sensitivity of BC patients to ICIs and chemotherapy drugs

As the immunophenotype score (IPS, ranging from 0 to 10) from the Cancer Immunome Atlas database correlates with immunogenicity and can be used to predict sensitivity to ICIs treatment in patients with BC, IPS values from patients in the high- and low-risk groups were compared. Additionally, to forecast the responsiveness of BC patients to chemotherapeutic agents, drug susceptibility data were sourced from the Genomics of Drug Sensitivity in Cancer (GDSC) database, and the 50% inhibitory concentrations (IC50) of various drugs were predicted using the R package oncoPredict. ²⁴

Prediction of BRCA molecular subtypes by cluster analysis

Clustering analysis was conducted utilizing the R package “ConsensusClusterPlus” to predict TGF-β related molecular subtypes in BC patients, based on candidate genes identified in previous studies. ²⁵ Subsequently, principal component analysis (PCA) was employed to validate the accuracy of the classifications, while Kaplan-Meier curves were utilized to evaluate the prognostic outcomes of the subgroups. Additionally, relevant immunoassays were performed to investigate the differences in the TME among the identified subgroups.

Single-cell RNA-seq data acquisition and processing

For the acquisition and processing of scRNA-seq data, the GSE161529 dataset from the GEO database was utilized to accurately assess the expression levels of five TSRGs at the cellular level. The “Seurat” software was employed for quality control, excluding cells with “nFeature” values below 200 and “percent.mt” values exceeding 20%. The data normalization was achieved using the “LogNormalize” method. Following normalization, PCA was applied for dimensionality reduction, and uniform manifold approximation and projection (UMAP) was utilized for unsupervised cell clustering. Finally, cell annotation was performed using the SingleR package.

Statistical analysis

Statistical analyses and data visualizations were executed using R software (version 4.2.0). Both univariate and multivariate Cox regression analyses were conducted to identify independent prognostic factors, with a two-tailed p-value of less than .05 deemed statistically significant.

Results

Identification of TSRGs in BC patients

First, the proportion of 225 TSRGs with somatic mutations in the TCGA-BRCA cohort was assessed, and a waterfall plot was used to display the top 20 mutant genes, including TP53, CDH1, and others (Figure 1(a)). Furthermore, Figure 1(b) depicted the difference in TIICs in the TME between normal breast and BC tissue, with mast cells, eosinophils, natural killer T cells, and effector memory CD8 T cells being higher in normal tissues, while activated CD4 T cells and CD56 dim natural killer cells were higher in BC tissues. Lower stromal and ESTIMATE scores in BC tissues indicate a worse prognosis in BC patients (Figure 1(c)).

Figure 1.

(a) Somatic mutation frequency of TGF-β in the TCGA-BRCA cohort, (b) Differences in the degree of invasion of 28 immune cells between normal and cancerous tissues in BC patients, (c) Variance of stromal score and ESTIMATE score in the different groups, (d) The LASSO Cox regression model was utilized to identify TSRGs, and (e) Select range of the optimal parameter (lambda) in the LASSO Cox regression model.

Construction of TGF-β risk score model for BC patients

An optimal prognostic model was established through the integration of univariate COX regression analysis, LASSO regression analysis, and multivariate COX regression analysis (Figure 1(d) and (e)). This model comprises five TSRGs: FUT8, IFNG, ID3, KLF10, and PARD6A. Among these, ID3 and KLF10 are identified as adverse prognostic factors, indicated by a hazard ratio (HR) greater than 1, while FUT8, IFNG, and PARD6A are classified as protective prognostic factors, as evidenced by an HR less than 1 (Figure 2(a)). The risk score is calculated using the formula: Risk score = −0.20464 * FUT8 - 0.65792 * IFNG + 0.21690 * ID3 + 0.37394 * KLF10 - 0.30315 * PARD6A. Based on the median risk score, subjects are categorized into high-risk and low-risk groups. Subsequent analysis revealed that ID3 and KLF10 exhibited higher expression levels in the high-risk group, whereas FUT8, IFNG, and PARD6A were more prominently expressed in the low-risk group (Figure 2(b)). The Kaplan-Meier survival curve analysis further assessed the independence of the five identified genes as predictive variables for BC. With the exception of FUT8 (P = .13), the other four genes were determined to be independent prognostic factors influencing BC outcomes (P < .05), as illustrated in Figure 2(c). Additionally, the differential expression of these five candidate genes in normal versus BC tissue is presented in Supplemental Figure S1 and S4.

Figure 2.

(a) The coefficient of the selected TSRGs, (b) Expression of 5 TSRGs in the high- and low-risk groups, and (c) The KM survival curves for 5 selected TSRGs based on expression level and OS.

Prognostic validation based on TGF-β risk score model

The same prognostic analysis was performed in the training (TCGA-BRCA) and three validation sets (GSE20685, GSE42568, and GSE58812) to further investigate the accuracy of the TGF-β risk score model in BC prediction. Initially, patients within each dataset were categorized into high- or low-risk groups based on their median risk scores (Figures 3(b) and 4(b)). The results aligned with expectations, indicating an increase in mortality among BC patients corresponding to elevated hazard scores (Figures 3(c) and 4(c)). Consequently, individuals classified in the high-risk group exhibited poorer prognoses compared to those in the low-risk group, with statistical significance observed across various datasets (TCGA-BRCA: P < .0001; GSE58812-OS: P = .003; GSE58812-MFS: P = .015; GSE20685: P = .0015; GSE42568-OS: P = .033; GSE42568-RFS: P = .033) (Figures 3(a) and (a)). Heat maps illustrating the distribution of the five TSRG genes across the high- and low-risk groups are presented in Figures 3(d) and 4(d). Additionally, the predictive capability of the TGF-β risk score model was further assessed through ROC analysis, demonstrating satisfactory performance across all cohorts (Figures 3(e) and 4(e)).

Figure 3.

Assessment and verification of the efficiency of the TGF-β risk score: (a) The KM analysis of the overall survival in the TCGA-BRCA, GSE58812-OS and GSE58812-MFS cohorts, (b-d) Risk scores, survival status, and expression distributions of TGF-β risk score in different cohorts, and (e) The ROC curves of the TGF-β risk score in predicting 3, 5, and 10-year survival state of BC patients.

Figure 4.

Additional validation set for the assessment and verification of the efficiency of the TGF-β risk score: (a) The KM analysis of the overall survival in the GSE20685-OS, GSE42568-OS and GSE42568-RFS cohorts, (b-d) Risk score, survival status, and expression distributions of TGF-β risk score in different cohorts, and (e) The ROC curves of the TGF-β risk score in predicting 3-, 5-, and 10-year OS in the GSE20685 cohort and 1,3,5 years survival state in the GSE42568 cohort.

Construction and validation of a nomogram for TGF-β risk score and clinical characteristics

The relationship between the TGF-β risk score and various clinicopathological factors within the TCGA-BRCA dataset is illustrated in Supplemental Figure S2. Initial assessments utilizing both univariate and multivariate Cox regression analyses aimed to determine whether the variables T, N, M, age, stage, and TGF-β risk score functioned as independent prognostic indicators within the TCGA-BRCA cohort. The analysis revealed that age, M classification, and the TGF-β risk score were indeed independent prognostic factors (Figures 5(a) and (b)). Consequently, a clinicopathological nomogram for BC was developed incorporating these three variables to predict OS for BC patients at 1, 3, and 5 years (Figure 5(c)). The calibration curves for 1, 3, and 5 years demonstrated a strong concordance with the 45-degree reference line, thereby affirming the reliability of the constructed nomogram (Figure 5(d)). Additionally, the area under the curve (AUC) for the ROC analyses at 1, 2, and 5 years for the TCGA cohort were recorded at 0.78, 0.72, and 0.75, respectively (Figure 5(e)). These values surpassed the AUCs for the TGF-β risk score model alone, which were 0.66, 0.67, and 0.71 for the 3, 5, and 10-year intervals (Figure 3(e)). Furthermore, the DCA indicated that the nomogram provided a superior net benefit in prognostic prediction for BC patients compared to the use of clinicopathological factors in isolation (Figure 5(f)).

Figure 5.

Establishment of a prognosis-related nomogram model based on TGF-β risk score: (a) Univariate Cox regression analysis of TGF-β risk score and clinicopathological characteristic, (b) Multivariate Cox regression analysis of TGF-β risk score and clinicopathological characteristic, (c) Development of the nomogram based on M, age and TGF-β risk score, (d) The Calibration plots for the nomogram, (e) The ROC curves and AUC values demonstrated favorable competence of the nomogram in predicting 1-, 2-, and 5-year OS of BC patients, and (f) The decision curve to evaluate the clinical decision effectiveness of the nomogram against other separate clinical parameters.

Gene set enrichment analysis based on TGF-β risk score model

The results of the GSEA, conducted to elucidate the relationship between TGF-β and biological signaling pathways, are presented in Figure 6(a) and (b). Subsequent GO analysis revealed that the pathways exhibiting significant upregulation included the immunoglobulin complex, cilium assembly, and cilium organization. Conversely, the pathways that demonstrated significant downregulation encompassed the T cell receptor complex, plasma membrane signaling receptor complex, regulation of cytosolic calcium ion concentration, and cellular calcium ion homeostasis, among others. Furthermore, the KEGG analysis indicated that the high-risk group was notably enriched in pathways related to Herpes simplex virus 1 infection, chemokine signaling, neutrophil extracellular trap formation, neuroactive ligand-receptor interactions, cytokine-cytokine receptor interactions, and natural killer cell-mediated cytotoxicity.

Figure 6.

(a-b) Results of GO and KEGG in TCGA-BRCA cohort, (c) Expression level of PDCD1, CTLA4, IDO1, TIGIT, LAG3 and PVR in TCGA-BRCA, (d) Correlation chord plot between immune checkpoints and riskscore, (e) Correlation of five TSRGs with six immune cells, and (f) Differences in TIDE and Exclusion scores between high and low risk groups.

Immune landscape analysis in BRCA TME

The differences in TME between high-risk and low-risk groups were further explored. Firstly, except for PVR, the expression of immune checkpoints was higher in the low-risk group (Figure 6(c)), which was also verified in the GEO dataset (Supplemental Figure S3). The chord plot showed that these six immune checkpoints are closely related to the patient’s risk score (Figure 6(d)). Furthermore, Figure 6(e) depicts the link of five TSRGs to six immunological checkpoints, with IFNG showing a substantial association, indicating that it plays a key role in the immune cells of the patient’s TME. Simultaneously, IFNG was found to have a substantial link with six immunoinfiltrating cells: CD4 T cells, B cells, CD8 T cells, neutrophils, macrophages, and dendritic cells (Supplemental Figure S4(a)). The differences in 28 types of immune cell infiltration in patients in the high and low risk groups are shown in the Supplemental Figure S4(b). In addition, patients in the high-risk group had higher TIDE scores, suggesting that they were more likely to have immune escape, which was associated with a poor prognosis (Figure 6(f)).

Drug sensitivity analysis in BC patients

IPS values were used to predict how patients would respond to anti-CTLA-4 and anti-PD-1/PD-L1 treatments. The findings are consistent with our immune checkpoint analysis, which found that the low-risk group had higher IPS values, indicating that they may be more sensitive to ICIs treatment (Figure 7(a) and (b)). Furthermore, IC50 was used to estimate the response of BC patients to 198 chemotherapy drugs. The lower the IC50 value, the better the drug’s effect. Chemotherapy drugs with a good treatment response in both high-risk and low-risk patients were chosen, and the results are shown in Figure 7(c) and (d). These medications include: Buparlisib, Dihydrorotenone, PD0325901, Rapamycin, Dactinomycin, AZD8055, CDK9.1, BMS.754807, Eg5, Epirubicin, GNE.317, Sabutoclax, BI.2536, and Rapamycin.

Figure 7.

(a) Relative probability of response to CTLA-4 treatment in low-risk and high-risk groups, (b) Relative probability of response to PD-1/PD-L2 treatment in low-risk and high-risk groups, (c) A venn diagram of 14 drugs sensitive to BC patients, and (d) A bar plot of IC50 values for 14 drugs sensitive to BC patients.

Description of TGF-β subtypes in BC

The R package “ConsensusClusterPlus” was used to create consistent clustering of five genes for the TGF-β risk score model. When BC patients were divided into two clusters, there was good consistency and stability within subgroups (Figure 8(a)-(d)). The PCA results revealed that the two clusters of patients were clustered separately, confirming the reliability of the clustering results (Figure 8(e)). The survival analysis revealed that cluster A had a better prognosis than cluster B (P = .0081) (Figure 8(f)). A heat map also depicted the expression of five TSRGs between the two clusters (Figure 9(a)).

Figure 8.

(a) The CDF value of consensus index, (b) Relative change in area under CDF curve for k = 2–6, (c) The tracking plot for k = 2 to k = 6, (d) Consensus matrix for k = 2, (e) Principal component analysis of the total RNA expression profile, and (f) KM curve of the survival difference between cluster A and cluster B.

Figure 9.

(a) Heatmap of the five genes between the two clusters, (b) Variance of immune score, stromal score, and ESTIMATE score in the two clusters, (c) Boxplots were used to depict the discrepancies in the infiltration extent of 28 immune cells between the two clusters, (d) GSEA analysis between Cluster A and Cluster B, and (e) Expression level of PDCD1, CTLA4, IDO1, TIGIT, LAG3 and PVR in the cluster A and cluster B groups.

The immunoinfiltration landscape in TME in two subgroups was also examined. Cluster B had higher matrix, immunity, and ESTIMATE scores as determined by the ESTIMATE algorithm (Figure 9(b)). Furthermore, with the exception of CD56 dim natural killer cells, which were more abundant in cluster A, the majority of the remaining immune cells were expressed more strongly in cluster B (Figure 9(c)). In addition, the pathways that GSEA analysis showed differences between the two clusters were mainly antigen binding, immunoglobulin complex immunoglobulin complex,circulating (Figure 9(d)). Surprisingly, the expression of immune checkpoints was higher in Cluster B than in Cluster A, indicating that they may be more suitable for ICIs (Figure 9(e)).

Single-cell RNA-seq analysis

Following quality control, 40,467 cells were kept for further investigation. UMAP is further subdivided into 14 cell clusters consisting of six different cells: HSCs, Erythroid cells, CD8+T cells, CD4+T cells, Eosinophils, Monocytes (Figure 10(a) and (b)). IFNG is mainly distributed in CD8+T cells, while KLF10, ID3, PARD6A are mainly distributed in HSCs and Erythroid cells (Figure 10(c)).

Figure 10.

(a) The results of the reduced dimensionality clustering analysis are shown in the UMAP plot, (b) Annotate cells into 2 different types of cells, and (c) Expression distribution of FUT8, IFNG, ID3, KLF10, and PARD6A.

Discussion

In recent years, the survival rates and prognoses of BC patients have significantly improved, largely due to the clinical implementation of ICIs. Nevertheless, it is important to note that not all patients exhibit favorable responses to ICIs, with a considerable proportion demonstrating either resistance or insensitivity to these therapies. ²⁶ Research indicates that the TGF-β signaling pathway within the TME plays a critical role in immune evasion and resistance to ICIs in cancer. Consequently, the inhibition of TGF-β has been proposed as a strategy to enhance the efficacy of ICIs.^27,28 Furthermore, TGF-β is implicated in the promotion of BC metastasis through the induction of epithelial-mesenchymal transition (EMT), which is another significant factor contributing to poor prognoses in patients. ²⁹ In summary, the TGF-β signaling pathway is integral to the initiation, progression, and treatment of BC; however, its interactions with other elements within the TME and its implications for the prognosis of BC patients warrant further exploration.

In this study, genes associated with prognostic OS were identified through univariate Cox regression analysis utilizing the TCGA-BRCA dataset. These prognostic OS-related genes were subsequently cross-referenced with 225 TSRGs. Following this, LASSO regression analysis and multivariate Cox regression analysis were conducted on the intersecting genes, ultimately leading to the identification of an optimal prognostic model comprising five genes: FUT8, IFNG, ID3, KLF10, and PARD6A.

Fucosyltransferase 8 (FUT8) is frequently upregulated in various cancers due to its role in promoting core fucosylation, a process that is crucial for cancer cell metastasis and immune evasion. ³⁰ Through core fucosylation, FUT8 not only facilitates the dissemination of BC cells but also enhances TGF-β signaling and EMT processes. ³¹ Numerous studies have demonstrated that the inhibition of FUT8 can decelerate the progression of BC and improve patient outcomes.^32,33 Similarly, PARD6A is associated with EMT during the invasion of BC cells. ³⁴ The gene IFNG, which encodes interferon-gamma (IFN-γ), is vital for anti-tumor immunity, as it can be released by tumor-infiltrating immune cells (TIICs) to exert anti-tumor effects and can also induce the expression of PD-L1, thereby improving patient outcomes.^35,36 KLF10, a Krüppel-like zinc finger transcription factor, has been shown to play a significant role in inducing apoptosis via the Smad signaling pathway in response to TGF-β. ³⁷ Additionally, ID3 has been implicated in promoting the proliferation and invasion of human MCF-7 breast cancer cells. ³⁸ Notably, in our study, KLF10 was identified as a risk factor, while FUT8 and PARD6A were classified as protective factors. To further elucidate the underlying reasons for these findings, the expression levels of these five genes were examined in both normal and BC tissues. The results corroborate previous studies, indicating that FUT8 and PARD6A are more highly expressed in tumor tissues, whereas KLF10 exhibits greater expression in normal tissues. Therefore, it can be concluded that the observed results across different BC risk groups may be influenced by the molecular subtypes of BC, leading to variations in outcomes.

This study examined the potential association between TGF-β risk scores and the TME of BC patients. Initially, it was observed that patients with low TGF-β risk scores exhibited increased infiltration of gamma delta T cells (γδT cells) and T follicular helper cells (Tfh cells) within the TME. This finding was consistent when comparing normal and tumor tissues. γδT cells play a pivotal role in both innate and adaptive immunity, ³⁹ as they can recognize tumor cells independently of major histocompatibility complex (MHC) antigens and induce tumor cell death through potent anti-inflammatory and cytotoxic mechanisms. ⁴⁰ Tfh cells enhance anti-tumor immunity in a CD8+-dependent manner and are significant effector cells in anti-PD-1/PD-L1 therapies. ⁴¹ Subsequently, the study assessed the correlation between TGF-β risk scores and six immune checkpoints within the TME, revealing a negative correlation across all datasets analyzed. Consequently, patients classified in the low TGF-β risk group were more likely to experience favorable outcomes from ICI treatments, a finding corroborated by comparisons using TIDE and IPS between the two groups. Furthermore, the prevailing understanding of IFNG suggests that it can augment the efficacy of anti-PD-1 therapy. ⁴² The study found a positive correlation between IFNG expression and several immune checkpoints, including It is important to note that IFNG, predominantly produced by CD8+ T cells, is integral to T cell immunity and cancer immunotherapy. The IFNγ produced by CD8+ T cells downregulates the expression of SLC3A2 and SLC7A11, leading to decreased cystine uptake by tumor cells, which in turn induces lipid peroxidation and ferroptosis. ⁴³

The analysis simulated the differential treatment responses to various drugs between high- and low-risk groups, suggesting that low-risk patients exhibit greater sensitivity to ICIs. Additionally, it was predicted that targeted chemotherapy agents would be beneficial for both high- and low-risk groups to enhance survival in BC. Bortezomib, a proteasome inhibitor that impedes myeloma proliferation and exerts anabolic effects on bone, is a critical treatment for BC bone metastases. ⁴⁴ Although the relationship between IFNG and bortezomib is unclear, bortezomib has been used in combination with other medicines in patients with amyloid light-chain (AL), and transcriptome analysis has revealed CD8 memory T cell activation and IFNG overexpression. ⁴⁵ Everolimus, a derivative of rapamycin, inhibits BC cell proliferation and aggressiveness through the PI3K/AKT/mTOR signaling pathway. ⁴⁶ Specifically, IFNG released by natural killer cells is inhibited by mTOR inhibitors. ⁴⁷ To address tamoxifen resistance, the mTORC1/2 dual inhibitor (AZD8055) can downregulate HSPB8. ⁴⁸ BMS.754807, by blocking the insulin-like growth factor 1 receptor (IGF1R), inhibits BC cell proliferation and enhances their sensitivity to the chemotherapeutic agent cisplatin. ⁴⁹ Furthermore, sabutoclax effectively targets cancer stem cells in BC by inhibiting the IL-6/STAT3 signaling pathway. ⁵⁰ Currently, there is insufficient evidence to support the efficacy of other drugs, such as CDK9.1, BI.2536, and Dihydrorotenone, in the treatment of BC, necessitating further investigation into their relationship with the disease.

The present study acknowledges several limitations. Firstly, this study was based on public databases, and the expression levels of five TSRG genes were limited by the study's objective experimental conditions and were not clinically verified, despite the fact that FUT8, IFNG, ID3, and KLF10 were all confirmed in other clinical breast cancer sample studies. Secondly, due to variations in the data within the existing database, the receiver operating characteristic (ROC) curve analysis in the GSE42568 validation set was limited to evaluating the model's predictive capabilities at 1, 3, and 5 years, rather than at 3, 5, and 10 years. Thirdly, further research is warranted to refine the classification and typing criteria for BC.

Conclusion

In conclusion, this study presents a robust risk model for BC patients based on TSRGs. The model indicates that TGF-β risk scores can effectively characterize the TME in BC patients and align them with various sensitive therapeutic agents. Moreover, TGF-β has the potential to delineate molecular subtypes of BC, thereby paving the way for enhanced clinical patient classification and treatment strategies.

Supplemental Material

sj-docx-1-cix-10.1177_11769351251316398 – Supplemental material for Identification of Molecular Subtypes and Prognostic Features of Breast Cancer Based on TGF-β Signaling-related Genes

Supplemental material, sj-docx-1-cix-10.1177_11769351251316398 for Identification of Molecular Subtypes and Prognostic Features of Breast Cancer Based on TGF-β Signaling-related Genes by Jia Qu, Mei-Huan Wang, Yue-Hua Gao and Hua-Wei Zhang in Cancer Informatics