Identification of molecular subtypes for breast cancer based on butyrate metabolism-related genes to assess prognosis and immune landscape

Introduction

Breast cancer (BRCA) ranks among the most prevalent malignant tumors globally, accounting for roughly 11.5% of all new cancer diagnoses (1-3). Based on the information supplied by the International Agency for Research on Cancer of the World Health Organization, BRCA has emerged as the most common malignant tumor among women globally, presenting a significant threat to their health (4). In 2022, the latest global cancer statistics suggest that there were 2.32 million new cases of BRCA (5). Drug resistance, recurrence, and metastasis have been identified as the three primary factors influencing the survival prognosis of individuals with BRCA (6). BRCA is marked by elevated rates of incidence, mortality, and recurrence (7). In the last twenty years, the treatment strategies for BRCA have advanced, integrating a range of methods including surgery, chemotherapy, radiation, targeted therapies, immunotherapy, and other innovative technologies (8-10). Neoadjuvant chemotherapy has shown potential in tailoring BRCA treatments, though its efficacy can differ between individuals and may result in the development of drug resistance (11-13). With immune therapies advance, they have demonstrated encouraging outcomes in the treatment of BRCA. However, a more thorough understanding of personalized immunotherapy selection for different BRCA subtypes still remains to deeper understand.

Butyrate is structurally classified as a short-chain fatty acid, a primary byproduct of gut microbial fermentation, and it is considered a key regulator of gut microbiota in maintaining overall energy balance in the body (14,15). Butyrate has demonstrated protective effects on immune function and intestinal disorders, including graft-versus-host disease, inflammatory bowel disease in the gastrointestinal tract, colon cancer, and these benefits are particularly valuable for conditions like ulcerative colitis and other ailments that impact the intestinal lining (16,17). Commonly, butyrate is crucial in the fight against cancer and anti-inflammation properties (18).

In 2021, researchers have found that butyrate generated by Escherichia coli from the healthy human gut microbiota and has been shown to suppress cancer cell growth (19). Furthermore, research has indicated that concentrations of 2–5 mM butyric acid can decrease the viability of BRCA cells (7). At present, comprehensive research remains limited focusing on genes involved in butyrate metabolism in BRCA.

In this study, we developed a predictive model by applying univariate Cox and least absolute shrinkage and selection operator (LASSO) regression analysis to genes involved in butyrate metabolism, which was validated using the BRCA dataset. We further examined the variations in immune cell infiltration, gene mutations, chemotherapeutic drug response, and immunotherapy effectiveness in high- and low-risk BRCA patients. We also identified eight key genes associated with butyrate metabolism, which could offer new potential markers for treating BRCA patients. Additionally, we conducted a detailed comparison of these eight genes, examining their relationship with survival outcomes, clinical features, and immune cell profiles. The prognostic model developed may act as an essential instrument for prognosis prediction and guide clinical treatment decisions for BRCA patients. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1035/rc).

Methods

Data collection

By data mining The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/), we downloaded the raw RNA sequencing (RNA-seq) data and collected relevant data (mutation and clinical data) for 1,113 BRCA and 113 standard breast tissue samples. Notably, samples with missing expression data, incomplete clinical information, or very short overall survival times (with survival time equal to zero) were excluded from this analysis. A total of 1,038 BRCA samples were examined in the research. Using the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), we retrieved the independent validation dataset comprising 327 BRCA samples (GSE20685). In addition, genes associated with butyrate metabolism were sourced from existing literature (20). This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Constructing butyrate metabolism-related genes subtypes

To construct a prognostic model using differentially expressed genes (DEGs) linked to butyrate metabolism, univariate Cox regression analysis was utilized, establishing a P value threshold of 0.05. The BRCA dataset was then clustered (K=2) into two subtypes using the K-means method. DEGs between the two subtypes were analyzed employing the “limma” package, followed by analyses of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways to further explore butyrate metabolism-related genes enrichment related signaling pathway.

Screening hub gene prognosis-related trait modeling

Through univariate Cox regression, a prognostic model was developed and derived from 3,071 DEGs (up =34, down =3,037) that differed between subtypes, setting the screening criteria, |log fold change (FC)| >0.585, P<0.05. To mitigate overfitting of the model, LASSO analysis was conducted on these selected candidate genes employing the “glmnet” package, using the penalty parameter (lambda) to prune genes with high correlation, thereby simplifying the model. In addition, candidate genes were screened using the “randomForestSRC” package. The “survival” package was utilized to develop a prognostic model by conducting univariate Cox regression analysis on the selected candidate genes. Risk scores for each sample were computed with the formula provided below:

Riskscore=∑Coef(i)∗Exp(i)[1]

According to the median risk score, all BRCA patients were classified into categories of low- and high-risk groups. Kaplan-Meier survival curves were generated employing the “survival” package in R. To assess the predictive accuracy of the model, a receiver operating characteristic (ROC) curve was plotted with the “timeROC” package and area under the curve (AUC) values computed for ROC curves at 1-, 3-, and 5-year. Further, visualizations were created to depict the score distribution, survival status, and expression levels were visualized as heatmaps for both high- and low-risk groups.

Independent prognostic analysis with the creation of a nomogram

We examined whether the BRCA risk score could serve as an independent factor for predicting of overall survival through both univariate and multivariate prognostic analyses. To create an intuitive predictive model relying on independent clinical parameters for assessing the overall survival of BRCA patients, we utilized the “rms” package in R. The association between risk scores and clinical factors was investigated using a nomogram, which displayed risk scores alongside variables like age, gender, disease stage, and tumor-node-metastasis (TNM) stage. Nomogram calculated the chances of 1- , 3- , and 5-year survival rates for patients with BRCA. Calibration plots and decision curve analysis (DCA) were utilized for assessing the predictive performance for this model.

Immune cell infiltration and immune checkpoint analysis

We assessed the types and functional immune responses of immune cells by calculating immune infiltration scores through the gene set variation analysis (GSVA). Subsequently, we determined the immune score, stromal score, ESTIMATE score, and tumor purity employing the “Estimate” package for analysis. Further, we performed immune infiltration analysis through single-sample gene set enrichment analysis (ssGSEA). The box plot was created to illustrate immune cell infiltration levels across high- and low-risk groups by using CIBERSORT. Moreover, we measured the expression profiles of immune checkpoints, and BRCA scoring data were sourced from The Cancer Immunome Atlas (TCIA; https://tcia.at), facilitating the classification of cancer tissue samples into six immune subtypes: C1 (wound healing), C2 [interferon-γ (IFN-γ) dominant], C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 [transforming growth factor-β (TGF-β) dominant-characterized].

Enrichment analysis for high- and low-risk groups

With the use of the prognostic scoring model, the data sets were categorized into high- and low-risk, followed by pathway enrichment analysis in each group through GSEA. By utilizing the “limma” R package, the DEGs were conducted across the high- and low-risk groups, with selection criteria of |logFC| >0.585 and P<0.05. Then, GO and KEGG analyses were carried out to examine the signaling pathways and biological processes related to the DEGs across high- and low-risk groups, with the “clusterProfiler” package.

Tumor mutation burden analysis

Using mutation data from TCGA-BRCA, each patient’s tumor mutational burden (TMB) was counted and contrasted in the high- and low-risk groups. Continuously, we contrasted survival outcomes based on different TMB scores. Collecting the mutation information of top 20 in high- and low-risk groups, waterfall plots were drawn employing the “GenVisR” package.

Drug sensitivity prediction

For the purpose of discovering novel promising therapeutic targets and more effective anti-cancer drugs, we conducted a screening of drugs targeting butyrate metabolism-related genes and their correlation with drug sensitivity and prognosis, utilizing the CellMiner database (https://discover.nci.nih.gov/cellminer/). Moreover, we employed the “pRRophetic” package for estimating the half-maximal inhibitory concentration (IC₅₀) values for a range of pharmaceutical compounds.

Cell lines and cell culture

Human mammary epithelial cells (MCF-10A) and BRCA cells (MDA-MB-231) were acquired from Pricella Biotechnology Co., Ltd. (Wuhan, China). The cells were maintained separately in DMEM/F12 and Leibovitz’s L-15 media (both from Pricella), each supplemented with 10% fetal bovine serum, and incubated at 37 ℃ with 5% CO₂.

Quantitative real-time polymerase chain reaction (qRT-PCR)

Total RNA was isolated from the cells, and cDNA synthesis was performed using the QuantiText Reverse Transcription Kit (QIAGEN, Hilden, Germany) according to the manufacturer’s protocol. qRT-PCR was subsequently carried out using Hieff UNICON Universal Blue qPCR SYBR Green Master Mix (Yeasen, Shanghai, China). The amplification protocol consisted of an initial denaturation at 95 ℃ for 30 seconds, followed by 40 cycles at 95 ℃ for 3 seconds and 60 ℃ for 20 seconds. Relative messenger RNA (mRNA) expression levels were quantified by the 2^−ΔΔCt method. The sequences of all primers are shown in Table 1.

Statistical analysis

Statistical analyses were performed using R software and GraphPad Prism 9. Univariate and multivariate Cox regression analyses were conducted with the “survival” package, reporting hazard ratios and 95% confidence intervals, with significance at P<0.05. Kaplan-Meier survival curves were compared using the log-rank test. ROC curve AUC values were calculated with the “timeROC” package, using DeLong’s test for significance. Differential expression analyses (“limma” package) used |logFC| >0.585 and P<0.05, with multiple testing correction via the Benjamini-Hochberg false discovery rate (FDR) method unless otherwise specified. Wilcoxon rank-sum tests or t-tests were used for group comparisons (immune infiltration, TMB, drug sensitivity), with FDR correction for multiple comparisons where applicable. qRT-PCR data were analyzed with paired t-tests (n=3 replicates), and P<0.05 was considered significant.

Results

Construction of butyrate metabolism-related molecular subtypes based on prognosis-related genes

In our study, we selected six butyrate metabolism-associated genes associated with prognosis with P value <0.05 through univariate Cox regression analysis (Figure 1A). The cumulative distribution function (CDF) curve was generated to evaluate the area under the distribution curve, which revealed that the change in area proportion was optimal when two subtypes were assumed (Figure 1B,1C). Principal component analysis (PCA) analysis reliably and consistently confirmed these distinctions (Figure 1D). To explore the expression levels of six butyrate metabolism-associated genes between the two subtypes, the comprehensive heatmap and box plot were performed, which illustrated all of genes were upregulated in cluster 1 compared with cluster 2 (Figure 1E,1F). Furthermore, we also found the six butyrate metabolism-associated genes clearly divided the normal and tumor samples into two directions in space by PCA analysis (Figure 1G), and most of genes were overexpressed within tumor samples related to non-tumor samples (Figure 1H).

Figure 1 Construction of butyrate metabolism-related molecular subtypes. (A) A forest plot depicting the butyrate metabolism-associated genes that are strongly linked to the prognosis of BRCA. (B) An observed rise within the area beneath the CDF curve as the assumed quantity of molecular subtypes increases. (C) The consensus clustering matrix reveals two distinct molecular subtypes. (D) PCA analysis on the two molecular subtypes. (E) A heatmap illustrating the expression profiles of genes linked to butyrate metabolism and prognosis across the two subtypes. (F) A box plot comparing the expression profiles of butyrate metabolism-related genes between the two subtypes. (G) PCA analysis on the two groups of BRCA and normal samples. (H) Box plot of differences in expression profiles of butyrate metabolism-related genes in tumorous and non-tumor tissues. ns, not significant (P>0.05); *, P<0.05; ****, P<0.0001. BRCA, breast cancer; CDF, cumulative distribution function; CI, confidence interval; Exp, expression; HR, hazard ratio; PCA, principal component analysis.

Identification of prognostic biomarkers based on BRCA molecular subtypes

To verify the survival outcomes between the two subtypes, the Kaplan-Meier method was performed, and cluster 2 was found to be better than cluster 1 (Figure 2A). Besides, we obtained 3,071 DEGs between cluster 1 (n=744) and cluster 2 (n=294) (Figure 2B). GO analysis showed that butyrate metabolism-related DEGs of upregulated and downregulated were mainly associated with biological process (including regulation of membrane potential, striated muscle tissue development, and muscle organ development) and cell component (including postsynaptic specialization, monoatomic ion channel complex, and asymmetric synapse) (Figure 2C,2D). KEGG analysis showed that DEGs of upregulated and downregulated are related to neuroactive ligand-receptor interaction and cytoskeleton in muscle cells pathways (Figure 2E,2F). Furthermore, to eliminate co-expressed butyrate metabolism-related prognostic genes and prevent overfitting, a total of 39 genes were selected for the establishment of the prognostic model through LASSO analysis (Figure 2G). Subsequently, we selected 17 feature genes using random forest and take the intersection of between 39 genes and 17 feature genes to get 12 signature genes (Figure 2H,2I). Additionally, a multivariate regression analysis was performed to enhance the gene identification and compute the regression coefficients (Figure 2J):

Riskscore=−0.1477∗CD79A+0.0527∗CLEC3A−0.0477∗CLIC6+0.3885∗DCTPP1−0.2392∗LDLRAD3−0.2774∗RBBP8+0.3524∗STXBP5+0.3444∗ZDHHC9[2]

Figure 2 Butyrate metabolism-associated genes that have substantial prognostic significance in BRCA. (A) The survival analysis of patients with two molecular subtypes of BRCA. (B) A volcano plot displaying DEGs across various molecular subtypes of BRCA. (C,D) GO analysis of DEGs that are (C) up-regulated and (D) down-regulated. (E,F) KEGG enrichment analysis for the (E) up-regulated and (F) down-regulated DEGs. (G) LASSO coefficient distribution of 39 butyrate metabolism-related genes. (H) The RF model of the number of trees and error, the lowest error rate occurred when the number of trees was 400, and the RF model of variable importance and 39 signature genes. (I) A Venn diagram highlighted 12 genes by intersecting features from the LASSO and RF algorithms. (J) Forest plot for multivariate regression analysis regarding OS in BRCA patients. BRCA, breast cancer; CI, confidence interval; DEGs, differentially expressed genes; GABA, gamma-aminobutyric acid; GO, Gene Ontology; HR, hazard ratio; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; OS, overall survival; RF, random forest.

Validation of butyrate metabolism-related prognostic models

To validate the robustness and generalizability of our model, we employed the TCGA-BRCA cohort for internal training and the GSE20685 cohort for external validation. ROC curves were generated, showing the AUC for BRCA patients: 0.758, 0.786, and 0.744 for 1-, 3-, and 5-year outcomes in the training cohort, and 0.79, 0.71, and 0.683 for 1-, 3-, and 5-year outcomes in the GSE20685 cohort (Figure 3A,3B). Further, both the training and GSE20685 cohorts showed that the patients in the low-risk group exhibited superior clinical prognosis than those in the high-risk group (Figure 3C,3D). Risk scores for each sample were computed independently for the training cohort and the GSE20685 validation cohort using the same risk formula. It was observed that higher risk scores in both cohorts were correlated with poorer OS and an increased risk of mortality among BRCA patients (Figure 3E,3F).

Figure 3 Butyrate metabolism-associated genes that exhibit significant prognostic importance in BRCA. (A,B) ROC curve for predicting 1-, 3-, and 5-year survival rates in the (A) training set and (B) GSE20685 set. (C,D) Kaplan-Meier survival curves showing the differences in survival rates between the high- and low-risk groups in the (D) training set and (D) GSE20685 set. (E,F) Risk score distribution and survival outcomes for BRCA patients in the (E) training set and (F) GSE20685 set. AUC, area under the curve; BRCA, breast cancer; ROC, receiver operating characteristic.

Development of nomogram for independent prognostic assessment in BRCA patients

To assess the additional influence of signature genes on the prognosis of BRCA patients, survival curves based on signature genes were performed, revealing that reduced expression of four genes was related to notably improved OS compared to higher expression levels, including DCTPP1, STXBP5, CLEC3A, and ZDHHC9 (Figure 4A-4G). Besides, the risk scores of the two molecular subtypes were compared, and our analysis revealed that cluster 1 exhibited a higher risk score compared to cluster 2 (Figure 4H). In order to determine whether risk score could function as an independent predictor of OS in BRCA patients, both univariate and multivariate Cox regression analyses were carried out. The findings showed that risk score and age are independent factors for forecasting the overall survival of BRCA patients (Figure 4I,4J). The DCA curves presented that the prognostic model exhibited excellent net benefit for 1-, 3-, and 5-year (Figure 4K). A nomogram integrating risk score, grade, TNM stage, age, and gender was developed to estimate the 1-, 3-, and 5-year survival rates for BRCA patients (Figure 4L). Further, a calibration curve for 1-, 3-, and 5-year survival outcomes suggested that the nomogram provided accurate predictions for OS of BRCA patients (Figure 4M). Additionally, we compared the risk scores between the high- and low-risk groups across four intrinsic subtypes: basal, luminal A (LumA), luminal B (LumB), and human epidermal growth factor receptor 2 (HER2), showing that the high-risk group exhibited significantly higher risk scores than the low-risk group in all four subtypes (Figure 4N).

Figure 4 Development of nomogram. (A-G) Kaplan-Meier curves of feature genes (A) DCTPP1, (B) STXBP5, (C) CD79A, (D) CLEC3A, (E) CLIC6, (F) LDLRAD3, and (G) ZDHHC9. (H) The box plot of risk score in cluster 1 and cluster 2. (I) Univariate and (J) Multivariate analysis was conducted to assess whether the nomogram could serve as an independent prognostic factor. (K) Analysis of decision curve of nomogram. (L) Nomogram to predict the OS of 1-, 3-, and 5-year. (M) Calibration curves to assess the precision of the nomogram. (N) Box plots compare the risk scores between high- and low-risk patients within each subtype: basal, LumA, LumB, and HER2. ****, P<0.0001. HER2, human epidermal growth factor receptor 2; LumA, luminal A; LumB, luminal B; M, metastasis; N, node; OS, overall survival; T, tumor.

Immune infiltration analysis of butyrate metabolism-related signature genes and prediction of efficacy

To assess what role immune cells play in the progression of BRCA, we employed the CIBERSORT algorithm to examine the variations in immune cell populations across the high- and low-risk groups (Figure 5A). The box plot revealed that nine immune cell types were elevated within the low-risk group, such as follicular helper T cells, activated CD4 memory T cells, M1 macrophages, activated natural killer (NK) cells, plasma cells, resting CD4 memory T cells, resting dendritic cells, naive B cells, and CD8 T cells (Figure 5B). Subsequently, using ESTIMATE, we observed markedly increased ESTIMATE, immune, and stromal scores, along with reduced tumor purity in the low-risk group (Figure 5C). To gain a deeper understanding of the impact of immune status on the prognosis of BRCA, ssGSEA was employed to assess immune cell infiltration and immune functions in high- and low-risk patients (Figure 5D,5E). The results revealed that all of the immune cells exhibited higher expression in the low-risk group, such as B cells, T cells, CD8⁺ T cells, and so on, and the immune function analysis revealed greater activity in the low-risk group relative to the high-risk group. Additional analysis of immune checkpoint gene expression in the two groups showed that most immune-related genes were expressed at higher levels in the low-risk group (Figure 5F). A supplementary analysis was carried out to evaluate the immunophenoscore (IPS) scores of BRCA patients across both the high- and low-risk groups (Figure 5G). Our findings demonstrated that the IPS scores were notably elevated in the low-risk group when contrasted with the high-risk group. This suggested that BRCA patients categorized as low-risk could potentially benefit from immunotherapy. It was observed that the high- and low-risk groups exhibited significant correlations with different subtypes, including wound healing (C1), IFN-γ-dominant (C2), and inflammatory (C3) (Figure 5H).

Figure 5 Immunological infiltration across risk groups. (A) Accumulation plots showing the distribution of immune cells in two groups. (B) Box plot displaying the distribution of tumor-associated immune cells in two groups. (C) The violin plots were created to visualize the estimate score, immune score, stromal score, and tumor purity. (D,E) Stack and box plots were employed to illustrate immune cell infiltration and to analyze immune functions in two groups. (F) The box plot graph illustrated the expression levels of immune checkpoint-related genes across two groups. (G) The IPS scores across the high- and low-risk groups. (H) A Sankey diagram depicting the connections between five immune subtypes, high- and low-risk groups, and the two clusters from the TCGA dataset. ns, not significant (P>0.05); *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. aDCs, artificial dendritic cells; APC, antigen-presenting cell; CCR, C-C chemokine receptor; DCs, dendritic cells; iDCs, interdigitating dendritic cells; IFN, interferon; IPS, immunophenoscore; MHC, major histocompatibility complex; NK, natural killer; pDCs, plasmacytoid dendritic cells; TCGA, The Cancer Genome Atlas; Tfh, follicular helper T; Th, T helper; TIL, tumor-infiltrating lymphocyte; Treg, T regulatory.

Enrichment analysis across risk groups

To evaluate the signaling pathways associated with butyrate metabolism genes more thoroughly, GSEA was applied to assess the enrichment of pathways in both groups (Figure 6A,6B). The findings revealed by the GSEA analysis revealed remarkable differences in the enrichment of gene sets within the high- and low-risk groups. Upregulated expression of butyrate metabolism-related genes was linked to processes such as the cell cycle, DNA replication, mismatch repair, and backbone biosynthesis. Conversely, reduced expression of these genes was associated with signaling pathways, receptor interactions, cell lineage, and primary immunodeficiency. The results of GO analysis showed that the upregulated DEGs were connected with cellular components such as the secretory granule lumen, cytoplasmic vesicle lumen, and vesicle lumen (Figure 6C), and the downregulated DEGs were related to lymphocyte differentiation, lymphocyte-mediated immunity, and immune receptors built from immunoglobulin superfamily in biological process (Figure 6D). KEGG pathway analysis highlighted that the upregulated DEGs were predominantly linked to the interleukin-17 (IL-17) signaling pathway and bladder cancer, while the downregulated DEGs were primarily connected to hematopoietic cell lineage, cytokine-cytokine receptor interactions, and primary immunodeficiency (Figure 6E,6F).

Figure 6 Function enrichment of GSEA analysis in both risk groups. (A) GSEA findings for the high-risk group. (B) GSEA findings for the low-risk group. (C) GO enrichment of up-regulated genes in both risk groups. (D) GO enrichment of down-regulated genes in both risk groups. (E) Pathway enrichment of up-regulated genes comparing the two groups in the KEGG database. (F) Pathway enrichment of down-regulated genes comparing the two groups in the KEGG database. BP, biological process; CC, cellular component; GO, Gene Ontology; GSEA, gene set enrichment analysis; IgM, immunoglobulin M; IL-17, interleukin-17; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; NF, nuclear factor.

Tumor mutation burden analysis

To examine the variations in mutation status and expression patterns between the two risk groups, TMB analysis was conducted. The results revealed that there was a notable difference across the high- and low-risk groups, with the high-risk group exhibiting a significantly higher TMB (Figure 7A). Further, we visualized and analyzed mutation data by focusing on the top 20 genes associated with both high- and low-risk groups (Figure 7B,7C). Within the high-risk group, the three genes with the most frequent mutations were TP53 (41%), PIK3CA (28%), and TTN (19%). In contrast, the low-risk group showed the highest mutation frequencies in PIK3CA (39%), TP53 (26%), and CDH1 (15%).

Figure 7 TMB analysis in high- and low-risk group. (A) Violin plot revealing the distinction between the two risk groups in TMB. (B) Mutation analysis of high-risk group. (C) Analysis of mutations in the low-risk group. ****, P<0.0001. TMB, tumor mutational burden.

Drug sensitivity analysis in patients with BRCA

In treating BRCA patients, potential drugs could serve as targeted therapies. Therefore, we further explored the difference in sensitivity to anticancer drugs across the two risk groups. Drug screening relying on gene expression levels in the prognostic model, using the CellMiner database, identified a significant negative correlation between the IC₅₀ values for BP-1-102 and the expression of CEBPD. Similarly, the expression of CD79A showed a significantly reverse correlation with IC₅₀ on pluripotin. In contrast, the expression of EGFL7 was strongly associated with the IC₅₀ of karenitecin and gemcitabine (Figure 8A). Further, we investigated the relationship between risk scores and IC₅₀ values of cyclopamine and docetaxel, showing that the high-risk group had a notably higher IC₅₀ values for cyclopamine and docetaxel (Figure 8B). These results hold substantial significance for the drug selection in the clinical management of BRCA.

Figure 8 Drug sensitivity analysis. (A) Relationship between signature gene expression levels and sensitivity to various drugs provided by the CellMiner database. (B) Drug sensitivity analysis on cyclopamine and docetaxel across the two risk groups. **, P<0.01. IC₅₀, half-maximal inhibitory concentration.

Validation of expression of characteristic genes that comprised the risk model by qRT-PCR

To validate the expression of the identified signature genes in vitro, qRT-PCR analysis was performed using human mammary epithelial cells (MCF-10A) and BRCA cells (MDA-MB-231). It was observed that CLEC3A, DCTPP1, LDLRAD3, and ZDHHC9 were significantly overexpressed in tumor samples relative to normal samples, while CLIC6, RBBP8, and STXBP5 were significantly downregulated (Figure 9). Therefore, we proposed that altered expression patterns of the identified genes could impact the malignant behavior of BRCA.

Figure 9 Validation of characteristic gene expression by qRT-PCR. qRT-PCR analysis of CLEC3A, CLIC6, DCTPP1, LDLRAD3, RBBP8, STXBP5, and ZDHHC9 in MCF-10A and MDA-MB-231 cells, showing relative mRNA expression. *, P<0.05. mRNA, messenger RNA; qRT-PCR, quantitative real-time polymerase chain reaction.

Discussion

BRCA, with its distinct epidemiological trends and considerable heterogeneity, is consistently one of the primary contributors to cancer-related mortality among women (21). The Chinese Women’s survey indicates that BRCA is the most prevalent malignancy among Chinese women, with its incidence rising annually (22). BRCA has a complex pathogenesis and a wide range of clinical presentations, creating major challenges for its treatment and prevention (23,24). Conventional treatment options for BRCA, such as surgery, chemotherapy, radiation therapy, endocrine treatment, targeted therapies, and other related methods, still have many problems and need a more comprehensive understanding of its molecular foundations and associated risk factors because of the complex etiology (25-27). Consequently, it is essential to develop accurate prognostic models for BRCA to predict patient outcomes and treatment responses. Butyrate is a prevalent short-chain fatty acid synthesized by the gut microbiota and has the ability to modulate epigenetic machinery by inhibition of histone deacetylase, which has been shown to influence the progression of BRCA (28). For example, propionate and butyrate have been found to suppress the growth of the MCF7 BRCA cell line and induce cell cycle arrest at the G1 phase (29). Besides, the sub-therapeutic doses of butyrate can impact the BRCA cell cultures. In 2024, another study revealed that butyrate is capable of activating hormone receptors, increasing the transcription of estrogen-responsive genes, and facilitating the migration of BRCA cells (30). Nevertheless, research on the role and regulatory mechanisms of butyrate in BRCA development remains sparse. This study seeks to investigate the potential connections across butyrate metabolism-related genes and BRCA, set the prognosis model, and provide guidance for the BRCA patients.

In this study, we classified two separate subtypes of BRCA by integrating butyrate metabolism-related genes from the GEO and TCGA databases. Further, we identified eight key signature genes to build a novel prognostic model through LASSO and univariate Cox regression analyses. The butyrate metabolism-related signature genes we identified were found to serve as independent prognostic indicators for BRCA and were split into two different prognostic groups using median risk scores as a reference. We then constructed ROC curves, a nomogram, and calibration curves. A thorough analysis revealed that the BRCA signature exhibited superior predictive accuracy in contrast to conventional clinical factors like age, gender, histological stage, and tumor grade. Meanwhile, the predicted values closely align with the observed data, offering a solid theoretical foundation to support clinicians’ decision-making.

The characteristic genes we identified included CD79A, CLEC3A, CLIC6, DCTPP1, LDLRAD3, RBBP8, STXBP5, and ZDHHC9. CD79A, a protein linked to the B cell receptor, is crucial for the development and function of B cells (31). CD79A is a key target in classical Hodgkin’s lymphoma as well (32). CLEC3A, a peptide originating from the human cartilage-specific C-type lectin domain family 3 member A, exhibited strong antimicrobial activity (33). A previous study has indicated that the reduced expression of CLEC3A was related to favorable survival outcomes in BRCA (34). A separate study also found that CLEC3A could function as an independent prognostic factor for the survival of BRCA individuals (35). CLIC6 plays an essential role in cancer development and has been extensively investigated as an important therapeutic target for various types of cancer. In 2025, Zhou et al. have revealed that CLIC6 exhibits strong antitumor properties in hepatocellular carcinoma by suppressing cell proliferation, enhancing apoptosis, regulating cytokine levels, modulating immune cell dynamics, and reducing oxidative stress pathways (36). Further, the immunohistochemistry results indicated that CLIC6 was expressed at significantly lower levels in BRCA tissues, and has been identified as a prognostic biomarker for BRAC patients (37). DCTPP1 is a key hydrolase involved in breaking down noncanonical dCTP and ensuring the balance of the deoxyribonucleoside triphosphate (dNTP) pool, making it a potential target for cancer therapies related to nucleotide metabolism (38). Several studies have indicated that DCTPP1 may serve as a promising biomarker and therapeutic target for a range of cancers, including colorectal, ovarian, and BRCA (39-41). LDLRAD3, a protein containing a class A domain of the low-density lipoprotein receptor, is widely recognized as the receptor for the venezuelan equine encephalitis virus (42). Li et al. proposed that LDLRAD3 promoted oncogenesis in non-small-cell lung cancer by regulating the miR-20a-5p-SLC7A5 pathway, which in turn activated the mTORC1 signaling cascade, which presents a potential therapeutic target for NSCLC treatment and management (43). RBBP8 has limited research focused on its role in BRCA. It has been suggested that hypomorphic mutations, rather than complete loss-of-function mutations, in RBBP8 are related to an increased risk of early-onset BRCA (44). STXBP5, also known as syntaxin binding protein 5, is frequently recognized as a key regulator of progerin expression, with its overexpression contributing to the early onset of cellular senescence, which associated with various cancers (45). Wang et al. found that the expression levels of mRNA and proteins for STXBP5 were notably reduced in glioma tissues, and these levels were positively associated with patient prognosis (46). Regarding BRCA, Park et al. demonstrated that ginsenoside Rh2 regulates the STXBP5-AS1/miR-4425/RNF217 pathway, which in turn inhibits the proliferation of BRCA cells (47). ZDHHC9 has been recognized as a crucial determinant of resistance to immune checkpoint blockade therapy, serving as an important marker for the classification, diagnosis, and prognosis of triple-negative BRCA, as well as for the development of targeted therapeutic strategies (48). Overall, these genes may function as biomarkers for patient’s prognosis with BRCA.

The role of tumor microenvironments in shaping various tumor phenotypes is now widely accepted. Our newly developed risk prediction model offers valuable insights for forecasting patient treatment outcomes and exploring the connections between target genes, immune cell infiltration, and clinical features. To investigate the role of immune cells in the development and progression of BRCA, we carried out an analysis of immune-related functionalities and infiltration scores of various immune cell types. Analysis of immune cell infiltration showed that a higher abundance of artificial dendritic cells (aDCs), B cells, CD8⁺ T cells, T cells, NK cells, mast cells, neutrophils, T helper (Th)1, Th2, T regulatory (Treg), and macrophages in the low-risk group. The approach of immunotherapy aims to hinder cancer progression by stimulating the innate immune molecules within the tumor microenvironment (49). A study has identified CD8⁺ cytotoxic T cells, Th1 helper cells, and their related cytokines, like interferon, as key components with antitumor activity within the tumor microenvironment (50). The low-risk group exhibits a higher prevalence of Th1 cells, CD8⁺ T cells, and effector memory T cells, which is associated with a more favorable prognosis (51). Treg cells secrete various cytokines that boost their capacity to combat tumors (52). Further, IPS analysis was performed between high- and low-risk groups, revealing that the low-risk group that had a higher IPS score implied an elevated level of immunological reactivity in contrast to the high-risk group. As a result, the evidence suggests that immunotherapy could enhance the treatment outcomes for BRCA patients in the low-risk group, potentially leading to outstanding efficacy.

To conclude, we have built a butyrate metabolism-related genes prognostic model, which helps to predict drug sensitivity and survival outcomes in BRCA patients. While our research offers significant clinical insights for prognostic evaluation and treatment selection in BRCA patients, it is not without its limitations. Initially, our research is based on retrospective data, and further validation through prospective studies is required. In addition, further biological studies are required to investigate the impact of butyrate metabolism-associated genes on immune function.

Conclusions

Overall, we created an eight-gene signature linked to butyrate, which can effectively forecast the prognosis of BRCA patients. This research enhances our comprehension of the function of butyrate metabolism-related genes in BRCA progression and provides valuable insights for guiding immunotherapy in BRCA patients.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1035/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1035/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1035/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1035/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Zhang X, Yang L, Liu S, et al. Interpretation on the report of global cancer statistics 2022. Zhonghua Zhong Liu Za Zhi 2024;46:710-21. [Crossref] [PubMed]
2. Qiu H, Cao S, Xu R. Cancer incidence, mortality, and burden in China: a time-trend analysis and comparison with the United States and United Kingdom based on the global epidemiological data released in 2020. Cancer Commun (Lond) 2021;41:1037-48. [Crossref] [PubMed]
3. Manzanedo I, Pereira F, Pérez-Viejo E, et al. Gastric Cancer with Peritoneal Metastases: Current Status and Prospects for Treatment. Cancers (Basel) 2023;15:1777. [Crossref] [PubMed]
4. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
5. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
6. Tufail M, Cui J, Wu C. Breast cancer: molecular mechanisms of underlying resistance and therapeutic approaches. Am J Cancer Res 2022;12:2920-49.
7. Zhang L, Huang S, Yuan Y. Butyrate inhibits the malignant biological behaviors of breast cancer cells by facilitating cuproptosis-associated gene expression. J Cancer Res Clin Oncol 2024;150:287. [Crossref] [PubMed]
8. Lau KH, Tan AM, Shi Y. New and Emerging Targeted Therapies for Advanced Breast Cancer. Int J Mol Sci 2022;23:2288. [Crossref] [PubMed]
9. Kerr AJ, Dodwell D, McGale P, et al. Adjuvant and neoadjuvant breast cancer treatments: A systematic review of their effects on mortality. Cancer Treat Rev 2022;105:102375. [Crossref] [PubMed]
10. Deng B, Kong W, Suo H, et al. Oleic Acid Exhibits Anti-Proliferative and Anti-Invasive Activities via the PTEN/AKT/mTOR Pathway in Endometrial Cancer. Cancers (Basel) 2023;15:5407. [Crossref] [PubMed]
11. Loi S, Salgado R, Curigliano G, et al. Neoadjuvant nivolumab and chemotherapy in early estrogen receptor-positive breast cancer: a randomized phase 3 trial. Nat Med 2025;31:433-41. [Crossref] [PubMed]
12. Villacampa G, Navarro V, Matikas A, et al. Neoadjuvant Immune Checkpoint Inhibitors Plus Chemotherapy in Early Breast Cancer: A Systematic Review and Meta-Analysis. JAMA Oncol 2024;10:1331-41. [Crossref] [PubMed]
13. Wu Y, Huang S, Wei Y, et al. Efficacy and safety of different regimens of neoadjuvant therapy in patients with hormone receptor-positive, her2-negative breast cancer: a network meta-analysis. Front Immunol 2024;15:1420214. [Crossref] [PubMed]
14. Sun J, Chen S, Zang D, et al. Butyrate as a promising therapeutic target in cancer: From pathogenesis to clinic Int J Oncol 2024;64:44. (Review). [Crossref] [PubMed]
15. Gu R, Mei K, Chen Z, et al. Development and validation of a prognosis prediction model for overall survival in correlation between butyrate metabolism and gastric cancer prognosis: Mendelian randomization and transcriptomics analysis. Transl Cancer Res 2025;14:743-60. [Crossref] [PubMed]
16. Cui X, Li C, Zhong J, et al. Gut microbiota - bidirectional modulator: role in inflammatory bowel disease and colorectal cancer. Front Immunol 2025;16:1523584. [Crossref] [PubMed]
17. Hodgkinson K, El Abbar F, Dobranowski P, et al. Butyrate's role in human health and the current progress towards its clinical application to treat gastrointestinal disease. Clin Nutr 2023;42:61-75. [Crossref] [PubMed]
18. Kang X, Liu C, Ding Y, et al. Roseburia intestinalis generated butyrate boosts anti-PD-1 efficacy in colorectal cancer by activating cytotoxic CD8(+) T cells. Gut 2023;72:2112-22. [Crossref] [PubMed]
19. Guo H. Interactions between the tumor microbiota and breast cancer. Front Cell Infect Microbiol 2024;14:1499203. [Crossref] [PubMed]
20. Yuan P, Li S, Xiong S, et al. Development of a butyrate metabolism-related gene-based molecular subtypes and scoring system for predicting prognosis and immunotherapy response in bladder cancer. J Cancer Res Clin Oncol 2023;149:12489-505. [Crossref] [PubMed]
21. Xiong X, Zheng LW, Ding Y, et al. Breast cancer: pathogenesis and treatments. Signal Transduct Target Ther 2025;10:49. [Crossref] [PubMed]
22. Sun K, Zhang B, Lei S, et al. Incidence, mortality, and disability-adjusted life years of female breast cancer in China, 2022. Chin Med J (Engl) 2024;137:2429-36. [Crossref] [PubMed]
23. Guo L, Kong D, Liu J, et al. Correction: Breast cancer heterogeneity and its implication in personalized precision therapy. Exp Hematol Oncol 2024;13:7. [Crossref] [PubMed]
24. Ishwar D, Premachandran S, Das S, et al. Profiling Breast Tumor Heterogeneity and Identifying Breast Cancer Subtypes Through Tumor-Associated Immune Cell Signatures and Immuno Nano Sensors. Small 2024;20:e2406475. [Crossref] [PubMed]
25. Hickey M, Basu P, Sassarini J, et al. Managing menopause after cancer. Lancet 2024;403:984-96. [Crossref] [PubMed]
26. Hu C, Hart SN, Gnanaolivu R, et al. A Population-Based Study of Genes Previously Implicated in Breast Cancer. N Engl J Med 2021;384:440-51. [Crossref] [PubMed]
27. Conner SJ, Borges HB, Guarin JR, et al. Obesity Induces Temporally Regulated Alterations in the Extracellular Matrix That Drive Breast Tumor Invasion and Metastasis. Cancer Res 2024;84:2761-75. [Crossref] [PubMed]
28. Soldado-Gordillo A, Álvarez-Mercado AI. Epigenetics, Microbiota, and Breast Cancer: A Systematic Review. Life (Basel) 2024;14:705. [Crossref] [PubMed]
29. Kharazi DM, Karam L, El Boustany C, et al. Sodium butyrate and sodium propionate inhibit breast cancer cell migration and invasion through regulation of epithelial-to-mesenchymal transition and suppression of MEK/ERK signaling pathway. Front Cell Dev Biol 2025;13:1535563. [Crossref] [PubMed]
30. Gutierrez-Martinez VD, León-Del-Río A, Camacho-Luis A, et al. Uncovering a novel mechanism: Butyrate induces estrogen receptor alpha activation independent of estrogen stimulation in MCF-7 breast cancer cells. Genet Mol Biol 2024;47:e20230110. [Crossref] [PubMed]
31. Huse K, Bai B, Hilden VI, et al. Mechanism of CD79A and CD79B Support for IgM+ B Cell Fitness through B Cell Receptor Surface Expression. J Immunol 2022;209:2042-53. [Crossref] [PubMed]
32. Sakatani A, Igawa T, Okatani T, et al. Clinicopathological significance of CD79a expression in classic Hodgkin lymphoma. J Clin Exp Hematop 2020;60:78-86. [Crossref] [PubMed]
33. Meinberger D, Hermes G, Brachvogel B, et al. CLEC3A-Derived Antimicrobial Peptides Reduce Staphylococcus aureus Bacterial Counts in an In Vivo Biomaterial-Associated Infection Mouse Model. Pharmaceutics 2025;17:234. [Crossref] [PubMed]
34. Chen X, Wang Y, Li Y, et al. Identification of immune-related cells and genes in the breast invasive carcinoma microenvironment. Aging (Albany NY) 2022;14:1374-88. [Crossref] [PubMed]
35. Chen X, Sun H, Yang C, et al. Bioinformatic analysis and experimental validation of six cuproptosis-associated genes as a prognostic signature of breast cancer. PeerJ 2024;12:e17419. [Crossref] [PubMed]
36. Zhou H, Xi Y, Chen X. Chloride intracellular channel 6 inhibits hepatocellular carcinoma progression by modulating immune cell balance and promoting tumor cell apoptosis. Cytojournal 2025;22:20. [Crossref] [PubMed]
37. Liu J, Wei T, Quan L, et al. Identification and multi-dimensional validation of mitochondrial permeability transition-driven necrosis-related model to assess the prognosis and immunotherapy value in breast cancer. Eur J Med Res 2025;30:113. [Crossref] [PubMed]
38. Liu S, Feng L, Wang Z. DCTPP1: A promising target in cancer therapy and prognosis through nucleotide metabolism. Drug Discov Today 2025;30:104348. [Crossref] [PubMed]
39. Feng L, Wang X, Guo X, et al. Identification of Novel Target DCTPP1 for Colorectal Cancer Therapy with the Natural Small-Molecule Inhibitors Regulating Metabolic Reprogramming. Angew Chem Int Ed Engl 2024;63:e202402543. [Crossref] [PubMed]
40. Fu R, Chen Y, Zhao J, et al. The signature of SARS-CoV-2-related genes predicts the immune therapeutic response and prognosis in breast cancer. BMC Med Genomics 2024;17:260. [Crossref] [PubMed]
41. Muñoz JP, Soto-Jiménez D, Calaf GM. DCTPP1 Expression as a Predictor of Chemotherapy Response in Luminal A Breast Cancer Patients. Biomedicines 2024;12:1732. [Crossref] [PubMed]
42. Kafai NM, Janova H, Cain MD, et al. Entry receptor LDLRAD3 is required for Venezuelan equine encephalitis virus peripheral infection and neurotropism leading to pathogenesis in mice. Cell Rep 2023;42:112946. [Crossref] [PubMed]
43. Li Y, Qin G, Du J, et al. circRNA LDLRAD3 Enhances the Malignant Behaviors of NSCLC Cells via the miR-20a-5p-SLC7A5 Axis Activating the mTORC1 Signaling Pathway. J Healthc Eng 2022;2022:2373580. Retraction in: J Healthc Eng 2023;2023:9756529. [Crossref] [PubMed]
44. Zarrizi R, Higgs MR, Voßgröne K, et al. Germline RBBP8 variants associated with early-onset breast cancer compromise replication fork stability. J Clin Invest 2020;130:4069-80. [Crossref] [PubMed]
45. Qi H, Wu Y, Zhang W, et al. The syntaxin-binding protein STXBP5 regulates progerin expression. Sci Rep 2024;14:23376. [Crossref] [PubMed]
46. Wang J, Yang L, Li Y, et al. Aberrant hypermethylation induced downregulation of antisense lncRNA STXBP5-AS1 and its sense gene STXBP5 correlate with tumorigenesis of glioma. Life Sci 2021;278:119590. [Crossref] [PubMed]
47. Park JE, Kim HW, Yun SH, et al. Ginsenoside Rh2 upregulates long noncoding RNA STXBP5-AS1 to sponge microRNA-4425 in suppressing breast cancer cell proliferation. J Ginseng Res 2021;45:754-62. [Crossref] [PubMed]
48. Xu C, Xie Y, Xie P, et al. ZDHHC9: a promising therapeutic target for triple-negative breast cancer through immune modulation and immune checkpoint blockade resistance. Discov Oncol 2023;14:191. [Crossref] [PubMed]
49. Sadeghi M, Dehnavi S, Sharifat M, et al. Innate immune cells: Key players of orchestra in modulating tumor microenvironment (TME). Heliyon 2024;10:e27480. [Crossref] [PubMed]
50. Yu J, Fu L, Wu R, et al. Immunocytes in the tumor microenvironment: recent updates and interconnections. Front Immunol 2025;16:1517959. [Crossref] [PubMed]
51. Li X, Qu X, Li S, et al. Development of a Novel CD8(+) T Cell-Associated Signature for Prognostic Assessment in Hepatocellular Carcinoma. Cancer Control 2024;31:10732748241270583. [Crossref] [PubMed]
52. Li Y, Zhang C, Jiang A, et al. Potential anti-tumor effects of regulatory T cells in the tumor microenvironment: a review. J Transl Med 2024;22:293. [Crossref] [PubMed]