Based on data from The Cancer Genome Atlas (https://portal.gdc.cancer.gov/), information on RNA-sequencing raw data of 1110 cancerous breast samples as well as 112 normal human breast samples that included therapeutically information, somatic mutation data and CNV data files, was obtained. It was necessary to download processed gene expression datasets, clinical samples collected from breast cancer patients (n=58), as well as normal breast tissue (n=4) from the Gene Expression Omnibus profile database (https://www.ncbi.nlm.nih.gov/geo/) (ID: GSE61304). These raw data were first standardized to fragments per kilobase million expression levels prior to comparison and figuring out the expression of CRGs. After that, CRG expression was determined using the limma program (20). We integrated the data once the data cleaning procedure was finished to get them ready for analysis. The study that followed did not include patients for whom there was inadequate data on their survival.

19 CRGs made up the signature that we were able to collect from earlier publications (9, 21–25), the list of genes is in Table S1 . We were able to classify individuals into discrete molecular clusters based on their CRG expression using the ConsensusClusterPlus R program (26). Through the use of unsupervised clustering, this was done. The clinical usefulness of CRGs in breast cancer was investigated using the Kaplan-Meier approach in a Kaplan-Meier study. We used the survival and survminer packages in R to examine the curves of survival as well as display the results. After that, the ggplot2 software was used to do a principal component analysis. The two subtypes' biological processes were maintained by using the Gene Set Variation Analysis tool (27). Malignant Tumour tissues employing expression  (28) and CIBERSORT (29) were also utilized to represent the percentage of immune and stromal cells in patients with breast cancer. The extent to which each immune cell within each sample carried an enrichment score was also assessed using an analysis of gene set enrichment on a single sample (30).

We associated the two cuproptosis-related subtypes with the primary clinical and pathological parameters of breast cancer patients, including their age, T phase, and N phase, as well as their prognosis, as part of our inquiry into the possible clinical functions of the two cuproptosis-related subtypes. Additionally, Kaplan-Meier survival analysis technique was utilized to look at differences in overall survival that were verified amongst the various subtypes. We discovered the differentially expressed genes between the cuproptosis-related subgroup using the limma R program. These genes required to possess an adjusted p-value < 0.05 and a fold change > 1.5. To clarify the pathways that were considerably enriched, gene ontology enrichment analysis and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis were also conducted. To further explore the hidden roles among the DEGs, the data were displayed using the ClusterProfiler program (31).

By computing the overall value of risk, the value of CRG was established in order to identify the cuproptosis patterns in specific patients. For a more thorough analysis, we utilized unsupervised consensus clustering to separate the breast cancer patients into four different subtype groups (cuproptosis-related gene subtype A-D). The train sets were then utilized to generate a CRG score for prognosis. A percentage of 1:1 was applied to all patient datasets in order to divide them into train and test sets. The glmnet package in R was used to perform least-squares regressions and selection operator regressions in order to minimize the possibility of overfitting the model (32). For the purposes of predicting the OS of the patients in the training set, a multivariate Cox regression with proportional hazards analysis was also utilized. Both the train set and the test set were split into groups of high-risk and low-risk based on their risk ratings. In each set, Kaplan-Meier analyses of survival and ROC curves were conducted.

The relationships between the CRG score and the clinically significant parameters, including age, T and N stage, were examined using chi-square tests. On both the train and test sets, univariable and multivariable analyses were conducted to see if the CRG score was influenced by any other easily accessible clinicopathological characteristics. Age, tumor grade, T and N stage were also taken into account in subgroup studies to see whether the CRG score still had the same predictive value it did in our model earlier.

A nomogram was created using the rms program to predict overall survival based on clinically significant characteristics and the CRG score. Each clinicopathologically significant characteristic was given a score using the nomogram model, and the overall score was obtained by summing all the individual scores. By contrasting the area under the time-dependent ROC curves of survival rates after one, three, and five years, the nomogram's accuracy in predicting survival rates was also validated. Additionally, model calibration was performed to compare the predicted likelihood of survival outcomes across the 1-, 3-, and 5-year periods with the actual survival occurrences.

To calculate the total number of tumor-infiltrating immune cells and subgroups of immune cells in each sample, we utilized the CIBERSORT method for comparison of 23 immune cells infiltrating the tumors between the high-risk and low-risk groups. This was done to assess how many immune cells altogether had invaded the tumor. The gene groups connected with 23 levels of immune cell infiltration were also found using the CRG score. Additionally, we looked at the connections between the CRG risk score and the cancer stem cell.

Depending on whether a sample was deemed high-risk or low-risk, its tumor mutation burden was assessed for each. Additionally, we used the maftools program to do a somatic variant analysis on patients with breast cancer in order to look at and analyse the somatic mutation data (33). Using the pRRophetic software, we calculated the semi-inhibitory concentrations (IC50) of frequently prescribed medications in breast cancer patients depending on their risk levels (34).

The R-based statistical analysis was conducted with a significance threshold of p 0.05 (version 4.0.2). 

On the TCGA dataset, 1110 breast cancer patients' expression levels of 19 cuproptosis-related genes were examined, along with 112 normal human breast tissues ( Figure 1A ). In the meanwhile, gene mutation analyses revealed that 55 out of the 976 samples (5.64%) had CRGs mutations, with ATP7A having the greatest gene mutation rates ( Figure 1B ). The majority of CRGs were accumulated on copy number loss or deletion, according to an examination of copy number variations ( Figure 1C ), and all 19 CRGs had frequent copy number alterations ( Figure 1D ). Additionally, it was discovered via research of the impact of gene expression patterns on overall survival in breast cancer that those expressing high levels of ATP7A, DBT, DLAT, DLD, GLS, PDHA1, and SLC31A1 had a bad prognosis. A higher level of ATP7B, LIPT1, and NLRP3 expression is linked to improved OS ( Figure 1E–N , and Table S2 ). The findings suggested that CNV alterations could modify the way CRGs are expressed. Additionally, a relationship between CRG expression levels and breast cancer prognostic variables was discovered, pointing to a potential involvement for CRGs in breast cancer. The biomarkers might be used as therapeutic targets or prediction biomarkers for breast cancer.

The correlation network picture showed the 19 CRGs' strong association with one another ( Figure 2A ). The cohort was subdivided into two groups, group A (n = 534) and group B (n = 605), based on a consensus cluster analysis of the 1168 breast cancer samples, which showed that a cluster of κ = 2 had the largest intragroup and lowest intergroup differences ( Figure 2B , and Figure S1 ). Differences in the transcription patterns of the two subtypes of cuproptosis were found using PCA ( Figure 2C ). Subtype B has a better prognosis than subtype A, according to Kaplan-Meier survival calculations (p = 0.001; Figure 2D ). A heatmap was created as a consequence of the relationship between features of clinical significance and patterns of CRG expression ( Figure 2E ). The bulk of CRGs expressed themselves more strongly in subtype A, whereas late-phase breast cancer was represented in subtype B.

We used gene set variation analysis enrichment analysis to look at the two subtypes' possible effects on biological behavior ( Figure 3A ). Compared to subtype B, subtype A had an enrichment in the pathways linked to immunological activation. According to a GSVA enrichment study, subtype A is considerably enriched in metabolic-activated pathways, such as the folate utilization of one carbon pool, lysine degradation, the citrate cycle, RNA metabolism, arachidonic acid metabolism, and N-glycan biosynthesis. In each breast cancer sample, we used the CIBERSORT method to assess the associations between two subtypes as well as the 23 other subtypes of immune cells in order to learn more about how CRGs work in the tumor microenvironment. According to our research, there are significant variations between the two subtypes in the quantity of immune cells that infiltrate ( Figure 3B ). As compared to subtype B, CD4 T cells, type 2 T helper cells, regulation T cells, gamma delta T cells, immature dendritic cells, and immature B cells were found to be more prevalent in subtype A. Subtype A, on the other hand, exhibited considerably reduced levels of neutrophil, eosinophil, mast cell, and CD56 dim natural killer cell infiltration. Then, we did a functional enrichment analysis to look into the two cuproptosis subtypes' possible biological roles after using the limma algorithm to identify 591 DEGs linked to them ( Figure S2 and Table S3 ). It was discovered that CRGs were mainly engaged in membrane protein targeting, membrane protein localization establishment, and pathway analysis using GO and KEGG ( Figures 3C, D , and Figure S3 ).

Using a consensus clustering technique, 1139 breast cancer patients were categorized into four molecular genetic categories based on prognostic genes. Subtypes A (n = 350), B (n = 502), C (n = 165), and D (n = 122) were found when κ = 4 indicated that the breast cancer instances may be separated into four subclasses ( Figure 4A and Figure S4 ). Additionally, the relationship between the clinical traits of breast cancer patients and the gene subtypes was investigated ( Figure 4B ). The genetic subtype D patients had the lowest OS, while patients with genetic cluster C had the greatest OS, according to Kaplan-Meier curves (p < 0.001; Figure 4C ). The four cuproptosis gene subtypes' expression of CRGs varied greatly, as expected by the cuproptosis patterns ( Figure 4D ).

Based on DEGs related to subtypes, a LASSO-Cox regression model was developed to provide a predictive CRG score for each patient. Figure 5A illustrates the proportion of patients among the two CRG score groups, the two cuproptosis subtypes, and the four gene subtypes. There was a statistically significant variation in CRG scores across cuproptosis subtypes. Subtype B had a much higher CRG score than subtype A. Figure 5B displays the risk score distributions for the two CRG subtypes. The highest CRG scores were for subtype D, while the lowest were for subtype C ( Figure 5C ).

Then, using R's caret package, patients were randomly assigned to training groups (n = 570) as well as testing groups (n = 569) at a ratio of 1:1 ( Tables S4, S5 ). Using LASSO and multivariate Cox analysis, 22 OS-related genes were selected using the least partial likelihood deviation from 591 cuproptosis subtype-related prognostic DEGs ( Figures 5D, E , and Table S6 ). Based on a Cox regression analysis involving several variables, Akaike information criteria value of 22 OS-associated genes was utilized to identify six genes (PGK1, RPL14, PRDX1, PSME1, MAL2, and SURF4) ( Table S7 ). These findings led to the following formula being chosen as the risk score formula: The risk score is calculated as follows: (0.00375* PGK1 expression) + (-0.00930*RPL14 expression) + (0.00278*PRDX1 expression) + (-0.00668*PSME1 expression) + (0.00147*MAL2 expression) + (0.00672*SURF4 expression). 13 out of 17 hallmark genes showed a significant variation in their expression of genes between high-risk individuals and low-risk individuals ( Figure 5F ). Based on their risk ratings, each theme was split into high- and low-risk patient groups, and the median scores were calculated for the training and test sets. According to their values of risk, patients were split into two groups: those at low risk and those at high risk ( Figure 6 ). In terms of survival rates and circumstances, there were significant differences among the two groups based on Kaplan-Meier curves. Patient survival rates and the distribution of CRG scores were analyzed independently for the train and test sets.

Using the data gathered, we created a nomogram using the rms program to forecast the life expectancy of breast cancer patients at the lifetime of 1, 3, and 5 ( Figure 7A , and Table S8 ). Each patient's total point values were determined based on prognostic characteristics such as their age, level of risk (low risk was indicated by a "low CRG score" and high risk was indicated by a "high CRG score"), as well as the T and N stage of their ailment. The harshness of the prognosis is directly correlated with the patient's overall score. The calibration plots showed that the nomogram performed better than an ideal model would have ( Figure 7B ). Additionally, ROC analysis indicated that the nomogram performed very well in terms of prediction ( Figures 7C–E ).

The CIBERSORT algorithm was used to assess the relationship between the CRG score ( Figure S5 ) and the number of immune cells. However, the CRG score was negatively correlated with naive B cells, resting dendritic cells, resting mast cells, monocytes, activated NK cells, plasma cells, CD8 + T cells, and follicular helper T cells. A correlation was found between the CRG score and activated memory CD4 + T cells, M0 macrophages, M2 macrophages, activated mast cells, and resting NK cells ( Figure 8A ). Additionally, our research looked at the association between six genes and the amount of immune cells. According to our study, the six genes affect the bulk of immune cells ( Figure 8B , and Table S9 ).

The TMB study revealed a significant association between anticipated TMB level and cuproptosis gene subtypes (R = 0.28, P < 0.001; Figure 8C ). To give further support, we looked at the variations in somatic mutation distribution across the cohort's two CRG score groups. The top 10 most changed genes in each of the two groups were PIK3CA, TP53, TTN, CDH1, GATA3, MUC16, MAP3K1, HMCN1, and FLG. The most often mutated genes in patients with a high CRG score are TP53 (46%) and PIK3CA (28%), while PIK3CA (41%) is the most frequently mutated gene in the low-risk category ( Figures 8D, E ).

Additionally, it was shown that there was a link between the CRG score and the CSC index that was positive (R = 0.22, P < 0.001), suggesting that cells from breast cancer with higher cell retention gene scores demonstrated more stem cell features and less differentiation ( Figure 8F ). Sensitivity analysis was done on a few medications presently being used to treat breast cancer among the two groups. For patients with high CRG scores, it was found that the IC50 values of drugs including paclitaxel, vinblastine, bleomycin, AUY922, ATRA, and AZD6244, among others, were considerably higher. It is evident from these results that CRGs are essential for the sensitivity of drugs ( Figures 9A–F ).