The expression of ten cuproptosis genes in BRCA
First, through the online site DAVID, we conducted the functional enrichment analysis of 10 cuproptosis genes (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS and CDKN2A). As shown in Fig. 1A-B, a total of 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 positions of these 10 cuproptosis genes on the chromosomes were exhibited in the circle diagram (Fig. 1C). Somatic mutation analysis of TCGA-BRCA cohort revealed the very low mutation rate of 10 cuproptosis genes in BRCA samples (Fig. 1D). Further analysis of transcriptomic data from TCGA-BRCA cohort revealed that the expression of CDKN2A and PDHB was significantly elevated in BRCA samples compared to normal samples, and the expression of DLAT, DLD, FDX1, GLS, LIAS, LIPT1, MTF1 and PDHA1 was significantly decreased in BRCA samples compared to normal samples (Fig. 1E).
Recognition Of Cuproptosis-associated Subtypes Of Brca
In the following, we categorized the 1091 BRCA patients in the TCGA-BRCA cohort into three cuproptosis-associated subtypes based on the expression of 10 cuproptosis genes and consensus clustering (Fig. 2A-C). The cluster 1 included 473 samples, cluster 2 included 296 samples, and cluster 3 included 322 samples. The results of UMAP reduced dimensional analysis demonstrated that the samples of three subtypes were distributed in different areas and could be clearly distinguished, which further illustrated the accuracy of the subtyping (Fig. 2D). Survival analysis revealed significant survival discrepancies between the three cuproptosis-associated subtypes, with more remarkable survival superiority for Cluster 1 (Fig. 2E). By means of chi-square tests, we explored the correlation between three cuproptosis-associated subtypes and clinical factors. The results showed that the three subtypes were significantly associated with age, race and N stage, but not with T stage, M stage and stage (Fig. 2F-K). To further explore the discrepancies in the TME of the three cuproptosis-associated subtypes, we first calculated and compared the immune scores, stromal scores, estimate scores and tumor purity of the three subtypes. As displayed in Fig. 3A-D, the stromal scores of all three subtypes were significantly different from each other. 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 1and cluster 2 was significantly lower than that of cluster 3. Through ssGSEA algorithm and Kruskal-Wallis test, we calculated and compared the fractions of 28 immune gene sets of the three subtypes to further speculate on their differential immune activity. As exhibited in Fig. 3E-F, the fraction of all 28 immune gene sets differed significantly among three subtypes. To further compare the differences in immune cell infiltration among the three subtypes, we adopted the CIBERSORT algorithm and the Kruskal-Wallis test. The results disclosed that 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).
The Cuproptosis-associated Degs In Brca
To authenticate the cuproptosis-associated DEGs in BRCA, we firstly 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 respectively. A total of 4878 DEGs (2441 up-regulated genes and 2437 down-regulated genes) between BRCA and normal samples, 581 DEGs (270 up-regulated genes and 311 down-regulated genes) between cluster 1 and cluster 2, 705 DEGs (433 up-regulated genes and 272 down-regulated genes) between cluster 2 and cluster 3, and 1221 DEGs (598 up-regulated genes and 623 down-regulated genes) between cluster 1 and cluster 3 were determined. By taking the intersection of the above 4 groups of DEGs, we acquired 38 intersecting genes, defined as cuproptosis-associated DEGs in BRCA (Fig. 4A, Supplementary Table 1). To further probe the function of these 38 genes in BRCA, functional enrichment analysis was executed. As displayed in Supplementary Table 2, 18 GO items (13 BP items, 2 CC items, and 3 MF items) and 1 KEGG pathway were derived. The top10 items under each classification were showcased in bar diagram (Fig. 4B). We observed that aforesaid 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’. Then, we categorized the 1091 BRCA patients in the TCGA-BRCA cohort into five subtypes based on the expression of 38 cuproptosis-associated DEGs in BRCA and consensus clustering (Fig. 4C-E). The cluster 1 included 225 samples, cluster 2 included 178 samples, cluster 3 included 210 samples, cluster 4 included 320 samples, and cluster 5 included 158 samples. The distribution of clinical features of the five subtypes and the expression of 38 cuproptosis-associated DEGs were illustrated in the heatmap (Fig. 4F). Survival analysis revealed significant survival discrepancies between the five subtypes (Fig. 4G).
The Cuproptosis-relevant Prognostic Signature In Brca
To establish a cuproptosis-relevant prognostic signature, we first detected 6 genes significantly associated with OS in BRCA patients based on 38 cuproptosis-associated DEGs in BRCA using univariate Cox analysis, namely SAA1, KRT17, VAV3, IGHG1, TFF1 and CLEC3A (Fig. 5A). Next, we performed principal component analysis based on the above six genes and calculated PC1 and PC2 values for each BRCA sample (Supplementary Table 3). The coefficients of PC1 values and PC2 values were obtained by multivariate Cox regression analysis. According to the formula of CuSig score, Cusig score = 0.087984*PC1 value + 0.226539*PC2 value, we calculated the Cusig score for each BRCA sample and divided the BRCA patients into high-Cusig score group (552 cases) and low-CuSig score group (517 cases) according to the optimal cut-off value (0.9889137). The K-M curve suggested that the survival of 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 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 (Fig. 4G, 5D). We further validated the Cusig score model in the external dataset GSE42568, and similarly, the high-Cusig score group had a significantly worse prognosis than the low-Cusig score group, consistent with the results of TCGA-cohort (Fig. 5E).
To further understand the mutation frequency of genes in high- and low-Cusig score samples, we performed somatic mutation analysis. The most frequently mutated gene in the high-Cusig score sample was TP53, and in the low-Cusig score sample, the most frequently mutated gene was PIK3CA (Fig. 5F-G). In addition, the top20 mutated genes were not identical in two Cusig score subgroups, indicating different somatic mutation patterns in the two Cusig score subgroups (Fig. 5F-G). We then analyzed the copy number variation (CNV) of genes in the two Cusig score subgroups and demonstrated that the CNV of genes in the low-Cusig score samples was significantly higher than that in the high-Cusig score samples (Fig. 5H).
Uncovering The Molecular Mechanisms Involved In The Cuproptosis-relevant Prognostic Signature In Brca
In order to uncover potential mechanisms for the disparity in prognosis between the two Cusig score subgroups, GSVA was conducted to analyze the enrichment differences in the terms of KEGG and hallmark between different Cusig score groups. As listed in Supplementary Table 4 and Fig. 6A-B, a total of 22 hallmark pathways and 73 KEGG pathways were enriched. The ‘E2F targets’, ‘G2M checkpoint’, ‘protein secretion’, ‘unfolden protein response’, ‘MTORC1 signaling’, ‘MYC targets V1’ and ‘DNA repair’ were enriched in high-Cusig score group (Fig. 6A). The ‘Citrate cycle’, ‘TCA cycle’, ‘RNA degradation’, ‘cell cycle’ and ‘DNA replication’ were also enriched in high-Cusig score group (Fig. 6B). Immune-related hallmark pathways, such as ‘IL2-STAT5 signaling’, ‘inflammatory response’, ‘IL6-JAK-STAT3 signaling’, ‘TNFA signaling via NFKB’, and KEGG pathways, such as ‘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 (Fig. 6A-B). Next, we screened DEGs from two Cusig score subgroups and carried out GO functional enrichment analysis. As shown in Supplementary Table 5, a grand total of 1289 GO entries (1152 BP entries, 67 CC entries, and 70 MF entries) were derived based on 590 DEGs between high- and low-Cusig score group. The top10 GO-BP entries were exhibited in Fig. 6C. Immune-related and cell adhesion-related biological processes were associated with these DEGs. The GSEA result indicated that 484 GO entries (NES > 0) were enriched in the high-Cusig score group and 1686 GO entries (NES < 0) were enriched in the low-Cusig score group (Supplementary Table 6). The top5 enriched entries for the two Cusig score subgroups were drawn into Fig. 6D. We noted that mitosis-related and DNA replication-related biological processes were closely associated with high-Cusig score group, and multiple immune-related biological processes were closely associated with the low-Cusig score group.
Association Of Cuproptosis-relevant Prognostic Signature With Tme
Since immune-related pathways and biological processes were associated with the low-Cusig score group, we next analyzed the relationship between Cusig score and TME and immune cell infiltration. The ESTIMATE algorithm demonstrated that low-Cusig score patients had higher immune scores, stromal scores, estimate scores and lower tumor purity (Fig. 7A-D). The CIBERSORT algorithm uncovered that the fraction of naïve B cells, Plasma cells, CD8 T cells, Monocytes, resting Dendritic cells, Eosinophils, and Neutrophils was elevated in the patients with lower Cusig scores, while the fraction of Macrophages M0 and Macrophages M2 were superior in the patients with higher Cusig scores (Fig. 7E). In addition, the PD-L1 expression and immunophenoscore (IPS) of low-Cusig score group were higher than high-Cusig score group (Fig. 7F-G). However, the results of TIDE analysis showed that the high-Cusig score group had a lower TIDE score than the low-Cusig score group and may be more sensitive to immune checkpoint blockade therapy (Fig. 7H).
The Expression Of Cuproptosis-relevant Prognostic Genes In Brca
As exhibited in Supplementary Fig. 1, CLEC3A, IGHG1, TFF1 and VAV3 were up-regulated in BRCA tissues, and KRT17 and SAA1 were down-regulated in BRCA tissues compared to normal tissues in the TCGA-BRCA cohort. We then check the expression of prognostic genes at the mRNA level in 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, CLEC3A and IGHG1 was up-regulated in breast cancer cell lines and KRT17 and SAA1 were down-regulated in BRCA cell lines (Fig. 8). However, inconsistent with the results of the tissue samples, TFF1 and VAV3 were down-regulated in breast cancer cell lines (Fig. 8), probably due to the complexity of the tumor tissue.