The necessary data were acquired from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/), comprising of RNA-sequencing (RNA-seq) data, clinical information, and simple nucleotide variation data from a total of 1182 BRCA patients. Fragments per kilobase million (FPKM) was selected for mRNA expression profiling. We also acquired RNA-seq data and clinical information from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/, ID: GSE20685, GSE42568, GSE58812). Specifically, GSE20685 contained RNA-seq and clinical data of 327 breast cancer samples, GSE42568 contained RNA-seq and clinical data of 104 breast cancer samples and 17 normal breast samples, and GSE58812 contained RNA-seq and clinical data of 107 breast cancer samples. To ensure uniformity, the mRNA expression data in the TCGA and GEO databases were transformed into transcripts per million (TPM) data by log2 transformation. We downloaded a total of 65 ferroptosis-related genes from the Molecular Signatures Database (MSigDB, http://www.gsea-msigdb.org/gsea/msigdb/) and an additional 19 cuproptosis-related genes from the literature (Tsvetkov et al., 2022) (Ren et al., 2021) (Tao et al., 2021) (Selvaraj et al., 2005).

GSVA is an algorithm used to analyze gene expression data and functional annotation. It can display the differences in samples on different gene sets by calculating gene set scores. The algorithm does not require standardization and is very fast, so it is widely used in biomedical research. A GSVA of the GSE42568 samples was performed on ferroptosis and cuproptosis genomes to obtain ferroptosis and cuproptosis scores, respectively. The hub genes were screened by weighted gene coexpression network analysis (WGCNA) based on their scores. Scale-free and average connectivity analyses were performed on modules with different power values using the PickSoftThreshold function to set the soft threshold power to 5. Then, the corresponding dissimilarity matrix (1-TOM) and topological overlap matrix (TOM) were obtained. Pearson correlation analysis was performed on coexpression modules with cuproptosis and ferroptosis scores after TOM.

X-tile software is a software designed by Yale University to calculate the optimal cutoff value for survival curves. We calculated the optimal cutoff values for the cuproptosis score and ferroptosis score in the TCGA cohort with this software and divided the scores into high- and low-score groups, respectively. Subgroups with both high and low scores for both cuproptosis and ferroptosis were selected. According to specific criteria (|log2FC| ≥ 1 and FDR<0.05), the “DeSeq2” R package was used to screen out the high and low subgroup differentially expressed genes (DEGs), and then the DEGs were cross-referenced with the genes screened by WGCNA, and the final genes were identified as CFRGs.

After screening out CFRGs related to prognosis, we used five machine learning methods—decision tree, random forest, least absolute shrinkage and selection operator (LASSO), extreme gradient boosting (XGBoost) and gradient boosting decision tree (GBDT)—to evaluate the survival weights of the CFRGs. GSE42568 was used as the training set, and the first 30 genes were selected to construct the LASSO Cox model. With the median risk score as the critical value, the model was divided into a high-risk group and a low-risk group. The overall survival (OS) of the high and low risk groups was analyzed using the KM curve. Subsequently, Kaplan‒Meier (KM) analysis was performed on the high- and low-risk groups using the “survival” and “survminer” R packages, and 1-, 3-, and 5-year receiver operating characteristic (ROC) analyses were performed with the “timeROC” R package to evaluate the sensitivity of the model. Finally, the TCGA and GEO cohorts (GSE20685, GSE58812) were used as the validation set, and the obtained prognoses were selected for validation. In addition, we constructed a nomogram by combining the risk score and clinical data.

We used seven algorithms from the Tumor Immune Estimation Resource (TIMER) database, including TIMER, CIBERSORT, CIBERSORT-ASB, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC, to analyze signatures and immunity. In addition, we evaluated the differences in immune cell infiltration and immune function between high-risk and low-risk subgroups using the single-sample gene set enrichment analysis (ssGSEA) algorithm. Furthermore, differences between high- and low-risk groups in immune therapy biomarkers were also analyzed through TIDE (http://tide.dfci.harvard.edu/)), TMB, and immune checkpoint.

Single-cell RNA-seq data for 12 samples were obtained from the GSE149655 dataset in the GEO database. The sequencing data were analyzed using the “Seurat” R package, and low-viability cells were removed. Using the CreateSeuratObject function, and based on specified criteria, low-quality cells were further removed (min.cells<3, min. features<50, percent. mt < 50). The filtered single-cell RNA sequencing (scRNA-seq) data were then normalized, and 1,500 genes with large coefficients of variation between cells were generated and subsequently subjected to principal component analysis (PCA) for dimensionality reduction of the scRNA-seq data. The top 4 PCs were selected for cell clustering analysis using the distributed stochastic neighbor embedding (t-SNE) algorithm. Through the DEGs of each cluster, the cell clusters were tool annotated. Finally, we analyzed the expression levels of OCFRGs in different clusters.

DEGs in high-risk and low-risk subgroups were screened and GSEA was performed using the GSEA software (https://www.gsea-msigdb.org/gsea/index.jsp). In addition, GSVA was used to perform GSVA on high-risk and low-risk subgroups using the “GSVA” R package. We collected some pathways closely related to the occurrence and development of tumors from the literature, including Hippo, Wnt, MAPK, PI3K/AKT, TGF-β, NF-kB, Notch, AMPK, JAK-STAT, PD-1/PD-L1, mTOR, Ras, TNF, HIF-1, ErbB, Nrf2, as well as functions and pathways related to cuproptosis and ferroptosis, including GLUTATHIONE PEROXIDASE ACTIVITY, P53 SIGNALING PATHWAY, GLUTATHIONE METABOLIC PROCESS, LIPID HOMEOSTASIS, PENTOSE PHOSPHATE PATHWAY, TRICARBOXYLIC ACID CYCLE, MITOCHONDRIAL TRICARBOXYLIC ACID CYCLE ENZYME COMPLEX, TRICARBOXYLIC ACID CYCLE ENZYME COMPLEX, REGULATION OF OXIDATIVE STRESS INDUCED CELL DEATH, RESPONSE TO OXIDATIVE STRESS, CELL MIGRATION INVOLVED IN SPROUTING ANGIOGENESIS, IRON METABOLISM. We used GSVA to calculate the enrichment score for each pathway and evaluate the correlation between risk score and sex enrichment score.

In order to assess the sensitivity of high and low risk groups to drugs, we used the “pRRophetic” R package to analyze the RNA-seq data of BRCA based on the Genomics of Drug Sensitibity in Cancer (GDSC) database. The level of half maximal inhibitory concentration (IC50) reflects the patient’s sensitivity to drugs.

Five breast invasive ductal carcinoma tissue chips were purchased from Shanghai Outdo Biotech Company (Shanghai, China). Each tissue chip includes 45 cancer tissues and 45 paracancerous tissues. ANKRD52 (rabbit polyclonal, catalog number: GTX32443, Genetex), HOXC10 (rabbit polyclonal, catalog number: ab153904, Abcam), KNOP1 (rabbit polyclonal, catalog number: ab126512, Abcam), SGPP1 (rabbit polyclonal, catalog number: ab126512, Abcam): ab126512, Abcam) and TRIM45 (rabbit polyclonal, catalog number: TA505920, origene). The results of the immunohistochemical staining were scored. Semiquantitative scoring was performed according to the staining intensity and the percentage of positive cells: No staining, pale yellow (light yellow particles), medium (brown yellow particles), and heavy (dark brown particles) were scored as 0, 1, 2, and 3, respectively. According to the percentage of positively stained cells in the total number of cells, 0% was scored as 0, 5%–25% was scored as 1; 26%–50% was scored as 2; 51%–75% was scored as 3; and >75% was scored as 4 points. The final score was the sum of the staining intensity and the percentage of positive cells. The sum of the staining intensity and the percentage of positive cells was less than 6 for the low expression group and ≥6 for the high expression group. Five 400x high-power fields were randomly selected for each sepction, the staining intensity and percentage of positive cells were scored in each field, and the average value was calculated. The immunohistochemical staining results were microscopically adjudicated by two pathologists in an independent, double-blind manner.

Normal breast epithelial MCF-10A cells and three human breast cancer cell lines, SK-BR-3, MDA-MB-231, and MCF-7, were obtained from the Central Laboratory of Shandong Provincial Hospital. Total RNA was extracted using TRIzol reagent (Invitrogen, United States). Complementary DNA (cDNA) was synthesized using the PrimeScript RT kit (Takara). Supplementary Table S2 shows the primer sequences for Quantitative real-time PCR analysis of OCFRGs.

Cells were lysed in cold Radioimmunoprecipitation assay (RIPA) buffer. The same amount of protein was subjected to SDS-PAGE, and then transferred to PVDF (polyvinylidene fluoride) membrane. Block with nonfat dry milk containing TBST for 1 h. The primary antibody (Western blot and IHC universal primary antibody) was diluted according to the instructions and incubated overnight at 4°C. After washing with TBST, the secondary antibody was added and incubated for 1 h at room temperature. After washing the membrane, it was developed using enhanced chemiluminescence (ECL) chromogenic solution.

The bioinformatics analysis, statistical analysis, and machine learning aspects of this study were all performed with R software (version 4.0.1). Statistical tests included the Pearson chi-square test and Wilcoxon rank-sum test. Kaplan-Meier analysis was used to assess OS in both groups of patients. p < 0.05 was considered to be statistically significant.

The specific process of this study is shown in Figure 1, which includes the acquisition of data, the acquisition of cuproptosis and ferroptosis-related genes, the construction and validation of the signature, the analysis of clinical data, the analysis of immune cell infiltration levels, and the analysis of pathways and functions.

We used the pickSoftThreshold function in the WGCNA R package to automatically select a soft threshold of 5 (Figure 2A). Multiple gene modules were divided by the dynamic cutting method, and the mergeCloseModules function was used to perform cluster analysis on each module into fewer modules and mark them with different colors (Figure 2B). To determine the correlation of coexpression modules with genetic differences between normal and tumor samples that may contribute to tumor development, cuprotosis score and ferroptosis score, we used Pearson correlation analysis, and the results showed that the module “turquoise” had the strongest correlation (R (Turmor) = 0.86, R (cuproptosis) = 0.77, R (ferroptosis) = 0.57) (Figure 2C).

X-tile software was used to determine the optimal cutoff values for the survival curves of the cuproptosis and ferroptosis scores in the TCGA BRCA cohort (Supplementary Figure S1). Patients with high cuproptosis and ferroptosis scores and patients with low cuproptosis and ferroptosis scores were divided into two groups, and DEGs were screened out with |log2FC| ≥ 1 and FDR<0.05 as criteria (Figure 2D). By analyzing the intersection of the genes screened by means of DEGs and WGCNA, a total of 541 CFRGs were identified (Figure 2E).

Survival-related CFRGs were screened by using Cox regression analysis (Figure 3A). The survival-related weights of the CFRGs were evaluated using five machine learning algorithms, including decision tree, random forest, LASSO, XGBoost, and GBDT. The top 30 CFRGs with average scores (Supplementary Table S1) were selected to screen genes and to construct models by means of LASSO Cox regression analysis. According to the best λ value, 5 prognostic genes (Figure 3B) were selected and signature was constructed (Figure 3C,D). Calculation of the risk score was conducted as follows: R i s k   S c o r e = 0.4024 * A N K R D 52   exp ⁡ . + 0.1025 * H O X C 10 ⁡ exp ⁡ . + 0.2146 * K N O P 1 ⁡ exp ⁡ . + − 0.3069 * S G P P 1   exp ⁡ . + − 0.4003 * T R I M 45   exp ⁡ .

Based on the median training set risk score, BRCA patients were divided into high- and low-risk groups and further analyzed (Figure 3F). In the training and validation sets, the KM curve showed that the OS (Overall survival) in the high-risk group was lower than that in the low-risk group (Figure 4A), which was statistically significant (p < 0.05). The ROC curve of the training set showed that the AUCs at 1 year, 3 years, and 5 years were 0.76, 0.726, and 0.728, respectively. The AUC of the TCGA cohort was 0.734, 0.713, and 0.756, and that of GSE58812 was 0.784, 0.714, and 0.698. This indicates that our model has reliable predictive performance. We constructed a nomogram, in which Age (p < 0.001), Stage (p < 0.001), and risk score (p < 0.01) were all statistically significant (Figure 4C). In the nomogram, the OS of high-scoring patients was significantly lower than that of low-scoring patients (Figure 4D), which was statistically significant (p < 0.05). The column chart had good predictive performance with AUCs of 0.876, 0.773, and 0.727 at 1, 3, and 5 years, respectively.

We analyzed the relationship between OCFRGs (ANKRD52, HOXC10, KNOP1, SGPP1, TRIM45) and clinical data. The results showed that the expression level of HOXC10 differed significantly in different T and N stages and was statistically significant. In addition, the higher the expression level of HOXC10, the later the T stage and N stage of the patient (see Figures 5A,B). The expression levels of KNOP1 and SGPP1 also differed significantly in different T stages (see Figures 5C,D) and were statistically significant. We also analyzed the correlation between the expression levels of OCFRGs, risk score, and the expression levels of ER, PR, and HER2 (see Figure 5E). The results showed that the expression level of ANKRD52 was positively correlated with ER (R = 0.28). TRIM45 was also positively correlated with ER (R = 0.38) and PR (R = 0.33). There were also differences in stage, T stage, and M stage between the high-risk and low-risk groups (Figure 6A).

We analyzed the relationship between several immune-related markers and the risk subtypes. The results showed that the expression of numerous immune checkpoints was different between high and low-risk groups (Figure 6A), especially CD274 (PDL1, Programmed Cell Death 1 Ligand 1). Additionally, TMB (tumor mutational burden) and TIDE (tumor immune dysfunction and exclusion) are considered important immune therapy markers. Our study found that TMB and TIDE were also different between high and low-risk groups. The results based on the ssGSEA algorithm showed that the expression levels of Mast cells, B cells, NK cells, and Neutrophils also differed between high and low-risk groups. Moreover, according to the results from TIMER database (Figure 6B), all immune-related markers had different expression levels between the high and low-risk groups. It is worth noting that all of the immune-related analyses had statistical significance (p < 0.05).

In terms of immune cell infiltration, the expression of OCFRGs was mainly correlated with the infiltration levels of neutrophils and regulatory T cells (Tregs). The expression of ANKRD52 was negatively correlated with the infiltration level of activated dendritic cells (aDCs) (R = −1.4, p < 0.001, Figure 6C), and the expression level of HOXC10 was negatively correlated with the infiltration level of Tregs (R = −0.076, p = 0.014, Figure 6D). The expression level of KNOP1 was negatively correlated with the infiltration level of neutrophils (R = −0.22, p < 0.001, Figure 6E), and the expression level of SGPP1 was positively correlated with the infiltration level of neutrophils (R = 0.34, p < 0.001, Figure 6F). The expression level of TRIM45 was positively correlated with the infiltration level of Tregs (R = −0.34, p < 0.001, Figure 6G).

We obtained 14592 cells from 12 samples of GSE168410. A total of 1,418 cells were obtained by screening the total cells according to the intracellular gene features, the percentage of chromosomal genes, etc. Afterward, the top 10 genes with large coefficients of variation were labeled. Latitude reduction was performed using PCA (Figure 7A). The first 4 PCs were retained for further t-SNE, and 13 cell subsets were obtained (Figure 7B). The top 5 genes in each cluster were visualized with a heatmap (Supplementary Figure S3). We used tools to annotate the cell subsets, namely, fibroblasts, tissue stem cells, epithelial cells, monocytes, endothelial cells, chondrocytes, and T cells (Figure 7C). To investigate the expression of marker genes in different cells, we visualized them with t-SNE and violin plots (Figures 7D,E). OCFRGs were highly expressed in epithelial cells; HOXC10, KNOP1, and SGPP1 were highly expressed in fibroblasts; and only SGPP1 was highly expressed in monocytes.

The biological functions of high- and low-risk groups were analyzed using GSEA software, which revealed that certain functions such as CELL CYCLE, CITRATE CYCLE/TCA CYCLE, DNA REPLICATION, GLYCOLYSIS/GLUCONEOGENESIS, MISMATCH REPAIR, PENTOSE PHOSPHATE PATHWAY, and PROTEASOME were more active in the high-risk group. On the other hand, the low-risk group exhibited significant enrichment in ARACHIDONIC ACID METABOLISM, MAPK SIGNALING PATHWAY, and P53 SIGNALING PATHWAY (Figure 8A). Additionally, GSVA analysis (Figure 8B) indicated that CELL CYCLE, HOMOLOGOUS RECOMBINATION, MISMATCH REPAIR, and DNA REPLICATION were more enriched in the high-risk group, whereas INTESTINAL IMMUNE NETWORK FOR IGA PRODUCTION, CYTOKINE CYTOKINE RECEPTOR INTERACTION, TYROSINE METABOLISM, and ARACHIDONIC ACID METABOLISM evidenced more enrichment in the low-risk group. Remarkably, several signaling pathways, including mTOR, HIF-1, and ErbB, were significantly activated in the high-risk group, whereas other signaling pathways such as MAPK, PI3K/AKT, TGF-β, and Ras were significantly inhibited (Figure 8C). Furthermore, the scores for cuproptosis/ferroptosis-related functional pathways, including Nrf2 SIGNALING PATHWAY, P53 SIGNALING PATHWAY, TRICARBOXYLIC ACID CYCLE, and REGULATION OF OXIDATIVE STRESS INDUCED CELL DEATH, were positively correlated with the risk score. In contrast, the scores for GLUTATHIONE PEROXIDASE ACTIVITY, GLUTATHIONE METABOLIC PROCESS, and IRON METABOLISM were negatively correlated with the risk score. All the aforementioned analyses exhibited statistical significance (p < 0.05).

Elesclomol is a potent copper ionophore that promotes cuproptosis. Our analysis shows that the IC50 of various anticancer drugs including elesclomol is different in the high- and low-risk groups, and the IC50 values of high-risk groups are generally higher than those of low-risk groups (Figure 9A–L).

To explore the differences in the protein expression of OCFRGs in normal tissues and adjacent tissues, we detected the expression levels of ANKRD52, HOXC10, KNOP1, SGPP1, and TRIM45 in tissues by using immunohistochemical staining. The results of ANKRD52 level expressed that the HOXC10, KNOP1, and TRIM45 in the tumor tissue were higher in the paracancerous tissue, and the opposite was true for SGPP1 (Figure 10A). qRT-PCR results showed that ANKRD52, HOXC10, KNOP1, and TRIM45 were significantly higher in MCF-10A cells than other cell lines, while SGPP1 was the opposite (Figure 10B). The results of Western blotting also verified the results of the database (Figure 10C).