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 ‘log₂’, the standard for the truncation P value <0.05 and |log₂ 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 log₂(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 (IC₅₀) 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).

Figure 1 Expression and mutation of CRGs in BC. (A) Differential expression of 19 CRGs in normal and breast cancer tissues from TCGA. (B) Mutation frequencies of CRGs in 981 samples from TCGA. (C) CNV frequencies of CRGs. (D) Locations of CNV alterations in CRGs on 23 chromosomes. (E) Relationship between CRGs and prognosis of breast cancer. *, P<0.05; **, P<0.01; ***, P<0.001. BC, breast cancer; CNV, copy number variation; CRG, cuprotosis-related gene; TCGA, The Cancer Genome Atlas; TMB, tumor mutational burden.

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).

Figure 2 Identification of CRG subtypes in BC. (A) Consensus matrix heatmap of the 2 CRG subtypes. (B) PCA showed significant differences between the two subtypes. (C) Differences in overall survival time between the 2 CRG subtypes. (D) Differences in clinicopathological features and CRG expression in different CRG subtypes. BC, breast cancer; CRG, cuprotosis-related gene; M, metastasis; N, node; PCA, principal component analysis; T, tumor; TCGA, The Cancer Genome Atlas.

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).

Figure 3 Characteristics of the tumor microenvironment in CRG subtypes. (A) GSVA of biological pathways between different CRG subtypes, in which red and blue represent activated and inhibited pathways, respectively. (B) Differences in infiltration of 22 immune cells between the two subtypes. *, P<0.05; **, P<0.01; ***, P<0.001. CRG, cuprotosis-related gene; GSVA, gene set variation analysis; MDSC, myeloid-derived suppressor cell; TCGA, The Cancer Genome Atlas.

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 |log₂FC| ≥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).

Figure 4 Analysis of DEG gene clusters based on cuprotosis-related DEGs. (A) Venn diagram showing 334 DEGs between the two CRG subtypes. (B,C) GO and KEGG enrichment analyses of DEGs between two CRG subtypes. (D) Consensus matrix heatmap of the three DEG gene clusters. (E) Kaplan-Meier curves for OS of the three DEG gene clusters. (F) Differences in the expression of 15 CRGs among the three DEG gene clusters. *, P<0.05; **, P<0.01; ***, P<0.001. BP, biological process; CC, cellular component; CRG, cuprotosis-related gene; DEG, differentially expressed gene; ECM, extracellular matrix; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; OS, overall survival.

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).

Figure 5 Construction of a prognostic risk assessment model. (A,B) LASSO regression analysis and Cox regression analysis of prognostic related genes. (C-E) Kaplan-Meier curves for OS of the three groups in all sample, training set, and testing set. (F) Differences in CRPRG scores between CRG subtypes. (G) Differences in CRPRG scores among DEG gene clusters. (H) Alluvial diagram illustrating the distribution of subtypes across different CRPRG score categories and survival outcomes. (I) Differential expression of CRGs at different high- and low-risk groups. **, P<0.01; ***, P<0.001. CRG, cuprotosis-related gene; CRPRG, cuprotosis-related pattern-related gene; DEG, differentially expressed gene; LASSO, least absolute shrinkage and selection operator; OS, overall survival.

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).

Figure 6 Development and validation of a nomogram. (A,B) Univariate and multivariate Cox regression analysis of the clinicopathological features and CRPRG score in BC. (C) Nomogram for predicting the 1-, 3-, and 5-year OS of BC patients in the all sample. (D) Calibration curves of the nomogram for predicting of 1-, 3-, and 5-year OS. (E-G) ROC curves at 1, 3, and 5 years were predicted for all sample, training set, testing set. *, P<0.05; ***, P<0.001. AUC, area under the curve; BC, breast cancer; CI, confidence interval; CRPRG, cuprotosis-related pattern-related gene; HR, hazard ratio; N, node; OS, overall survival; ROC, receiver operating characteristic; T, tumor.

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).

Figure 7 Relationship between CRPRG score and immune cells. Heat map of correlation between feature genes and immune cells, positive correlation in red and negative correlation in blue. CRPRG, cuprotosis-related pattern-related gene; NK, natural killer.

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.

Figure 8 Relationship between CRPRG score and immune cells, mutations, CSC index. (A) TMB in different risk groups. (B,C) The waterfall plot of somatic mutation features established with high- and low-risk groups. (D) Relationships between CRPRG score and CSC index. CRPRG, cuprotosis-related pattern-related gene; CSC, cancer stem cell; TMB, tumor mutational burden.

Pearson analysis revealed a positive correlation between CRPRG score and CSC index (R=0.13, P=1.5×10⁻⁵). 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 IC₅₀ of epothilone-B, ABT-263, AKT-inhibitor-VIII, and vinorelbine were lower in the high-risk group, whereas the IC₅₀ 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).

Figure 9 Validation of the protein expression of the prognostic genes based on the HPA database. (A) ADIRF, normal sample (image available from https://www.proteinatlas.org/ENSG00000148671-ADIRF/tissue/breast#img). (B) ADIRF, tumor sample (image available from https://www.proteinatlas.org/ENSG00000148671-ADIRF/cancer/breast+cancer#img). (C) CD24, normal sample (image available from https://www.proteinatlas.org/ENSG00000272398-CD24/tissue/breast#img). (D) CD24, tumor sample (image available from https://www.proteinatlas.org/ENSG00000272398-CD24/cancer/breast+cancer#img). (E) CRISP3, normal sample (image available from https://www.proteinatlas.org/ENSG00000096006-CRISP3/tissue/breast#img). (F) CRISP3, tumor sample (image available from https://www.proteinatlas.org/ENSG00000096006-CRISP3/cancer/breast+cancer#img). (G) FKBP5, normal sample (image available from https://www.proteinatlas.org/ENSG00000096060-FKBP5/tissue/breast#img). (H) FKBP5, tumor sample (image available from https://www.proteinatlas.org/ENSG00000096060-FKBP5/cancer/breast+cancer#img). (I) GPR107, normal sample (image available from https://www.proteinatlas.org/ENSG00000148358-GPR107/tissue/breast#img). (J) GPR107, tumor sample (image available from https://www.proteinatlas.org/ENSG00000148358-GPR107/cancer/breast+cancer#img). (K) MTFR1, normal sample (image available from https://www.proteinatlas.org/ENSG00000066855-MTFR1/tissue/breast#img). (L) MTFR1, tumor sample (image available from https://www.proteinatlas.org/ENSG00000066855-MTFR1/cancer/breast+cancer#img). (M) PANX1, normal sample (image available from https://www.proteinatlas.org/ENSG00000110218-PANX1/tissue/breast#img). (N) PANX1, tumor sample (image available from https://www.proteinatlas.org/ENSG00000110218-PANX1/cancer/breast+cancer#img). HPA, Human Protein Atlas.

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.

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-aw-2403/rc

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

Funding: This study was supported by Research Project Fund of The Fourth People’s Hospital of Nanning (No. NNSY2023003).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2403/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. The 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. 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]
2. Tsvetkov P, Coy S, Petrova B, et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science 2022;375:1254-61. [Crossref] [PubMed]
3. Lu K, Wijaya CS, Yao Q, et al. Cuproplasia and cuproptosis, two sides of the coin. Cancer Commun (Lond) 2025;45:505-24. [Crossref] [PubMed]
4. Basu S, Singh MK, Singh TB, et al. Heavy and trace metals in carcinoma of the gallbladder. World J Surg 2013;37:2641-6. [Crossref] [PubMed]
5. Lener MR, Scott RJ, Wiechowska-Kozłowska A, et al. Serum Concentrations of Selenium and Copper in Patients Diagnosed with Pancreatic Cancer. Cancer Res Treat 2016;48:1056-64. [Crossref] [PubMed]
6. Zhang T, Fu Y. Copper Metabolism-Related Genes as Biomarkers in Colon Adenoma and Cancer. Int J Gen Med 2025;18:3021-43. [Crossref] [PubMed]
7. Gao D, Zhou L, Bao Y, et al. Novel cuprotosis-related gene signature: a prognostic indicator and regulator of the glioma immune microenvironment. Transl Cancer Res 2024;13:6282-97. [Crossref] [PubMed]
8. Tang J, Zhao L, Qiu X, et al. Copper ionophore enhanced cisplatin efficiency through DLAT-cuprotosis. Toxicol Appl Pharmacol 2025;503:117496. [Crossref] [PubMed]
9. Zhu Z, Zhu K, Zhang J, et al. Elucidating the evolving role of cuproptosis in breast cancer progression. Int J Biol Sci 2024;20:4872-87. [Crossref] [PubMed]
10. Angell H, Galon J. From the immune contexture to the Immunoscore: the role of prognostic and predictive immune markers in cancer. Curr Opin Immunol 2013;25:261-7. [Crossref] [PubMed]
11. Gentles AJ, Newman AM, Liu CL, et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat Med 2015;21:938-45. [Crossref] [PubMed]
12. Bruni D, Angell HK, Galon J. The immune contexture and Immunoscore in cancer prognosis and therapeutic efficacy. Nat Rev Cancer 2020;20:662-80. [Crossref] [PubMed]
13. Wu P, Cai J, Fan S, et al. A novel risk score predicts prognosis in melanoma: The combination of three tumor-infiltrating immune cells and four immune-related genes. Clin Immunol 2021;228:108751. [Crossref] [PubMed]
14. Qiu X, Chen D, Huang S, et al. Identification and verification of m6A-related miRNAs correlated with prognosis and immune microenvironment in colorectal cancer. Medicine (Baltimore) 2023;102:e35984. [Crossref] [PubMed]
15. Xiao J, Liu Z, Wang J, et al. Identification of cuprotosis-mediated subtypes, the development of a prognosis model, and influence immune microenvironment in hepatocellular carcinoma. Front Oncol 2022;12:941211. [Crossref] [PubMed]
16. Liu R, Shen Y, Hu J, et al. Comprehensive Analysis of m6A RNA Methylation Regulators in the Prognosis and Immune Microenvironment of Multiple Myeloma. Front Oncol 2021;11:731957. [Crossref] [PubMed]
17. Shu S, Li CM, You YL, et al. Electroacupuncture Ameliorates Cerebral Ischemia-Reperfusion Injury by Regulation of Autophagy and Apoptosis. Evid Based Complement Alternat Med 2016;2016:7297425. [Crossref] [PubMed]
18. Wang Z, Yu J, Wu J, et al. Scutellarin protects cardiomyocyte ischemia-reperfusion injury by reducing apoptosis and oxidative stress. Life Sci 2016;157:200-7. [Crossref] [PubMed]
19. Liu S, Ge J, Chu Y, et al. Identification of hub cuproptosis related genes and immune cell infiltration characteristics in periodontitis. Front Immunol 2023;14:1164667. [Crossref] [PubMed]
20. Zhang G, Wang N, Ma S, et al. Comprehensive analysis of the effects of the cuprotosis-associated gene SLC31A1 on patient prognosis and tumor microenvironment in human cancer. Transl Cancer Res 2024;13:714-37. [Crossref] [PubMed]
21. Yue Q, Zhang M, Jiang W, et al. Prognostic value of FDX1, the cuprotosis key gene, and its prediction models across imaging modalities and histology. BMC Cancer 2024;24:1381. [Crossref] [PubMed]
22. Bao L, Locovei S, Dahl G. Pannexin membrane channels are mechanosensitive conduits for ATP. FEBS Lett 2004;572:65-8. [Crossref] [PubMed]
23. Chen W, Li B, Jia F, et al. High PANX1 Expression Leads to Neutrophil Recruitment and the Formation of a High Adenosine Immunosuppressive Tumor Microenvironment in Basal-like Breast Cancer. Cancers (Basel) 2022;14:3369. [Crossref] [PubMed]
24. Jalaleddine N, El-Hajjar L, Dakik H, et al. Pannexin1 Is Associated with Enhanced Epithelial-To-Mesenchymal Transition in Human Patient Breast Cancer Tissues and in Breast Cancer Cell Lines. Cancers (Basel) 2019;11:1967. [Crossref] [PubMed]
25. Edgar AJ. Human GPR107 and murine Gpr108 are members of the LUSTR family of proteins found in both plants and animals, having similar topology to G-protein coupled receptors. DNA Seq 2007;18:235-41. [Crossref] [PubMed]
26. Sáez-Martínez P, Jiménez-Vacas JM, León-González AJ, et al. Unleashing the Diagnostic, Prognostic and Therapeutic Potential of the Neuronostatin/GPR107 System in Prostate Cancer. J Clin Med 2020;9:1703. [Crossref] [PubMed]
27. Tonachini L, Monticone M, Puri C, et al. Chondrocyte protein with a poly-proline region (CHPPR) is a novel mitochondrial protein and promotes mitochondrial fission. J Cell Physiol 2004;201:470-82. [Crossref] [PubMed]
28. Zhang N, Dong L, Liu S, et al. MTFR1 phosphorylation-activated adaptive mitochondrial fusion is essential for colon cancer cell survival during glucose deprivation. Neoplasia 2025;63:101159. [Crossref] [PubMed]
29. Barkal AA, Brewer RE, Markovic M, et al. CD24 signalling through macrophage Siglec-10 is a target for cancer immunotherapy. Nature 2019;572:392-6. [Crossref] [PubMed]
30. Ren X, Chen X, Zhang X, et al. Immune Microenvironment and Response in Prostate Cancer Using Large Population Cohorts. Front Immunol 2021;12:686809. [Crossref] [PubMed]
31. Tong F, Lu G, Zang J, et al. FKBP5 associated CD8 T cell infiltration is a novel prognostic biomarker in luminal B breast cancer. J Int Med Res 2023;51:3000605231211771. [Crossref] [PubMed]
32. Balasubramaniam AA. Neuropeptide Y family of hormones: receptor subtypes and antagonists. Peptides 1997;18:445-57. [Crossref] [PubMed]
33. Bhat R, Thangavel H, Abdulkareem NM, et al. NPY1R exerts inhibitory action on estradiol-stimulated growth and predicts endocrine sensitivity and better survival in ER-positive breast cancer. Sci Rep 2022;12:1972. [Crossref] [PubMed]
34. Meng J, Wang LH, Zou CL, et al. C10orf116 Gene Copy Number Loss in Prostate Cancer: Clinicopathological Correlations and Prognostic Significance. Med Sci Monit 2017;23:5176-83. [Crossref] [PubMed]
35. Medeiros B, Allan AL. Molecular Mechanisms of Breast Cancer Metastasis to the Lung: Clinical and Experimental Perspectives. Int J Mol Sci 2019;20:2272. [Crossref] [PubMed]
36. Dvir K, Giordano S, Leone JP. Immunotherapy in Breast Cancer. Int J Mol Sci 2024;25:7517. [Crossref] [PubMed]
37. Mantovani A, Sozzani S, Locati M, et al. Macrophage polarization: tumor-associated macrophages as a paradigm for polarized M2 mononuclear phagocytes. Trends Immunol 2002;23:549-55. [Crossref] [PubMed]
38. Wang W, Chen S, Duan Y, et al. M2 macrophages inhibit CD8(+) T cells and facilitate the malignant biological behavior of esophageal cancer cells. Xi Bao Yu Fen Zi Mian Yi Xue Za Zhi 2024;40:887-93.
39. Li W, Wu F, Zhao S, et al. Correlation between PD-1/PD-L1 expression and polarization in tumor-associated macrophages: A key player in tumor immunotherapy. Cytokine Growth Factor Rev 2022;67:49-57. [Crossref] [PubMed]
40. Lin Y, Choukrani G, Dubbel L, et al. VISTA drives macrophages towards a pro-tumoral phenotype that promotes cancer cell phagocytosis yet down-regulates T cell responses. Exp Hematol Oncol 2024;13:35. [Crossref] [PubMed]
41. Lu Y, Han G, Zhang Y, et al. M2 macrophage-secreted exosomes promote metastasis and increase vascular permeability in hepatocellular carcinoma. Cell Commun Signal 2023;21:299. [Crossref] [PubMed]
42. Katoh M. FGFR inhibitors: Effects on cancer cells, tumor microenvironment and whole-body homeostasis Int J Mol Med 2016;38:3-15. (Review). [Crossref] [PubMed]
43. Shapouri-Moghaddam A, Mohammadian S, Vazini H, et al. Macrophage plasticity, polarization, and function in health and disease. J Cell Physiol 2018;233:6425-40. [Crossref] [PubMed]
44. Watanabe T, Ishino T, Ueda Y, et al. Activated CTLA-4-independent immunosuppression of Treg cells disturbs CTLA-4 blockade-mediated antitumor immunity. Cancer Sci 2023;114:1859-70. [Crossref] [PubMed]