Cuprotosis-related gene subtypes, prognostic modeling, and tumor microenvironment remodeling in breast cancer.

Introduction
The most common malignancy in women is breast cancer (BC). According to the Global Cancer Statistics in 2022, there are 2.3 million new cases of female BC and 680,000 deaths from BC, making it ranking first among female malignant tumors (1). Although surgery, chemotherapy, and radiotherapy are effective treatments for BC, the morbidity and mortality of BC remain high. Therefore, it is imperative to explore new prognostic predictors and potential therapeutic targets. Tsvetkov et al. (2) identified a novel mode of cell death regulated by copper, which is different from the known mechanism of cell death and dependent on mitochondrial respiration, and indicated that this new mode of death occurs through direct binding of copper to lipoylated components of the tricarboxylic acid (TCA) cycle. This process leads to aliphatic protein aggregation and subsequent loss of iron-sulfur cluster protein stability, causing toxic-stressed protein, eventually leading to cell death. This completely new mode of copper-dependent cell death is termed ‘cuprotosis’. Previous studies have shown copper overload or deficiency can lead to impaired cell function and eventually cell death (3). Previous studies have reported elevated copper levels in gallbladder, pancreatic, and colon cancers compared with normal tissues (4-7). At present, there are a few comprehensive reports on the mechanism of cuprotosis-related pattern-related genes (CRPRGs) in BC (8,9).
Current studies suggest that the immune microenvironment plays an important role in the occurrence and development of cancer (10-12). Prognostic models based on the immune microenvironment have been shown to be useful in evaluating the overall survival (OS) of a variety of cancers, including melanoma, colon cancer and hepatocellular carcinoma (13-15). However, the mechanism of the complex interaction between CRPRGs and BC immune microenvironment still needs to be further explored.
Therefore, in this study, we explored the potential biological roles and multi-omics analysis of CRPRGs in BC. Finally, a model that can predict the prognosis and efficacy of treatments for BC was constructed and validated. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2403/rc).

Methods

Data collection and processing
The gene expression data and corresponding clinical information of BC patients were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) and the GSE20685 dataset of the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). All mRNA expression data were normalized using quantile normalization. The intersection genes were obtained after looping the two sets of gene data using the ‘limma’ and ’sva’ R packages. Tumor tissue samples data to merge, through FPKM (fragments per kilobase of transcript per million) calculation, and transform it into transcripts per million (TPM) value, data values for ‘log2’, the standard for the truncation P value <0.05 and |log2 fold change (FC)| ≥1, followed by a correction of data in batches. Samples lacking survival data, with follow-up time less than 30 days, or with incomplete key clinical variables were excluded. Finally, 1,409 samples with necessary clinical information and transcriptome data were included in the study. We treated missing values as an “unknown” category in multivariable analyses rather than performing imputation. Clinical variables included age at diagnosis, follow-up time, survival status, and tumor stage. To further validate the protein expression levels of the eight identified prognostic CRPRG score, immunohistochemistry (IHC) staining images of BC tissues and normal breast tissues were retrieved and analyzed from the Human Protein Atlas (HPA). The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Differential analysis of CRGs in BC and tumor mutational burden (TMB) analysis
The differential analysis of cuprotosis-related genes (CRGs) was performed on the BC transcriptome gene expression data obtained from the TCGA database. The ‘limma’, ‘reshape2’ and ‘ggpubr’ R packages were used to analyze the differences between normal breast tissue samples and BC tissue samples. CRG expression differences were visualized using box plots. The gene mutation expression levels of CRGs were extracted from the BC gene mutation data in the TCGA database. Then, the expression of gene mutations related to cuprotosis in BC was visualized by drawing a Waterfall Plot through the ‘maftools’ R package.

Copy-number analysis of CRGs in BC
The frequency of gene copy number variation (CNV) is the distribution of copy number gain or loss based on the analysis of CRGs copy number matrix data. Then the ‘RCircos’ R package was used to draw the chromosome circle diagram of CRGs from the CRGs expression file, copy-number annotation file and human chromosome reference file.

Survival and prognosis analysis of CRGs in BC
Clinical survival data from TCGA and GEO databases were merged. The packages ‘limma’, ’survival’ and ’survminer’ were used to loop CRGs with the combined clinical survival data to find prognostic CRGs of BC. Then univariate Cox analysis and Kaplan-Meier analysis were used to draw a visual survival analysis curve to compare the survival difference between the high- and low-risk expression groups. Finally, the network prognosis relationship diagram was drawn.

Analysis of CRG subtypes in BC
Univariate Cox regression analysis was used to screen genes related to cuprotosis associated with prognosis, and ‘ConsensuClusterPlus’ R package was used to cycle all samples into 9 subtypes (Samples can be divided into up to 9 subtypes, all set to a maximum K value of 9). Survival analysis, principal component analysis (PCA), gene set variation analysis (GSVA) analysis and single-sample gene set enrichment analysis (ssGSEA) analysis were performed on the patients with different types. R-packages ‘limma’, ‘GSVA’, ‘GSEABase’, and ‘pheatmap’ were used for GSVA analysis and heat mapping to assess pathway enrichment for CRG subtypes. ssGSEA analysis was used to assess the degree of tumor immune cell infiltration (16). The R packages ‘reshape2’, ‘ggpubr’, ‘limma’, ‘GSEABase’ and ‘GSVA’ were cited to compare the composition of different immune cells between the different CRG subtypes.

Differentially expressed gene (DEG) identification and function notes
The ‘limma’ R package was used to extract intersection genes by cycle comparison of CRG subtypes, and the differences of intersection genes were analyzed. Variance analysis of the set filter conditions to logFC =0.585, difference ratio of 1.5 times, and the adjust P value is 0.05. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were used to explore the biological functions of DEGs in the CRG subtypes. Samples were genotyped on the basis of clinical survival data. Then, the DEG gene clusters were subjected to survival analysis and CRGs differential analysis.

Construction and validation of CRPRG score
All BC patients were randomly divided into training set and testing set in a 1:1 ratio. The least absolute shrinkage and selection operator (LASSO) regression analysis was performed on the training set to construct the LASSO regression model, and cross-validation was performed to find out the characteristic genes. Through the characteristics of the genetic coefficient of building model of Cox model to generate the optimal formula, and the formula was subsequently validated in the testing set. CRPRG score was calculated as following: CRPRG score = ∑ (Normalized Expression * Coefi) n. Gene expression values were normalized using log2(TPM+1) followed by z-score transformation. Detailed model coefficients and baseline hazard functions are provided in Appendix 1 to facilitate external validation. According to the median risk value, patients were divided into high-risk group and low-risk group in the total sample, training set and testing set, and Kaplan-Meier survival analysis was performed. The ‘ggalluvial’ package was used to visualize the relationships among CRG subtypes, DEG gene clusters, CRPRG score and survival status. To further explore the predictive value of prognostic risk score in BC. The ‘limma’ R package was used to analyze differences in risk scores for CRG subtypes and DEG gene clusters. The expression differences of CRGs in BC between high- and low-risk groups were analyzed to identify CRGs affecting prognosis. Risk scores and clinical information such as age, clinical stage, T stage, and N stage were included in univariate and multivariate Cox analyses. The ‘rms’ package was then used to construct a nomogram to predict the 1-, 3-, and 5-year OS of BC patients based on the magnitude of the regression coefficients of each variable. The receiver operating characteristic (ROC) curve was used to evaluate the accuracy of the nomogram, and the calibration curve was used to estimate the consistency.

Analysis of risk curves of prognostic model
Survival information of samples from high- and low-risk groups was extracted, and the samples were sorted according to the patients’ risk scores to draw the risk curve, and the survival state scatter plot was drawn according to the patients’ survival state and survival time. Characteristic genes and gene expression were extracted from the prognostic model to draw the risk heat map of the prognostic model.

Immunoinfiltration analysis and stem cell index of prognostic model
To evaluate the correlation between risk score and immune infiltration, the CIBERSORT algorithm was used to calculate the percentage of different immune cells and the variability of immune cells between high- and low-risk groups of BC. All risk data and stem cell data were interposed and plotted to investigate the stem cell relevance of the prognostic model.

TMB and drug sensitivity analysis
In order to analyze the difference of tumor mutation between high- and low-risk groups, the tumor mutation frequency of high- and low-risk group was analyzed by ‘maftools’ package in R software, and the waterfall diagram was drawn. To explore the differences in chemosensitivity of patients with different risk groups, we calculated the half-maximal inhibitory concentration (IC50) of commonly used chemotherapeutic drugs for BC using the ‘pRRophetic’ package.

Statistical analysis
All data in this study were statistically analyzed using R software (version 4.4.0). When the P value was less than 0.05, statistical tests were considered statistically significant. For multiple comparisons, P values were adjusted using the Benjamini-Hochberg false discovery rate (FDR) method. Hazard ratios (HRs) and area under the curve (AUC) values are reported with 95% confidence intervals (CIs).

Results

Expression and mutation of CRGs in BC
The BC transcriptome gene expression data downloaded from the TCGA database included tumor (n=1,106) and normal (n=112) samples. To compare the expression of 19 CRGs between tumor and normal tissues, differential expression analysis was performed by R software (Figure S1). The results of the difference analysis are shown in Figure 1A. The expressions of ATP7B, SLC31A1, PDHB and CDKN2A were significantly up-regulated in BC tissues. While NFE2L2, NLRP3, ATP7A, FDX1, LIAS, LIPT1, LIPT2, DLD, PDHA1, MTF1, GLS, DBT, GCSH, DLST is significantly downregulated. Next, we analyzed the relationship between tumor gene mutations and CRGs in BC. The results showed that 55 of the 981 samples had mutations in the CRGs. ATP7A had the highest mutation frequency at 2%. However, no mutations of FDX1, LIPT2, CDKN2A, DBT, or GCSH were found in all BC samples (Figure 1B). By analyzing the variation frequency of CRGs gene copy number, we found that the increase frequency of NLRP3, LIPT2, GLS, MTF1, NFE2L2, LIPT1, SLC31A1, PDHA1 and LIAS copy number was greater than the loss frequency, while DLD, DBT, GCSH, FDX1, ATP7B, CDKN2A, DLAT, DLST, PDHB showed decreased CNV (Figure 1C) and their CNV locations on chromatin were shown in Figure 1D. Through univariate Cox analysis and Kaplan-Meier analysis of clinical survival data, 17 CRGs related to BC survival prognosis were screened. Then, according to the correlation between CRGs and BC prognosis, the network relationship prognosis map was drawn (Figure 1E).

Identification of CRG subtypes in BC
Unsupervised cluster analysis was performed based on the expression profiles of 17 genes. According to the cumulative distribution function (CDF) curve of the consensus score, we chose k=2 as the best division to classify BC patients into subtype A (n=788) and subtype B (n=622) (Figure 2A). PCA showed significant differences between the two subtypes (Figure 2B). Survival analysis showed that the prognosis of subtype A was significantly better than that of subtype B (Figure 2C). In addition, we visualized the relationship between CRG subtypes, CRGs expression, and clinicopathological features using heat maps (Figure 2D).

Characteristics of the tumor microenvironment (TME) in CRG subtypes
GSVA analysis showed that subtype B was mainly enriched in arachidonic acid metabolism, metabolism of xenobiotics by cytochrome p450, drug metabolism-cytochrome p450 pathways. Subtype A was enriched primarily in the ErbB signaling pathway, neurotrophin signaling pathway, basal transcription factors, RNA degradation (Figure 3A). To further explore the relationship between different CRG subtypes and immune infiltration, the differences in immune cell content in BC samples of different subtypes were analyzed using ssGSEA. The results showed that thirteen immune cell types differed significantly between subtypes. Among them are activated CD4+
T cell, activated dendritic cell, gamma delta T cell, immature B cell, regulatory T cell, type 2 T helper cell in subtype A had significantly higher infiltration. While activated CD8+ T cell, CD56bright natural killer (NK) cell, CD56dim NK cell, eosinophil, monocyte, NK cell, neutrophil had significantly higher infiltration in the B subtype (Figure 3B).

Analysis of DEG gene clusters based on cuprotosis-related DEGs
In order to further explore the potential biological functions of CRG subtypes, ‘limma’ package was used to screen differential genes between subtype A and subtype B. The screening criteria were |log2FC| ≥1, FDR <0.05, and 334 DEGs were finally screened (Figure 4A). Then, GO and KEGG analysis were performed based on DEGs, and the results showed that DEGs were enriched in cell cycle, proteoglycans in cancer, cellular senescence, and efferocytosis pathway (Figure 4B,4C). Through unsupervised cluster analysis of the differential genes between the two CRG subtypes, three DEG gene clusters A-C were identified (Figure 4D). Survival analysis revealed significant differences in prognosis among the three subtypes, with gene subtype A having the worst prognosis and gene subtype B having the best prognosis (Figure 4E); 15 CRGs were differentially expressed in three different DEG gene clusters (Figure 4F).

Construction of a prognostic risk assessment model for BC
Clinical characteristics of the training set and testing set were shown in Table S1. Based on the DEGs of the three DEG gene clusters, LASSO regression analysis and Cox proportional hazards regression analysis were conducted to identify characteristic genes associated with BC prognosis (Figure 5A,5B). Ultimately, eight signature genes and their corresponding risk coefficients were determined and analyzed. Formulation of a prognostic model for risk score estimation: CRPRG score = 0.4449 × Expr(PANX1) + 0.4443 × Expr(GPR107) + 0.2276 × Expr(MTFR1) − 0.1919 × Expr(FKBP5) + 0.1577 × Expr(CD24) + 0.0647 × Expr(CRISP3) − 0.0753 × Expr(NPY1R) − 0.0578 × Expr(ADIRF).
The Kaplan-Meier survival analysis demonstrated that low-risk patients exhibited significantly higher survival rates compared to high-risk patients across all samples, including the training cohort and the validation cohort (Figure 5C-5E). In addition, the results revealed substantial differences in CRPRG scores between CRG subtypes, with subtype A exhibiting a significantly higher CRPRG score compared to subtype B (Figure 5F). Furthermore, significant variations in CRPRG scores were observed among different DEG gene clusters, where gene subtype B demonstrated the highest CRPRG score, while gene subtype C exhibited the lowest CRPRG score (Figure 5G). We utilized Sankey diagrams to visualize the relationships and variations among CRG subtypes, DEG gene clusters, CRPRG scores, and survival statuses (Figure 5H). Thirteen CRGs exhibited significant differences in expression levels between the high-risk and low-risk groups (Figure 5I).

Clinical correlation analysis of CRPRG score, development and validation of a nomogram
To examine the correlation between CRPRG score and other clinical characteristics, univariate Cox analysis revealed that age, clinical stage, tumor size, lymph node stage, and risk score were significantly associated with OS in BC patients (Figure 6A). Subsequent multivariate Cox analysis further demonstrated that age, clinical stage, and risk score served as independent prognostic factors for OS in BC patients (Figure 6B).
To facilitate the clinical prediction of the prognosis of BC patients using the CRPRG score, we developed a nomogram based on five key factors: age, clinical stage, tumor size, lymph node status, and risk score. This nomogram is designed to predict the 1-, 3-, and 5-year OS rates of BC patients (Figure 6C). The calibration plot showed that the prediction of OS by the nomogram was highly consistent with the actual results (Figure 6D). ROC curve analysis demonstrated that the nomogram model exhibited satisfactory AUC values (with 95% CIs) in the overall sample, training set, and validation set (Figure 6E-6G).

Relationship between CRPRG score and immune cells, mutations, cancer stem cell (CSC) index, drug sensitivity and histology analysis
To investigate the impact of immune cells on prognosis, the correlation between CRPRG and immune cells were systematically analyzed. The resulting correlation heatmap revealed that certain immune cell types were correlated with CRPRG (Figure 7). Specifically, B memory cells, monocytes, activated NK cells, and CD8+ T cells exhibited significant negative correlations with the risk score, whereas M0 macrophages and activated CD4 memory T cells demonstrated significant positive correlations (Figure S2).
By analyzing TCGA samples, it was demonstrated that the TMB value in the high-risk group was significantly greater than that in the low-risk group (Figure 8A). Subsequently, we conducted an analysis of the distribution of somatic mutations in both the high-risk and low-risk groups within the TCGA cohort. The top ten mutated genes in both the high-risk group and the low-risk group were identified as follows: PIK3CA, TP53, TTN, CDH1, GATA3, MUC16, KMT2C, MAP3K1, HMCN1, FLG (Figure 8B,8C). The frequency of TP53 mutation was higher in the high-risk group than in the low-risk group, while the frequency of PIK3CA mutation was higher in the low-risk group.
Pearson analysis revealed a positive correlation between CRPRG score and CSC index (R=0.13, P=1.5×10−5). This suggests that higher CRPRG score values are associated with more pronounced stem cell characteristics and a greater degree of cell differentiation (Figure 8D).
We subsequently evaluated the sensitivity of patients in the high-risk and low-risk groups to commonly used chemotherapy drugs for BC. The results demonstrated that the IC50 of epothilone-B, ABT-263, AKT-inhibitor-VIII, and vinorelbine were lower in the high-risk group, whereas the IC50 values of methotrexate, rapamycin, AG014699, mitomycin-C, temsirolimus, and elesclomol were lower in the low-risk group. These exploratory findings suggest that the CRPRG score might be associated with predicted sensitivity to specific chemotherapeutic agents, warranting further experimental validation (Figure S3).
In this study, the expression levels of CRPRG in normal and BC tissues were systematically analyzed using the IHC data from histological and pathological profiles available in the HPA database. The findings revealed that CD24, CRISP3, GPR107, MTFR1, and PANX1 exhibited significantly higher expression in BC tissues, whereas ADIRF and FKBP5 demonstrated markedly lower expression levels (Figure 9).

Discussion
Studies have demonstrated that modes of cell death encompass apoptosis, necrosis, autophagy, ferroptosis, and others (17,18). Furthermore, copper-dependent programmed cell death, termed “cuprotosis”, was reported in 2022 (2). This mode of cell death is implicated in various diseases and holds promise as a potential therapeutic strategy for tumors in the future. Additional studies have indicated a significant association between cuprotosis and immune cells (19). However, the role of CRPRG score in BC has not been comprehensively investigated. In this study, we utilized an open-access database to conduct bioinformatics multi-omics analyses of BC, elucidating its heterogeneity, constructing a gene prognostic model associated with cuprotosis, and exploring the relationship between cuprotosis and both clinical characteristics and the immune microenvironment.
Most prior studies have primarily concentrated on the impact of individual CRG on oncological processes (20,21). Consequently, the combined effects of multiple CRGs on disease prognosis and the immune microenvironment have not been systematically elucidated. This study comprehensively analyzed the global alterations of CRPRG score in BC at the genomic and transcriptomic levels, with supportive protein-level evidence from HPA IHC images. Two CRG subtypes were identified based on 17 differentially expressed CRGs. Subtype A exhibited a significantly better survival prognosis compared to subtype B. Furthermore, substantial differences were observed in the immune microenvironment and immune-related signaling pathways between the two subtypes. Based on the differential genes of the two CRG subtypes, three genetic subtypes were further delineated. Our findings indicate that CRPRG score can serve as an independent predictor for BC prognosis. Distinct patient immune statuses corresponded to varying CRPRG scores. BC patients with high and low CRPRG scores demonstrated significant differences in clinical characteristics, immune infiltration, mutations, CSC index, drug sensitivity, and histology. Finally, a nomogram was constructed to predict the prognosis of BC patients, incorporating age, tumor stage, and CRPRG score. The ROC curve and calibration curve demonstrated the nomogram’s robust predictive performance.
Among the eight characteristic genes comprising the CRPRG score, PANX1, GPR107, MTFR1, CD24, and CRISP3 were identified as risk factors, whereas FKBP5, NPY1R, and ADIRF exhibited protective roles. PANX1 is a transmembrane channel-forming protein involved in intercellular signaling, ATP release, and inflammatory response regulation (22). Recent studies have demonstrated that PANX1-mediated ATP release may contribute to the formation of an immunosuppressive microenvironment (23). Overexpression of PANX1 in BC correlates with epithelial-mesenchymal transition (EMT) phenotype shifts, which subsequently lead to adverse clinical outcomes (24). GPR107, a member of the G protein-coupled receptor family, participates in cellular endocytosis, signal transduction, and tumor-associated biological processes (25). Elevated expression of GPR107 has been observed in prostate cancer tissues and is associated with poor prognosis (26). MTFR1 serves as a key regulator of mitochondrial fission and contributes to apoptosis, energy metabolism, and cell migration by maintaining mitochondrial dynamic balance (27). High levels of MTFR1 in colon cancer cells interact with NEK1 to promote mitochondrial fusion and sustain mitochondrial function under glucose deficiency, thereby enhancing cellular tolerance to glucose deprivation (28). CD24 is a highly glycosylated cell surface protein. In breast tumor tissue, CD24 interacts with Siglec-10 expressed in tumor-associated macrophages, facilitating tumor immune escape (29). CRISP3, a secreted protein belonging to the cysteine-rich secretory protein family, regulates tumor progression through mechanisms involving cell adhesion, immune microenvironment modulation, and cell signaling pathways. Elevated CRISP3 expression has been linked to weakened anti-tumor immune responses (30). FKBP5, a member of the immunophilin family that binds to the immunosuppressive drug FK506, exhibits positive correlation with CD8+ T cell-related marker expression in luminal B BC and is associated with favorable patient prognosis (31). NPY1R, a major receptor of neuropeptide Y and a member of the G protein-coupled receptor family, plays a role in regulating cell proliferation, apoptosis, and energy metabolism (32). High NPY1R expression correlates with increased sensitivity to endocrine therapy and prolonged OS in hormone receptor-positive BC patients (33). ADIRF is a critical factor in adipogenesis and metabolic regulation, its deficiency in prostate cancer is associated with lymph node metastasis and poor survival, suggesting its potential as a tumor suppressor gene in prostate cancer (34). Collectively, these findings indicate that these genes influence cancer development. However, the underlying mechanisms of their synergistic action in BC remain unclear. In our study, these genes were significantly associated with BC prognosis and may represent influential elements in BC progression.
Despite the availability of numerous treatments for BC, including surgery, chemotherapy, radiotherapy, and targeted therapy, the prognosis for recurrent and metastatic advanced BC remains unfavorable (35). Although immunotherapy has gained increasing attention in the field of BC, significant heterogeneity persists in the prognostic outcomes of BC patients, underscoring the critical role of the TME in the initiation and progression of BC (36). In the correlation analysis between CRPRG score and immune cells, CRPRG score exhibited a significant association with M0 macrophages and activated CD4 memory T cells. M0 macrophages can be readily polarized into M2 macrophages within tumor tissues under the influence of cytokines derived from both tumor cells and T cells (37). M2 macrophages highly express IL-10 and TGF-β, which directly suppress the cytotoxic function of CD8+ T cells (38). Furthermore, M2 macrophages overexpress PD-L1 and VISTA (V-domain Ig suppressor of T cell activation) on their surface, which bind to PD-1 and other receptors on T cells, thereby inhibiting T cell activation and suppressing the immune response while promoting tumor progression (39,40). Additional studies have demonstrated that M2 macrophages support tumor angiogenesis and matrix remodeling by secreting factors such as VEGF and FGF, indirectly inhibiting anti-tumor immunity (41,42). Increased infiltration of M2 macrophages in the TME correlates with poor prognosis in BC patients. Recent research indicates that Th2 subsets within activated CD4 memory T cells secrete cytokines like IL-4, IL-5, and IL-13, directly contributing to the polarization of M2 macrophages (43). However, some activated memory CD4 T cells differentiate into Treg cells, which inhibit effector T cell function through the secretion of IL-10, TGF-β, and overexpression of CTLA-4 (44). M2 macrophages and activated CD4 memory T cells establish a positive feedback loop via cytokine signaling, collectively suppressing anti-tumor immune responses. These findings align with our results, indicating that patients with high CRPRG scores exhibit poor responses to immunotherapy and unfavorable prognoses.
Furthermore, although the 8 genes in our model are not core cuprotosis regulators, correlation analysis revealed significant associations between the CRPRG score and the expression of cuprotosis genes, suggesting that our signature captures the downstream biological consequences of disrupted copper homeostasis. While our prognostic model was validated in a split-cohort design, we further corroborated the expression patterns of the signature genes at the protein level using IHC data from the HPA, confirming their aberrant expression in BC tissues.
Due to the inconsistent recording of clinical information in TCGA and GEO, not all samples included in this study have complete molecular subtype [estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER2) status] or PAM50 classification data. Therefore, we were unable to adjust for molecular subtypes in the multivariate Cox regression analysis, nor could we analyze the distribution of receptor status in different risk groups. Although our model shows that it is independent of standard clinicopathological features such as age, clinical stage, T stage, and N stage, in future studies, it is still necessary to use prospective cohorts with complete molecular data to verify the prognostic independence of the CRPRG score in different intrinsic BC subtypes. The prognosis prediction model for BC developed in this study is entirely based on data from public databases. However, our study relies on retrospective bioinformatic analyses. The identified signature is correlative. Future experimental studies, such as FDX1 knockdown or copper ionophore sensitivity assays in co-culture systems, are essential to elucidate the precise mechanistic links between the CRPRG signature, cuprotosis induction, and immune microenvironment remodeling.

Conclusions
In our integrative analysis, we identified two CRG subtypes and three cuprotosis-related DEG gene clusters in BC, and the associated CRGs uncovered multi-dimensional regulatory mechanisms that shape the tumor immune landscape, drive clinicopathological heterogeneity, and dictate patient outcomes. The CRPRG score established through these subtypes demonstrated robust predictive accuracy, offering a novel framework for stratifying immune checkpoint responsiveness, guiding targeted therapeutic strategies, and optimizing personalized survival prediction in BC patients.

Supplementary
The article’s supplementary files as