Abstract

Introduction

Breast cancer is a major global health issue, with approximately 2.3 million new cases in 2022, making it the second most common cancer after lung cancer[1] . In China, it ranks second in incidence and fifth in mortality among female cancers[2] . Despite improved diagnostics leading to lower mortality rates, breast cancer incidence continues to rise, particularly among younger women, largely due to aging populations and lifestyle changes[3] . This underscores the need for deeper understanding of its underlying molecular mechanisms.

Although treatment methods have advanced, chemotherapy resistance remains a significant clinical challenge[4] . 5-Fu, an essential chemotherapy drug, inhibits thymidylate synthase to disrupt DNA synthesis[5] . However, resistance to 5-Fu severely limits its effectiveness and contributes substantially to cancer-related mortality[6] . Thus, clarifying the mechanisms of 5-Fu resistance and identifying reliable biomarkers are crucial to improving patient prognosis.

Lipid metabolism reprogramming is a key driver of chemotherapy resistance, particularly to 5-Fu treatment, where cancer cells adapt by accumulating lipid droplets under chemotherapeutic stress[7] . For example, 5-Fu-resistant colorectal cancer cells show elevated sphingolipids and reduced ceramides mediated by SMPD1[8] . Additionally, lipid metabolic alterations can influence cell membrane fluidity, drug uptake, oxidative stress, and apoptosis pathways, collectively forming a protective mechanism against 5-Fu treatment[9] .Understanding these metabolic adaptations not only clarifies resistance mechanisms but also guides potential therapeutic strategies targeting lipid metabolism to enhance chemotherapy efficacy.

In this study, we integrated lipid metabolism-related genes from public databases to construct a prognostic model predicting 5-Fu resistance in breast cancer. Through differential expression analysis and Cox regression, we identified key prognostic genes and validated our risk-scoring model using ROC curve and Kaplan-Meier analyses in external cohorts. We further developed a nomogram to facilitate clinical use. Additionally, functional enrichment, protein-protein interaction analyses, and immune microenvironment comparisons were conducted to clarify the biological mechanisms of identified genes. Our model provides a valuable tool for individualized prognosis and therapeutic decisions in breast cancer.

Materials and Methods

Data Sources and Collection

The GSE243526 single-cell RNA sequencing dataset was obtained from the Gene Expression Omnibus (GEO)[10] to identify 5-Fu-resistant and -sensitive tumor subpopulations. Lipid metabolism-related genes (n = 16,338) were retrieved from GeneCards[11] . Bulk RNA-seq data, along with clinical features and survival information of breast cancer (BRCA) patients, were downloaded from The Cancer Genome Atlas (TCGA) using the easyTCGA package and served as the primary dataset for model construction. The external validation cohort, GSE103091, was also retrieved from GEO and included expression profiles and follow-up data. Detailed sample information is presented in Table 1.

Table 1. Sample Information from TCGA-BRCA and GEO Databases

Database	Information
GSE243526	Single-cell sequencing samples (n=16): breast cancer tissues (n=12) and normal breast tissues (n=4)
TCGA-BRCA	RNA-bulk-Seq samples (n=1231): tumor samples (n=1118) and normal samples (n=113) Clinical characteristics of tumor samples（n=906）
GSE103091	Expression profiling microarray samples with survival information (n=107)

Single-Cell RNA Sequencing Dataset Analysis

The GSE243526 single-cell RNA sequencing dataset was analyzed to identify gene expression differences associated with 5-Fu resistance. Quality control was conducted using the Seurat V5 framework. Cells were retained if they had between 300 and 7000 detected genes, less than 20% mitochondrial gene expression, less than 3% hemoglobin gene expression, and sequencing depth above 1000. Batch effects were corrected using the Harmony algorithm.

Cell type annotation was based on cluster-specific marker genes. Malignant epithelial cells were identified by comparison with normal myoepithelial cells. Chromosomal copy number variations (CNVs) were evaluated using the inferCNV algorithm. CNV scores were calculated by normalizing gene-level signals to the [−1,1] interval and summing the squared CNV values across all genes in each cell. Uniform Manifold Approximation and Projection (UMAP) was used for dimensionality reduction and visualization of cell subpopulations.

Following reclustering of malignant epithelial cells, the Scissor algorithm was applied to classify 5-Fu-resistant and -sensitive subpopulations. Scissor integrates single-cell expression profiles, bulk transcriptomic data, and phenotypic responses such as drug sensitivity to identify phenotype-associated cell populations[12] . Specifically, the single-cell expression matrix of malignant epithelial cells, bulk RNA-seq data from the Genomics of Drug Sensitivity in Cancer (GDSC) database, and IC50 values for 5-Fu were jointly used. Cells identified as Scissor+ were labeled as resistant (epi_5-fu_Res), whereas Scissor− and unclassified cells were grouped as sensitive (epi_5-fu_Sen).

Differential gene expression between the resistant and sensitive groups was analyzed using the FindMarkers function in Seurat. Genes were retained if they showed |log2 fold change| > 0.25, were expressed in at least 25% of cells, had an intergroup expression proportion difference >10%, and a p-value < 0.05 by Wilcoxon rank-sum test.

DEGs Identification and Core Gene Set Screening

Differential expression analysis was performed on the TCGA-BRCA dataset using the DESeq2 package, with significant genes defined as those with |log₂FC| > 0.5 and p < 0.05. A multi-set intersection was then conducted among three gene sets: differentially expressed genes from TCGA-BRCA, 5-Fu resistance–associated genes identified from single-cell analysis, and lipid metabolism–related genes. Overlapping genes were considered candidates related to both 5-Fu resistance and lipid metabolism. The intersections were visualized using Venn diagrams generated with the VennDiagram package.

Subsequently, functional enrichment analysis was conducted on the candidate genes using the clusterProfiler package. GO terms and KEGG pathways with p-values < 0.05 were considered significant, and the top 10 terms were visualized using ggplot2. Protein-protein interaction (PPI) analysis was conducted via the STRING database[13] with an interaction score cutoff of ≥ 0.4. The PPI network was visualized using Cytoscape, with node size and color intensity representing protein importance..

Development and Validation of a Prognostic Risk Signature

Univariate Cox regression was conducted using the survival package to identify genes significantly associated with 5-Fu resistance and lipid metabolism (p < 0.05). Genes meeting the proportional hazards (PH) assumption, assessed via Schoenfeld residuals (p > 0.05), were further included in multivariate Cox regression. PH assumption was re-evaluated in the multivariate model, and independent prognostic factors (p < 0.05) were visualized in forest plots using the forestploter package, showing hazard ratios (HR) and 95% confidence intervals.

A prognostic risk score was constructed using coefficients from multivariate Cox regression: Risk score=(−0.224×PDLIM4)+(0.166×SDC1)+(0.204×EMP1).Patients in both the TCGA-BRCA training set and GSE103091 validation set were divided into high- and low-risk groups based on the median risk score. Overall survival (OS) differences between groups were assessed using Kaplan-Meier analysis and log-rank tests via the survminer package. The predictive accuracy of the signature at 1-, 3-, and 5-year intervals was evaluated using time-dependent ROC curves generated with the survivalROC package.

Gene Expression and Survival Risk Analysis

To evaluate the prognostic model and risk stratification effectiveness, several visualization methods were applied. Risk score distribution and survival status were visualized using the survminer and forestploter packages, allowing assessment of patient distribution and corresponding outcomes. Heatmaps of the key genes in the model were generated with the pheatmap package, highlighting differential expression patterns between high- and low-risk groups.

Additionally, principal component analysis (PCA) was performed for three-dimensional visualization. Samples were projected onto the space defined by the top three principal components to illustrate the separation and clustering of risk groups in multidimensional gene expression space.

Construction of a Nomogram Model

Patients with complete clinical and survival data were selected from the TCGA-BRCA cohort. Risk scores and clinical variables were analyzed using univariate Cox regression to identify factors significantly associated with overall survival (p < 0.05), and results were visualized in forest plots via the survminer and forestploter packages. Variables satisfying the proportional hazards (PH) assumption, assessed using Schoenfeld residuals, were included in a multivariate Cox regression model. PH assumptions were reassessed to ensure model stability, and independent prognostic factors were identified.

A nomogram was constructed using the rms and regplot packages based on these independent predictors. Model performance was evaluated using time-dependent ROC curves (via survivalROC) at 1-, 3-, and 5-year intervals, calibration plots, and decision curve analysis (DCA) with the ggDCA package, confirming its potential for clinical application in breast cancer prognosis.

Immune Correlation Analysis

To evaluate the relationship between risk stratification and the tumor immune microenvironment, immune infiltration profiles were compared between high- and low-risk groups in the TCGA-BRCA cohort. Using the IOBR package [14] , the relative proportions of 22 immune cell types were quantified via the CIBERSORT algorithm. Immune cell compositions across samples were visualized using stacked bar plots to display overall differences between risk groups.

Spearman correlation analysis was performed to assess interactions among immune cells, and the results were visualized as correlation heatmaps using the corrplot package. Correlation strength (|r| ≥ 0.3) and significance levels (p < 0.05, 0.01, 0.001) were annotated. Differences in immune cell abundance between risk groups were evaluated using the Wilcoxon rank-sum test, with results presented in box plots.

The associations between prognostic genes (PDLIM4, SDC1, and EMP1) and immune cell infiltration levels were further analyzed via Spearman correlation, and results were illustrated in heatmaps. Additionally, correlation coefficients between risk scores and immune cell proportions were calculated to identify immune cell types significantly associated with the risk score (|r| > 0.3, p < 0.05), suggesting potential biological links between the risk model and specific immune features.

Functional Enrichment and Pathway Activity Analysis

To explore molecular pathway activity across different risk groups, functional enrichment analysis was performed. Patients from the TCGA-BRCA cohort were stratified into high- and low-risk groups based on the prognostic risk score. Differentially expressed genes (DEGs) between the two groups were identified using DESeq2, with cutoffs of |log2FC| > 0.5 and p < 0.05. Results were visualized as volcano plots using ggplot2 and ggrepel.

KEGG enrichment analysis was conducted separately for upregulated and downregulated DEGs using clusterProfiler and org.Hs.eg.db. The top 10 enriched pathways and their associated genes were visualized using the GOChord function from the GOplot package.

Additionally, Gene Set Enrichment Analysis (GSEA) was performed based on ranked genes calculated as −log10(p-value) × sign(log2FC). The MsigDB KEGG gene set (c2.cp.kegg.v7.4.symbols.gmt) was used for reference. GSEA results were visualized via the gseaNb function in the GseaVis package, and significantly enriched pathways in the high-risk group were defined as those with NES > 1, p < 0.05, and FDR < 0.25.

Protein Expression Validation and Functional Network Analysis

To validate the expression of prognostic genes at the protein level, immunohistochemical data were obtained from the Human Protein Atlas (HPA) database. Protein expression patterns in normal and breast cancer tissues were compared to evaluate their potential roles in tumor development. Functional network analysis was performed using GeneMANIA [15] to investigate biological associations of the prognostic genes. This analysis incorporated multiple data types, including physical interactions, co-expression, co-localization, pathway co-occurrence, and protein domain similarity, to construct a functional association network. Key functional partners were identified, along with their involvement in major biological processes, molecular functions, and cellular components..

Statistical Analysis

All statistical analyses and visualizations were performed using R (version 4.4.0) and RStudio (version 2024.12.1+563), along with relevant online databases. A two-sided p-value < 0.05 was considered statistically significant. Significance levels were indicated as follows: p < 0.05 (*), p < 0.01 (**), p < 0.001 (***), and p < 0.0001 (****).

Results

Identification of Differential Genes Associated with 5-Fu Resistance through Single-Cell Data Analysis

Differential gene expression associated with 5-FU resistance was investigated using the GSE243526 single-cell RNA sequencing dataset. After rigorous quality control, 123,901 high-quality cells were retained for analysis (Figure 1A). Batch effects were corrected through multi-sample integration, confirmed by UMAP visualization (Figure 1B).

Cells were clustered using a resolution of 0.2, resulting in 15 initial clusters (Figure 1C). Based on marker gene expression, 10 distinct cell types were annotated (Figure 1D–E), including B cells, T cells, macrophages, mast cells, fibroblasts, endothelial cells, pericytes, luminal epithelial cells, myoepithelial cells, and cancer cells. Cancer cells were identified via chromosomal copy number variation (CNV) analysis using inferCNV.

Subsequent analysis of epithelial subpopulations (Supplementary Figure 1A–F) refined luminal epithelial clusters into 16 subgroups. Subgroups 0, 2, 7, 9, and 15 exhibited CNV profiles similar to normal epithelial tissue and were classified as normal luminal cells. Remaining clusters were deemed malignant, while mixed subgroups (5 and 13) were excluded to improve clustering accuracy.

Reclustering of cancer cells revealed 11 subgroups (Figure 2A–B). Using the Scissor algorithm, cells were classified into 5-FU-resistant and -sensitive subpopulations based on similarity to bulk expression and drug sensitivity data (Figure 2C). Their spatial distribution was visualized via UMAP (Figure 2D). Differential expression analysis between resistant and sensitive clusters identified 84 5-FU resistance–associated genes (Figure 2E).

Identification of Candidate Genes Linking 5-FU Resistance to Lipid Metabolism

Differential expression analysis of the TCGA-BRCA dataset identified 5,985 upregulated and 3,591 downregulated genes in breast cancer samples. These results were visualized using a volcano plot (Figure 3A) and a heatmap of the top 10 upregulated and downregulated genes (Figure 3B).

A multi-set intersection analysis was performed among three gene sets: differentially expressed genes from TCGA-BRCA, 5-FU resistance–associated genes identified from single-cell analysis, and lipid metabolism–related genes from GeneCards. This integrative approach identified 42 candidate genes associated with both 5-FU resistance and lipid metabolism (Figure 3C).

Protein-protein interaction (PPI) analysis of these 42 genes was conducted using the STRING database (Figure 3D). The resulting network revealed HSPA1A and FOS as central hub proteins, interacting with multiple partners, suggesting their key regulatory roles in lipid metabolism and 5-FU resistance.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were then performed. GO analysis revealed 23 enriched pathways, with the top 10 involving processes such as response to metal ions, protein folding, and cAMP signaling (Figure 3E). KEGG analysis identified 12 pathways, with the top 10 enriched in estrogen signaling, proteoglycans in cancer, and cell adhesion molecules (Figure 3F).

Identification of Prognostic Genes and Construction of Risk Prediction Model

Univariate Cox regression analysis identified seven prognostic gene candidates (Figure 4A). Following proportional hazards (PH) assumption testing, HSPB8 and EFHD2 were excluded due to non-compliance (p < 0.05; Figure 4B). Multivariate Cox regression on the remaining five genes identified PDLIM4, SDC1, and EMP1 as independent prognostic factors (igure 4C), all of which met the PH assumption (Figure 4D).

A prognostic risk model was developed using the expression levels and regression coefficients of these three genes. Risk scores were calculated for each patient, and individuals were stratified into high- and low-risk groups based on the median score. Receiver operating characteristic (ROC) curve analysis showed good predictive performance in the TCGA-BRCA training set, with AUCs of 0.631, 0.662, and 0.641 at 1, 3, and 5 years, respectively (Figures 4E–F). In the external validation cohort GSE103091, AUCs were 0.822, 0.65, and 0.69 at corresponding time points. These results demonstrate the model’s stable discriminative power across short- and long-term survival. Kaplan–Meier survival analysis confirmed significantly shorter overall survival in the high-risk group compared to the low-risk group in both the training (p = 0.00027) and validation (p = 0.047) sets (Figures 4G–H), as determined by log-rank test.

Further visualization of the model included risk score distribution, survival status, gene expression heatmaps, and principal component analysis (PCA). Risk score plots (Figures 5A–B) showed a gradient distribution from low- to high-risk. Survival scatter plots (Figures 5C–D) indicated a higher concentration of death events in the high-risk group. Heatmaps (Figures 5E–F) revealed distinct expression patterns of the three prognostic genes. Three-dimensional PCA (Figures 5G–H) confirmed clear separation between high- and low-risk groups, validating the model’s discriminative capacity and generalizability across independent datasets.

Construction and Validation of a Nomogram Model

To evaluate the risk score as an independent prognostic factor in breast cancer, univariate and multivariate Cox regression analyses were performed using the complete TCGA cohort. Univariate analysis (Figure 6A) identified T stage, N stage, M stage, age, pathological stage, and risk score as significant predictors of overall survival (p < 0.001), while molecular subtype was not statistically significant (p = 0.063). Proportional hazards (PH) assumption testing (Figure 6B) indicated that T stage and pathological stage violated the assumption (p < 0.05). Multivariate Cox analysis (Figure 6C) confirmed that N stage, M stage, age, and risk score were independent prognostic factors, all meeting the PH assumption (Figure 6D).

A Nomogram was constructed incorporating age, N stage, M stage, and risk score to predict 1-, 3-, and 5-year overall survival (Figure 6E). Calibration plots demonstrated good agreement between predicted and observed survival outcomes (Figure 6F). Time-dependent ROC curves showed AUCs of 0.797, 0.757, and 0.766 for 1-, 3-, and 5-year survival, respectively, indicating favorable discrimination performance (Figure 6G). Decision curve analysis (Figure 6H) confirmed that the Nomogram achieved stable net clinical benefit across all time points, outperforming treat-all or treat-none strategies within clinically relevant risk thresholds.

In summary, the Nomogram model exhibited robust performance in both short-term and long-term survival prediction, providing a valuable tool for individualized prognostic assessment in breast cancer patients.

Molecular Pathway Signatures in Different Risk Groups

To explore molecular differences between risk groups, differential expression and pathway enrichment analyses were performed. A total of 407 genes were significantly upregulated and 133 downregulated in the high-risk group. The top 20 differentially expressed genes were visualized in a volcano plot (Figure 7A).

KEGG pathway enrichment analysis revealed that upregulated genes in the high-risk group were primarily involved in neuroactive ligand–receptor interaction, bile secretion, and cytochrome P450–mediated drug metabolism (Figure 7B). In contrast, downregulated genes were enriched in neuroactive ligand–receptor interaction and cytokine–cytokine receptor interaction pathways (Figure 7C).

GSEA further showed significant activation of amino acid metabolism pathways, including tyrosine, histidine, and phenylalanine metabolism, in the high-risk group (Figure 7D). These findings suggest that metabolic pathway reprogramming, particularly involving amino acid metabolism, may contribute to the aggressive phenotype and poorer prognosis observed in high-risk breast cancer patients.

Association Between Risk Score and Tumor Immune Microenvironment

The relationship between the risk score and tumor immune microenvironment was analyzed. Immune cell infiltration was quantified using the CIBERSORT, revealing the proportions of 22 immune cell types (Figure 8A). Comparative analysis (Figure 8B) showed that resting dendritic cells and M1 macrophages were significantly more abundant in the low-risk group (p < 0.05).

Immune cell interaction analysis (Figure 8C) revealed a strong positive correlation between activated CD8⁺ T cells and CD4⁺ memory activated T cells (r = 0.439, p < 0.001), while M0 macrophages exhibited significant negative correlations with CD4⁺ memory resting T cells (r = −0.404) and CD8⁺ T cells (r = −0.368), suggesting a dynamic balance between pro- and anti-tumor immune components.

Correlation analysis between prognostic genes and immune cell types (Figure 8D) showed that SDC1 was positively correlated with M0 macrophages and neutrophils; PDLIM4 was negatively correlated with CD4⁺ memory activated T cells and M2 macrophages; and EMP1 was positively correlated with CD4⁺ memory resting T cells but negatively correlated with regulatory T cells. Although statistically significant, the correlation coefficients were all below 0.3, indicating modest but potentially biologically relevant immune associations of the prognostic genes.

Protein Expression Validation and Functional Network Analysis

Immunohistochemistry data from the Human Protein Atlas (HPA) were used to evaluate the protein-level expression of prognostic genes in breast cancer. As shown in Figure 9A, PDLIM4 exhibited higher expression in tumor tissues compared to normal tissues, whereas SDC1 showed a notable downregulation. No immunohistochemical data were available for EMP1 in the HPA database.

A functional interaction network was constructed using GeneMANIA (Figure 9B), revealing that the three prognostic genes were functionally connected with proteins involved in cell adhesion (e.g., SDCBP, CD9), immune regulation (e.g., ANXA1, CXCL2), and signal transduction (e.g., TRIP6, NQO1). These interactions were supported by multiple mechanisms, including physical binding, co-expression, and pathway co-occurrence. This network suggests potential roles for the prognostic genes in breast cancer progression through diverse molecular pathways.

Discussion

Chemotherapy resistance remains a major obstacle in the effective management of breast cancer, particularly in aggressive subtypes such as triple-negative breast cancer (TNBC)[16] . This resistance is driven by multifactorial mechanisms, including genomic instability, epigenetic alterations, metabolic reprogramming, and changes in the tumor microenvironment [17] . Immune-related signaling pathways have also been implicated in chemoresistance, with previous studies highlighting the roles of genes such as PRC1, GGTLC1, and IRS1[18] . Moreover, chemoresistance has been associated with shortened patient survival and remodeling of the immune landscape, which together contribute to reduced treatment efficacy [19] . Given these challenges, elucidating the molecular mechanisms underlying 5-FU resistance and identifying robust predictive biomarkers are critical for advancing precision medicine in breast cancer.

Recent studies have highlighted lipid metabolism reprogramming as a key contributor to chemotherapy resistance. Breast cancer cells exhibit unique lipid metabolic profiles influenced by their tissue microenvironment and endocrine signaling [20] . For instance, estrogen receptor–positive breast cancer promotes fatty acid synthesis to modulate membrane fluidity [21] , whereas triple-negative breast cancer relies on lipophagy and lipid droplet dynamics to maintain energy homeostasis under stress [7] . These distinct metabolic adaptations underscore the pivotal role of lipid metabolism in shaping chemotherapeutic responses and influencing patient prognosis.

In this study, a lipid metabolism–related prognostic model for 5-FU resistance in breast cancer was developed by integrating single-cell transcriptomic and bulk multi-omics data. A total of 42 candidate genes associated with both 5-FU resistance and lipid metabolism were identified through differential expression analysis combined with the Scissor algorithm and TCGA-BRCA dataset. Based on Cox regression analysis, a three-gene risk signature comprising PDLIM4, SDC1, and EMP1 was constructed and validated in independent cohorts. The risk score model effectively stratified patients by prognosis, with high-risk patients exhibiting significantly shorter overall survival. Importantly, the risk score remained an independent prognostic factor after adjustment for clinicopathological variables. A nomogram incorporating the risk score and clinical features was also established, demonstrating good predictive performance and clinical applicability. This integrative approach provides a novel framework for predicting chemoresistance and offers a potentially valuable tool for individualized treatment decision-making in breast cancer.

The three genes identified in our risk model—PDLIM4, SDC1, and EMP1—have all been previously implicated in tumor progression and therapy resistance.

PDLIM4, a member of the PDZ and LIM domain protein family, regulates cytoskeletal organization and signal transduction[22] . Its aberrant methylation has been associated with drug resistance in multiple malignancies. In chronic myeloid leukemia, hypermethylation of PDLIM4 correlates with imatinib resistance and disease progression[23] .In prostate cancer, it has been reported as a prognostic marker linked to aberrant promoter methylation [24] . In breast cancer, PDLIM4 overexpression enhances migratory and metastatic capacity, potentially acting as an oncogene independent of Src signaling [25] . These findings align with our results, where PDLIM4 was upregulated in high-risk patients and associated with poor prognosis.

SDC1 encodes a transmembrane heparan sulfate proteoglycan involved in cell proliferation, adhesion, migration, and angiogenesis [26] . Immunohistochemical data from the Human Protein Atlas confirmed reduced SDC1 expression in breast cancer tissues in our study. Interestingly, SDC1 exhibits context-dependent functions: in ER-positive breast cancer, it promotes cell migration and correlates with poor prognosis [27] ; in pancreatic cancer, it facilitates epithelial–mesenchymal transition (EMT) and tumor progression[28] . In contrast, SDC1 downregulation has been observed in gastric[29] and colorectal cancers[30] , suggesting tumor-type specific roles. Our data support a tumor-suppressive role of SDC1 in breast cancer, consistent with its lower expression in the high-risk group.

EMP1 encodes an epithelial membrane protein associated with inflammation and epithelial pathologies[31] . It has been shown to promote triple-negative breast cancer progression by enhancing IL-6 secretion via NF-κB signaling, which stimulates tumor-associated fibroblast proliferation and metastatic spread[32] . The inclusion of EMP1 in our prognostic model and its correlation with poor survival underscore its potential as a therapeutic target in aggressive breast cancer subtypes.

Pathway enrichment analysis revealed that high-risk patients exhibited activation of bile secretion and cytochrome P450–mediated drug metabolism pathways, both of which are involved in lipid processing and drug detoxification [33,34] . These findings are consistent with previous studies suggesting that altered lipid metabolism contributes to chemoresistance through enhanced metabolic flexibility and reduced drug efficacy[35] . Additionally, enrichment of amino acid metabolism pathways in the high-risk group indicates that metabolic reprogramming may provide survival advantages to tumor cells under chemotherapeutic stress, further reinforcing the role of metabolic adaptation in treatment resistance.

Our analysis revealed significant associations between the risk score and components of the tumor immune microenvironment. CIBERSORT analysis demonstrated that resting dendritic cells and M1 macrophages were significantly enriched in the low-risk group. Dendritic cells play a key role in antigen presentation and activation of immune effector cells, thereby contributing to tumor immune surveillance [36] . However, their function is often suppressed within the tumor microenvironment, a state that can be reversed through immune-stimulatory interventions such as CpG oligodeoxynucleotide treatment[37] . M1 macrophages support anti-tumor immunity by promoting Th1-type responses and secreting pro-inflammatory cytokines such as IL-12 and TNF-α[38] . Their enrichment in the low-risk group may indicate a more active immune surveillance environment, which could contribute to the improved prognosis observed in these patients.

Although the correlations between prognostic gene expression and immune cell infiltration were modest, they were statistically significant. These findings suggest that the identified genes may influence the tumor immune landscape through complex regulatory networks involving cytokine signaling, cell adhesion, and metabolic pathways. Further investigations using single-cell sequencing and spatial transcriptomics are warranted to delineate these interactions and guide the development of precise immunotherapeutic strategies.

This study has several limitations that should be acknowledged. First, it is a retrospective analysis based on publicly available datasets, which may introduce selection bias and limit the generalizability of the findings to broader patient populations. The heterogeneity in sequencing platforms, sample processing protocols, and clinical annotations across datasets could influence the robustness of the integrated analysis. Second, the predictive model has not yet been validated in prospective clinical trials. Experimental validation is essential to confirm the biological roles of the three identified prognostic genes—PDLIM4, SDC1, and EMP1—in mediating 5-FU resistance and lipid metabolism. Future studies should incorporate in vitro and in vivo assays to elucidate the functional relevance of these genes. Third, the immune microenvironment analysis in this study was primarily based on computational inference using CIBERSORT. While useful, such approaches rely on transcriptomic data and may not fully capture spatial or functional heterogeneity. Experimental validation, such as immunohistochemistry, multiplex immunofluorescence, or spatial transcriptomics, is needed to corroborate these findings. Lastly, the prognostic model was developed based on overall survival as the sole clinical endpoint. Future research should evaluate additional endpoints such as progression-free survival, treatment response, and recurrence risk, and test the model's predictive power across different breast cancer subtypes and treatment regimens.Together, these limitations highlight the need for prospective, multi-dimensional validation to ensure the clinical applicability and translational potential of the proposed risk model.

Conclusion

This study highlights the critical role of lipid metabolism in breast cancer resistance to 5-FU. We constructed a prognostic model based on PDLIM4, SDC1, and EMP1, which effectively stratified patient outcomes and demonstrated independent prognostic value. These findings offer a clinically relevant tool for personalized treatment and suggest that targeting lipid metabolism may enhance therapeutic efficacy in resistant breast cancer.

Supplementary Material

Supplementary methods, results, spectra, figures.

Acknowledgements

The authors are grateful for the valuable contribution of the TCGA and GEO project.

Author contributions

Zhihao Wei, Bin Ye and Zezhong Zhang conceptualized and designed the study, performed the majority of the bioinformatic analyses and data interpretation, and drafted the manuscript. Zhuojun Luo conducted investigation procedures and managed data curation. Jie Zhou supervised the project, provided resources, administered project execution, and critically reviewed and edited the manuscript. All authors read and approved the final manuscript.

Funding information

This work was supported by Guangzhou Science and Technology Project（2024A03J0649). 

Competing Interests

The authors declare that they have no existing or potential commercial or financial relationships that could create a conflict of interest at the time of conducting this study.

Data Availability

All data needed to evaluate the conclusions in the paper are present in the paper or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

References