Bioinformatic analysis and experimental validation of six cuproptosis-associated genes as a prognostic signature of breast cancer

Data sources

The Cancer Genome Atlas (TCGA)-BRCA cohort comprising sequencing data and clinical information of 1,091 BRCA samples and 113 normal samples, was acquired from TCGA database (https://xenabrowser.net). Of these, 1,069 BRCA samples with survival information (after filtering out samples without details of age, M, N, T, and stage) were incorporated into the survival analysis. GSE42568 and GSE20711 cohorts were retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/) and employed as external validation sets for assessing the cuproptosis-relevant prognostic signature. GSE42568 and GSE20711 cohorts included 104 and 88 BRCA samples with corresponding survival information and data, respectively. The RNA-seq FPKM data of TCGA-BRCA cohort and GSE42568 dataset were normalized with log₂ (FPKM+1) values. Ten cuproptosis-relevant genes (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A) were extracted from a previous study (Tsvetkov et al., 2022).

Additionally, human epithelial cell line from the mammary gland, MCF-10A (normal group), and three breast cancer cell lines, HCC1937, MCF-7, and MDA-MB-231 (BRCA group), were acquired from iCell Bioscience Inc (Shanghai, China), they were employed to perform RT-qPCR. Specifically, in the environment of 37 °C with 5% CO₂, MCF-10A in MEGM Kit medium (Lonza/Clonetics, CC-3150), HCC1937 cells in RPMI-1640 medium (iCell-0002), MCF-7 cells in MEM basic medium (iCell-0012), and MDA-MB-231 cells in L15 medium (iCell-0009) were separately cultured.

Functional enrichment analysis

The online website DAVID (https://david.ncifcrf.gov/) (Sherman et al., 2022; Huang, Sherman & Lempicki, 2009) and the R package ‘clusterProfiler’ (Wu et al., 2021) were employed for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses (Kanehisa et al., 2017). GO terms were annotated as cellular component (CC), molecular function (MF), and biological process (BP).

Identifying BRCA-related subtypes

BRCA cases in TCGA-BRCA cohort were classified based on the expression of corresponding genes by consensus clustering using the ‘ConsensusClusterPlus’ package (Wilkerson & Hayes, 2010). The consensus matrix and consensus CDF curve were generated for selecting the optimal typing. Uniform manifold approximation and projection (UMAP) maximized the retention of features of the original data while reducing the feature dimensionality and was used to conduct dimensional reduction analysis on different subtypes.

Differential expression analysis

Depending on p-value < 0.05 and |log₂FoldChange(FC)| > 0.5, we obtained differentially expressed genes (DEGs) using the ‘limma’ package (Ritchie et al., 2015).

Tumor microenvironment analysis

The ‘estimate’ package in R was utilized to compute the immune, stromal, and ESTIMATE scores, along with tumor purity for each BRCA sample (Yoshihara et al., 2013). Using the ssGSEA (Hänzelmann, Castelo & Guinney, 2013) and CIBERSORT algorithms (Newman et al., 2015), immune gene sets and the fraction of immune infiltrating cells were calculated for each BRCA sample.

Establishing cuproptosis-relevant prognostic signature in BRCA

First, we acquired genes that were significantly associated with overall survival (OS) of BRCA patients by univariate Cox analysis. Subsequently, we calculated the two principal components, PC1 and PC2, for each sample. The PC1 and PC2 values were subjected to multivariate Cox regression analysis to obtain corresponding coefficients; according to the formula, $\mathrm{C}\mathrm{u}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\sum \left(\mathrm{P}\mathrm{C}{1}_{\mathrm{i}}+\mathrm{P}\mathrm{C}{2}_{\mathrm{i}}\right)$, where i represents the expression of genes, we calculated the Cusig score of each BRCA sample. The surv cutpoint function of the ‘survminer’ package (He et al., 2024) was used to obtain the optimal cutoff value to classify BRCA patients into the high-CuSig score and low-CuSig score groups. Kaplan-Meier (KM) curves were generated using the ‘survminer’ and the ‘survival’ packages (He et al., 2024; Liu et al., 2021; Shen et al., 2023).

Enrichment analysis of pre-defined gene sets based on the prognostic signature

Gene set variation analysis (GSVA) was implemented through the ‘gsva’ package in R (Hänzelmann, Castelo & Guinney, 2013; Subramanian et al., 2005) and ‘c2.cp.kegg.v7.2.symbols.gmt’ and ‘h.all.v7.5.1.symbols.gmt’ were the reference gene sets. Differential HALLMARK/KEGG entries were then obtained using the ‘limma’ package, and the filtering criteria were |log₂FC| > 0.1 and p-value < 0.05. Gene set enrichment analysis (GSEA) (He et al., 2024) was also conducted by setting the GO gene set as the reference. The threshold values for significant entries were SIZE > 20 and NOM p-value < 0.05.

Immunotherapy efficacy analysis based on the prognostic signature

The sensitivity of the two Cusig score subgroups to immune checkpoint blockade (ICB) therapy was inferred and assessed on the Tumor Immune Dysfunction and Exclusion (TIDE) website (Jiang et al., 2018). The calculated TIDE value was used for evaluating the efficacy of immunotherapy.

Assessing mRNA expressions of cuproptosis-relevant prognostic genes in cell lines

The RT-qPCR validation was performed and reported according MIQE guidelines (Bustin et al., 2009). Total RNA from the one normal and three BRCA cell lines in the logarithmic phase was extracted using the TRIzol Reagent following the manufacturer’s instructions (Ambion, Austin, TX, USA). Chloroform (Chengdu Guerda Rubber Industry Co., LTD, Chengdu, China) was employed to remove proteins and fat-soluble magazines, ice isopropanol (Chengdu Guerda Rubber Industry Co., LTD, Chengdu, China) was utilized to precipitate RNA, and 75% ethanol (Chengdu Colon Chemical Co., LTD, Chengdu, China) was applied to further remove impurities. Following this, the RNA was solubilized by adding 20–50 μL of RNase-free water (Servicebio, Guangzhou, China) to the obtained RNA precipitate and the RNA concentration was detected with NanoPhotometer N50. Subsequently, total RNA was reverse transcribed into cDNA using the SweScript-First-strand-cDNA-synthesis-kit (Servicebio, Guangzhou, China) and the reaction system was made up of 4 μL of 5 × Reaction Buffer, 1 μL of primer, 1 μL of SweScript RT I Enzyme Mix, 0.1 ng-5 μg of total RNA, and nuclease-free water replenished to 20 μL. Afterthat, the qPCR was performed using the 2 × Universal Blue SYBR Green qPCR Master Mix following the manufacturer’s direction (Servicebio, Guangzhou, China) with the reaction system of 3 μL cDNA, 5 μL 2 × Universal Blue SYBR Green qPCR Master Mix and 1 μL each upstream and downstream primers. Finally, the reactions were performed on a CFX96 real-time quantitative fluorescence PCR instrument. The amplification reactions were programmed with pre-denaturation at 95 °C for 1 min, followed by 40 cycles, each cycle consisting of 95 °C for 20 s, 55 °C for 20 s, and 72 °C for 30 s. The relative expression of genes was calculated by the 2^−ΔΔCt method using GADPH as the internal reference gene. p-value was calculated by Graphpad Prism 5, and p less than 0.05 was considered to have statistical difference. The sequences of the primers used in this study are presented in Table 1, and the concentration of primer was 10 μM. The amplicon size was 80–350 bp.

Statistical analysis

The somatic mutation analysis and visualization were using the ‘maftools’ package in R. All bioinformatics analyses were undertaken in R language. The Wilcoxon test (comparison between two groups), Kruskal-Wallis test (comparison between three or more groups), and chi-square test (correlation analysis between cuproptosis-relevant subtypes and clinical factors) were employed to compare the data from different groups. Unless specified, a p-value < 0.05 was considered statistically significant.

Expression of ten cuproptosis-relevant genes in BRCA

First, we conducted the functional enrichment analysis of 10 cuproptosis genes (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A) using the DAVID web tool. As shown in Figs. 1A and 1B, 11 GO entries (5 BP, 4 CC, and 2 MF) and 8 KEGG pathways were enriched. These genes were involved in ‘acetyl-CoA biosynthetic process from pyruvate’, ‘mitochondrial acetyl-CoA biosynthetic process from pyruvate’, ‘tricarboxylic acid cycle’, ‘glucose metabolic process’, ‘protein lipoylation’, ‘pyruvate metabolism’, ‘glycolysis/gluconeogenesis’, ‘carbon metabolism’, ‘metabolic pathways’, ‘central carbon metabolism in cancer’, ‘lipoic acid metabolism’ and ‘biosynthesis of cofactors’. The locations of these 10 cuproptosis genes on the chromosomes are exhibited in the circle diagram (Fig. S1). Somatic mutation analysis of 986 samples with somatic mutation information from TCGA-BRCA cohort revealed a very low mutation rate of the 10 cuproptosis genes in BRCA from the index of tumbor mutational burden (TMB), thus there was a lower probability for natural immune response (NA) (Fig. 1C). Further analysis of transcriptomic data from TCGA-BRCA cohort revealed that the expressions of CDKN2A and PDHB were significantly elevated, and those of DLAT, DLD, FDX1, GLS, LIAS, LIPT1, MTF1, and PDHA1 was significantly lowered in BRCA samples compared to normal samples (p < 0.05) (Fig. 1D).

Identification of cuproptosis-associated subtypes of BRCA

Consensus clustering was implemented to identify BRCA-related subtypes, and 1,091 BRCA patients in TCGA-BRCA cohort were categorized into cuproptosis-associated subtypes based on the expression of 10 cuproptosis genes. Optimal clustering stability was confirmed at K = 3 (Figs. S2A–S2C). Clusters 1, 2, and 3 included 473, 296, and 322 samples, respectively. UMAP reduced dimensional analysis demonstrated that the samples of cluster 1 and cluster 2 were distributed separately and could be well distinguished (Fig. 2A). Survival analysis revealed significantly higher OS for patients in cluster 1 than in cluster 2 (p = 0.03688) (Fig. 2B). By chi-square tests, we examined the correlation between three cuproptosis-associated subtypes and clinical factors. The results showed that the three subtypes were significantly associated with age, race, and Node stage (N stage) but not with Tumor stage (T stage), Metastasis stage (M stage), and stage (Figs. 2C–2H). To further explore the discrepancies in the TMEs of the three cuproptosis-associated subtypes, we first calculated and compared the immune scores, stromal scores, estimate scores, and tumor purity across the three subtypes. As shown in Figs. 3A–3D, the stromal scores of all three subtypes were significantly different. The immune score, stromal score, and estimate score of cluster 1 and cluster 2 were significantly higher than those of cluster 3, and the tumor purity of cluster 1 and cluster 2 was significantly lower than that of cluster 3. Using the ssGSEA algorithm and Kruskal-Wallis test, we calculated and compared 28 immune gene sets across the three subtypes to further assess their differential immune activities. As shown in Figs. 3E and 3F, all 28 immune gene sets differed significantly among the three subtypes. To further compare the differences in fractions of immune cell infiltrates among the three subtypes, we employed the CIBERSORT algorithm and performed the Kruskal-Wallis test. The infiltration levels of a total of 14 immune cell types differed significantly among the three subtypes, including naïve B cells, memory B cells, activated dendritic cells, macrophages M1, macrophages M2, resting mast cells, neutrophils, activated NK cells, resting NK cells, resting CD4 memory T cells, activated CD4 memory T cells, CD8 T cells, follicular helper T cells, and regulatory T cells (Tregs) (Fig. 3G).

Cuproptosis-associated DEGs in BRCA

To authenticate the cuproptosis-associated DEGs in BRCA, we first identified the DEGs between BRCA and normal samples, between cluster 1 and cluster 2, between cluster 2 and cluster 3, and between cluster 1 and cluster 3. A total of 4,878 DEGs (2,441 up-regulated and 2,437 down-regulated genes) between BRCA and normal samples, 581 DEGs (270 up-regulated and 311 down-regulated genes) between cluster 1 and cluster 2, 705 DEGs (433 up-regulated and 272 down-regulated genes) between cluster 2 and cluster 3, and 1,221 DEGs (598 up-regulated and 623 down-regulated genes) between cluster 1 and cluster 3 were identified. By taking the intersection of the above four groups of DEGs using a Venn diagram, 38 common genes were defined as cuproptosis-associated DEGs in BRCA (Fig. 4A, Table S1). The results of Pearson correlation analysis between the 38 intersecting genes and 10 cuproptosis genes are shown in Fig. 4B. To further probe the function of these 38 genes in BRCA, a functional enrichment analysis was conducted. As shown in Table S2, 18 GO terms (13 BP, 2 CC, and 3 MF) and 1 KEGG pathway were significantly enriched. The top 10 items under each classification are shown in a bar diagram (Fig. 4C). The abovementioned genes were mainly linked to ‘response to iron ion’, ‘response to estradiol’, ‘response to estrogen’, ‘response to xenobiotic stimulus’, ‘neutrophil chemotaxis’, ‘response to drug’, ‘hormone metabolic process’, ‘mammary gland alveolus development’, and ‘estrogen signaling pathway’. We categorized the 1,091 BRCA patients in TCGA-BRCA cohort into five subtypes based on the expression of 38 cuproptosis-associated DEGs in BRCA and consensus clustering (Figs. S3A–S3C). Subtypes 1, 2, 3, 4, and 5 consisted of 225, 178, 210, 320, and 158 samples, respectively. The distribution of clinical features of the five subtypes and the expression of 38 cuproptosis-associated DEGs are shown as a heatmap (Fig. 4D). Survival analysis revealed no significant survival discrepancies among the five subtypes (Chi-square test = 6, df = 4, p = 0.2) (Fig. 4E).

Cuproptosis-relevant prognostic signature for BRCA

To establish a cuproptosis-relevant prognostic signature, we first identified six genes significantly associated with OS in BRCA patients among the 38 cuproptosis-associated DEGs in BRCA using univariate Cox analysis, namely SAA1, KRT17, VAV3, IGHG1, TFF1, and CLEC3A (Fig. 5A). Next, we performed a principal component analysis for the above six genes and calculated the PC1 and PC2 values for each BRCA sample (Table S3). The coefficients of PC1 and PC2 values were obtained by multivariate Cox regression analysis. According to the formula, Cusig score = 0.087984 * PC1 value + 0.226539 * PC2 value, we calculated the Cusig score for each BRCA sample and subsequently categorized the BRCA patients into the high-Cusig score group (552 cases) and low-CuSig score group (517 cases) according to the optimal cut-off value (0.9889137). The KM curve suggested that the survival of the patients in the high-Cusig score group was significantly worse than that of the low-Cusig score group (Fig. 5B). The Sankey diagram exhibited the correlation between the different subtypes and high- and low-Cusig scores (Fig. 5C). The Cusig scores were also significantly different between the five BRCA subtypes based on 38 cuproptosis-associated DEGs in BRCA (Fig. 5D). Cluster 2, which had a better prognosis, had a lower Cusig score, compared to the other clusters (Figs. 4E, 5D). We further validated the Cusig score model in external datasets, GSE42568 and GSE20711, and similarly, the high-Cusig score group showed a significantly worse prognosis than the low-Cusig score group, consistent with the results of TCGA-BRCA cohort (Figs. 5E, 5F).

Uncovering the molecular mechanisms of cuproptosis-relevant prognostic signature underlying BRCA

In order to uncover potential mechanisms underlying the differential prognoses between the two Cusig score subgroups, GSVA was conducted to analyze the enrichment differences in the terms of KEGG and hallmark pathways between different Cusig score groups. As shown in Table S4 and Figs. 6A, 6B, 22 hallmark pathways and 73 KEGG pathways were enriched. ‘E2F targets’, ‘G2M checkpoint’, ‘protein secretion’, ‘unfolded protein response’, ‘MTORC1 signaling’, ‘MYC targets V1’, and ‘DNA repair’ were enriched in the high-Cusig score group (Fig. 6A). ‘Citrate cycle’, ‘TCA cycle’, ‘RNA degradation’, ‘cell cycle’ and ‘DNA replication’ were also enriched in the high-Cusig score group (Fig. 6B). Immune-related hallmark pathways, like ‘IL2-STAT5 signaling’, ‘inflammatory response’, ‘IL6-JAK-STAT3 signaling’, ‘TNFA signaling via NFKB’, and KEGG pathways, including ‘chemokine signaling pathway’, ‘NOD-like receptor signaling pathway’, ‘TGF-β signaling pathway’, ‘complement and coagulation cascades’, ‘natural killer cell-mediated cytotoxicity’, ‘cytokine-cytokine receptor interaction’, and ‘B cell receptor signaling pathway’ were enriched in low-Cusig score group (Figs. 6A, 6B). Next, we screened DEGs of the Cusig score subgroups and conducted GO functional enrichment analysis. As shown in Table S5, 1,289 GO entries (1,152 BP, 67 CC, and 70 MF) were derived based on 590 DEGs between the high- and low-Cusig score groups. The top 10 GO-BP terms are shown in Fig. 6C. Immune-related and cell adhesion-related BPs were associated with these DEGs. The results of GSEA indicated that 484 and 1,686 GO entries (NES > 0 and NES < 0, respectively) were enriched in the high-Cusig and low-Cusig score groups, respectively (Table S6). The top five enriched entries of the two Cusig score subgroups are shown in Fig. 6D. Notably, mitosis-related and DNA replication-related BPs were closely associated with the high-Cusig score group, and multiple immune-related BPs were closely associated with the low-Cusig score group.

Association of cuproptosis-relevant prognostic signature with TME

Since immune-related pathways and BPs were found to be associated with the low-Cusig score group, we next analyzed the relationship between the Cusig score and TME and immune cell infiltration. The ESTIMATE algorithm demonstrated that patients with low Cusig scores had higher immune scores, stromal scores, estimate scores, and lower tumor purity (Figs. 7A–7D). The results of the CIBERSORT algorithm showed that the fractions of naïve B cells, plasma cells, CD8 T cells, monocytes, resting dendritic cells, eosinophils, and neutrophils were elevated in the patients with lower Cusig scores, while the fractions of macrophages M0 and macrophages M2 were high in the patients with high Cusig scores (Fig. 7E). Moreover, PD-L1 expression and immunophenoscore (IPS) of the low-Cusig score group were higher than those of the high-Cusig score group (Figs. 7F, 7G). However, the results of the TIDE analysis showed that the high-Cusig score group had a lower TIDE score than the low-Cusig score group and these patients may be more sensitive to ICB therapy (Fig. 7H).

The mutational landscape of high- and low-Cusig score samples

To further assess the mutation frequency of genes in high- and low-Cusig score samples, we performed somatic mutation analysis. The most frequently mutated genes in the high-Cusig score and low-Cusig score samples were TP53 and PIK3CA, respectively (Figs. 8A, 8B). The top 20 mutated genes were not identical between the two Cusig score subgroups, indicating different somatic mutation patterns (Figs. 8A, 8B). We then analyzed the copy number variation (CNV) in genes of the two Cusig score subgroups and found that the CNV of genes in the low-Cusig score samples were significantly higher than that in the high-Cusig score samples (Fig. 8C).

Expression of cuproptosis-relevant prognostic genes in BRCA

As shown in Fig. S4, CLEC3A, IGHG1, TFF1, and VAV3 were up-regulated, while KRT17 and SAA1 were down-regulated in BRCA tissues compared to normal tissues in TCGA-BRCA cohort. We then checked the expression of prognostic genes at the mRNA level in the human epithelial cell line from the mammary gland, MCF-10A, and three breast cancer cell lines HCC1937, MCF7, and MDA-MB-231. Consistent with the trend of results from public databases, levels of CLEC3A and IGHG1 were up-regulated, while those of KRT17 and SAA1 were down-regulated in breast cancer cell lines (Fig. 9). The amplification and dissolution curves of prognostic genes were displayed in Table S7. However, inconsistent with the results from the tissue samples, TFF1 and VAV3 were down-regulated in breast cancer cell lines (Fig. 9, Table S7). The amplification and dissolution curves of prognostic genes were displayed in Table S8. However, inconsistent with the results from the tissue samples, TFF1 and VAV3 were down-regulated in breast cancer cell lines (Fig. 9), probably due to the complexity of the tumor tissue.