Shaoran Song1,2†
Yaochun Wang- 1Center for Translational Medicine, The First Affiliated Hospital of Xi’an Jiaotong University, Xi’an, China
- 2The Key Laboratory for Tumor Precision Medicine of Shaanxi Province, The First Affiliated Hospital, Xi’an Jiaotong University, Xi’an, China
- 3Department of Breast Surgery, The First Affiliated Hospital of Xi’an Jiaotong University, Xi’an, China
- 4Department of Oncology, The First Affiliated Hospital of Xi’an Jiaotong University, Xi’an, China
Background: Cuproptosis is a newly discovered programmed cell death dependent on overload copper-induced mitochondrial respiration dysregulation. The positive response to immunotherapy, one of the most important treatments for invasive breast cancer, depends on the dynamic balance between tumor cells and infiltrating lymphocytes in the tumor microenvironment (TME). However, cuproptosis-related genes (CRGs) in clinical prognosis, immune cell infiltration, and immunotherapy response remain unclear in breast cancer progression.
Methods: The expression and mutation patterns of 12 cuproptosis-related genes were systematically evaluated in the BRCA training group. Through unsupervised clustering analysis and developing a cuproptosis-related scoring system, we further explored the relationship between cuproptosis and breast cancer progression, prognosis, immune cell infiltration, and immunotherapy.
Results: We identified two distinct CuproptosisClusters, which were correlated with the different patterns between clinicopathological features, prognosis, and immune cell infiltration. Moreover, the differences of the three cuproptosis-related gene subtypes were evaluated based on the CuproptosisCluster-related DEGs. Then, a cuproptosis-related gene signature (PGK1, SLC52A2, SEC14L2, RAD23B, SLC16A6, CCL5, and MAL2) and the scoring system were constructed to quantify the cuproptosis pattern of BRCA patients in the training cohort, and the testing cohorts validated them. Specifically, patients from the low-CRG_score group were characterized by higher immune cell infiltration, immune checkpoint expression, immune checkpoint inhibitor (ICI) scores, and greater sensitivity to immunotherapy. Finally, we screened out RAD23B as a favorable target and indicated its expression was associated with breast cancer progression, drug resistance, and poor prognosis in BRCA patients by performing real-time RT-PCR, cell viability, and IC50 assay.
Conclusions: Our results confirmed the essential function of cuproptosis in regulating the progression, prognosis, immune cell infiltration, and response to breast cancer immunotherapy. Quantifying cuproptosis patterns and constructing a CRG_score could help explore the potential molecular mechanisms of cuproptosis regulating BRCA advancement and provide more effective immunotherapy and chemotherapy targets.
Introduction
Breast cancer (BRCA) is one of the most common malignant cancers among women worldwide, with a high incidence and recurrence rate (1, 2). According to the latest statistics from 2022, breast cancer alone accounts for nearly one-third of all new cancer diagnoses (287,850 new cases) in women in the United States (3). Despite diagnostic and therapeutic strategies that have taken into consideration the heterogeneity of breast cancer (4), there are presently insufficient strategies to improve the prognosis of recurrence-free survival (RFS) and overall survival (OS) in breast cancer. At the same time, the resistance of breast cancer patients to chemotherapy, radiotherapy, or endocrine therapy makes maintaining their long-term survival an urgent challenge. In recent years, attention has been focused on the role of the tumor microenvironment (TME) in regulating breast cancer progression and prognosis and the effect of immunotherapy in breast cancer treatment (5–7). However, there is still a lack of sensitive immune-related diagnostic and therapeutic targets for breast cancer.
Cuproptosis, unlike apoptosis, ferroptosis, pyroptosis, and necroptosis, is a kind of non-apoptotic programmed cell death induced by the accumulation of intracellular copper (8, 9). Direct binding to lipid-acylated mitochondrial proteins of the tricarboxylic acid (TCA) cycle to aggregate them, followed by proteotoxic stress, is the crucial mechanism for the initiation of cuproptosis. Previous studies have illustrated the relationship between copper homeostasis and human diseases, including Wilson disease and other neurological copper disorders (10, 11), cancers (12–15), abnormal fetal development (16), and so on. However, there are no studies on the association between the newly defined cuproptosis and breast cancer oncogenesis, immune microenvironment, or immunotherapy. Therefore, exploring the physiological and pathological activities associated with cuproptosis, elucidating its underlying mechanisms affecting breast cancer progression, and identifying sensitive and effective targets for breast cancer diagnosis and treatment are crucial for early detection, diagnosis, and treatment of breast cancer.
In this study, a comprehensive assessment of the expression profile of 12 cuproptosis-related genes in breast cancer was performed to comprehensively analyze the role of CRGs on TME and immunotherapy. First, the BRCA patients in the training cohort were stratified into two cuproptosis-related clusters based on CRGs expression levels. The differentially expressed genes (DEGs) in these two clusters were then used to classify the patients into three cuproptosis-related gene subtypes. Further, differentially expressed genes with prognostic significance were used to construct a cuproptosis-related gene signature (PGK1, SLC52A2, SEC14L2, RAD23B, SLC16A6, CCL5, and MAL2) and scoring system. Three independent external testing cohorts also confirmed the stability and reliability of the scoring system. We used the cuproptosis-related gene score (CRG_score) to classify patients into high and low-CRG_score subgroups to predict overall survival (OS) and the immune landscape in BRCA, thus accurately predicting the patient long-term prognosis and response to immunotherapy. Finally, RAD23B was selected as a valuable target for in vitro experimental validation.
Materials and methods
Data acquisition and preprocessing
The gene expression profile cohort and its corresponding clinical data of BRCA patients were obtained from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) and the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) databases. Specifically, the training cohort consisted of the BRCA dataset of TCGA (113 standard samples and 1091 BRCA samples) and GSE20685 (327 BRCA samples) (Table S1). The testing cohorts consisted of the GSE7390 (198 BRCA samples), the GSE58812 (107 BRCA samples), and the GSE42568 (104 BRCA samples). All the gene expression data were fragments per kilobase million (FPKM) and transformed into transcripts per kilobase million (TPM) values for further analysis using the R language (version 4.1.2), the “edge” R package (Storey et al., 2021, R package version 2.26.0). Reducing the batch impact induced by non-biotechnological variations, was achieved by using the “ComBat” method in the “SVA” R package (Leek et al., 2021, R package version 3.42.0.).
Construction of CuproptosisClusters and PCA analysis
The 12 cuproptosis-related genes, including FDX1, LIPT1, LIAS, DLD, DBT, GCSH, DLAT, PDHA1, PDHB, SLC31A1, ATP7A, ATP7B, were retrieved from previously published literature (8). To identify different cuproptosis patterns in BRCA, we performed consensus classification using the “ConsensusClusterPlus” R package (17). The tendency and smoothness of the cumulative distribution function (CDF) curve were used to figure out the clustering number (17). Principal component analysis (PCA) was conducted with the help of the function “prcomp” in the R package “stats” (R Core Team, 2021).
Clinical characteristics and prognosis analysis in different CuproptosisClusters
Utilizing the “survival” (Therneau et al., 2021, R package version 3.2-13) and “survminer” R packages (Kassambara et al., 2021, R package version 0.4.9.), we conducted the Kaplan–Meier plot to estimate the prognostic values of BRCA patients in different molecular subtypes (18, 19). Clinical features (age, stage T, and stage N) were also compared among molecular subtypes.
Gene set variation analysis in different CuproptosisClusters
Using the “GSVA” R package (20), the gene set variation analysis (GSVA) was conducted to estimate the differences in biological processes responsible for the characteristic patterns of cuproptosis (20, 21).
Estimation of TME in different CuproptosisClusters
By performing the “estimate” R package (Yoshihara et al., 2016, R package version 1.0.13/r21.), the immune score of every BRCA sample was evaluated in the ESTIMATE algorithm (22). Furthermore, the fractions of 23 human immune cell subsets in the TME of different BRCA molecular subtypes were also estimated using a single-sample gene set enrichment analysis (ssGSEA) algorithm (23–26). The expression of 32 critical immune checkpoints retrieved from previous research was compared in different BRCA molecular subtypes (27).
Identification of DEGs in different CuproptosisClusters, functional analysis, and construction of gene subtypes
We identified 968 DEGs among the different CuproptosisClusters using the “limma” package (28) in R with a fold-change of 6 and an adjusted p < 0.001. The functional analysis (GO and KEGG) was conducted on the DEGs using the “clusterprofiler” R package (17, 29, 30). The gene set file (c2.cp.kegg.v7.2.symbols.gmt) was obtained from the MSigDB database (https://www.gsea-msigdb.org) (30, 31). To investigate the molecular function of these cuproptosis-related DEGs mentioned above, we performed survival analysis and picked out 25 DEGs with significant prognostic values (p < 0.001) for further study (Table S8) (28, 32). Then, consensus classification was performed using the ConsensusClusterPlus” R package based on the 25 prognostic genes to divide patients into three gene subtypes (Figure S2, k = 3; gene subtypes A, B, and C).
Development and validation of the cuproptosis-related gene model and CRG_score
We performed the LASSO Cox regression model using the “glmnet” R package (33, 34) to filter down the candidate cuproptosis-related genes. Finally, the 5 genes and their coefficients were kept. We obtained the penalty parameter (λ) according to the minimum criteria. After standardizing the data from the training cohort by the “scale” R package (Hadley Wickham and Dana Seidel, 2022, R package version 1.2.0.), the CRG_score was calculated as follows:
CRG_score = (0.00523335734652904 * PGK1) + (0.0185186220360186 * SLC52A2) − (0.0261649532232335 * SEC14L2) + (0.0129284730002406 * RAD23B) − (0.0297061289632435 * SLC16A6) − (0.0142610013224485 * CCL5) + (0.0014557255170858 * MAL2).
We also calculated the CRG_score of the testing cohorts using the same formula. We divided the patients from the training and testing cohorts based on the median CRG_score, and the high- and low-CRG_score groups. We performed the Kaplan–Meier analysis of overall survival using “survival” and “survminer” R packages, and ROC curve analysis using the “timeROC” R package (35). Then the calibration plots of the nomogram were executed to predict the prognosis value between the predicted 3-, 5-, and 8- or 10-year survival events and the virtually observed outcomes. Lastly, a stratified analysis was done to see if the CRG_score could still predict in different subgroups of age (65 or 65), stage T (T1-2 or T3-4), and stage N (N0 or N1-3).
Finally, we reduced the batch impact of testing cohorts (GSE7390, GSE58812, and GSE42568; n = 409) using “ComBat” method. By performing the PAM50 algorithm (36) with the “genefu” R package (37), we then assessed the molecular subtypes for each patient from testing cohorts. Samples were classified into the normal-like (n = 13), basal-like (n = 126), HER2+ (n = 57), luminal A (n = 111), and luminal B (n = 102) subtypes (Table S2). Each subtype was then used as an independent external validation cohort and CRG_score was calculated respectively. Survival and ROC curve analysis were then performed, as mentioned above.
Analysis of chemotherapeutic drugs effects in high-and low-CRG_score groups
The semi-inhibitory concentration (IC50) values of chemotherapeutic drugs in high- and low-CRG_score groups were calculated using the “pRRophetic” R package (Paul Geeleher, 2014, R package version 0.5.).
Analysis of protein expression in clinical specimen
The Human Protein Atlas (HPA) database contains sections from 46 normal human tissues and over 20 human cancers labeled with antibodies targeting more than 11000 human proteins (38). Based on the laser power and detector gain parameters used for image acquisition, combined with the image’s visual appearance, the staining intensity is rated as negative, weak, moderate and strong (39). The scoring method of protein expression is the same as described previously (40).
Cell culture and transfection
Our human breast cancer cell lines were obtained from the Shanghai Institute of Biochemistry and Cell Biology, including MCF10A, SUM-159, MDA-MB-231, BT549, and MCF7 cells. MCF10A and SUM159 cells were cultured in DMEM/F12 (1:1) medium, MDA-MB-231 and MCF7 cells were cultured in a DMEM medium, and BT549 cells were cultured in an RPMI-1640 medium, with all recommended supplements, respectively. All cells were cultured at 37°C in a humidified incubator in an atmosphere of 5% CO2.
SUM-159 and MCF7 cells were transfected with corresponding siRNAs using Lipofectamine 8000 (#C0533, Beyotime, Nanjing, China) following the manufacturer’s protocol. The RAD23B siRNA constructs and a negative control siRNA were as follows: RAD23B-1, 5’- CUCCAGCAUCAGCGACAGCAUTT −3’ and 5’- AUGCUGUCGCUGAUGCUGGAGTT −3’; RAD23B-2, 5’- AGAAGCUGGAAGUGGUCAUAUTT −3’ and 5’- AUAUGACCACUUCCAGCUUCUTT −3’; NC siRNA, 5’- UUCUCCGAACGUGUCACGUdTdT −3’ and 5’- ACGUGACACGUUCGGAGAAdTdT −3’.
Cuproptosis cell model construction and mRNA expression analyses
To promote the occurrence of Cuproptosis (8), SUM-159 and MCF7 cells were treated with 100 nM elesclomol (+1 µM CuCl2 in medium) for a 2-hour-pulse. After 24 h, cells were harvested and lysed. The real-time RT-PCR was performed as previously described (40). All primers were synthesized (Sangon Biotech, Shanghai, China) and listed in Table S3.
Cell viability and IC50 assay
After transfection for 48h, SUM-159 and MCF7 cells (5 × 103 cells/well) were loaded on a 96-well plate and cultivated for 0 h, 24 h, 48 h, and 72 h for cell viability assays or 48 h for the IC50 assay of Paclitaxel (0, 0.125, 2.5, 5, 10, 20 µM/L). After incubation with 20 μL of 3-(4.5-Dimethylghiazol-2-yl)-2,5-diphenyltetrazolium Bromide (MTT; 5 mg/mL; Absin Bioscience, Shanghai, China; Catalog no. abs50010) for 4 h at 37°C, the culture medium was removed, and 150 μL dimethyl sulfoxide (DMSO; Sigma, St. Louis, MO, USA) was added. Afterward, the cells were shaken for 15 min in the dark, and the optical density (OD) at 490 nm was measured using a Benchmark microplate reader (Bio-Rad, Hercules, CA, USA).
Statistical analysis
All the data analysis was conducted by R (version 4.1.2) and in vitro experimental data were analyzed with the GraphPad Prism 5 (GraphPad Software, Inc., La Jolla, CA, USA). The Log-rank test, Spearman test, Wilcoxon test, Student’s t-test, and Two-way ANOVA tests were applied in this study. p < 0.05 was considered as significant. P-values were adjusted to control for the false discovery rate (FDR) using the Benjamini-Hochberg method (41). Each experiment was done in triplicate and repeated at least three times.
Results
The landscape of cuproptosis-related genes in BRCA
Figure 1 shows the overall design and flow chart of this study. The expression of twelve cuproptosis-related genes (CRGs) was obtained from previous studies (8). By performing the “limma” package, we analyzed the mRNA expression of CRGs based on the data of 113 normal and 1091 BRCA tissues from TCGA. As shown in Figure 2A, the transcriptional levels of DLAT, PDHB, SLC31A1, and ATP7B were significantly higher in the BRCA tissues than in the normal tissues. At the same time, FDX1, LIPT1, LIAS, DLD, DBT, GCSH, PDHA1, and ATP7A were significantly lower. Furthermore, we developed the correlation network containing twelve CRGs in Figure 2B (red: positive correlations; blue: negative correlations). Then, univariable Cox regression analysis also showed significant differences in overall survival between patients with high or low expression of the CRGs (Table S4). Specifically, LIPT1 (HR = 0.79, 95% CI: 0.65–0.95, p < 0.01), PDHB (HR = 0.97, 95% CI: 0.94–1.00, p < 0.01), and ATP7B (HR = 0.96, 95% CI: 0.93–1.00, p < 0.01) were “protective” factors for BRCA patients with HR < 1, while the DLAT (HR = 1.06, 95% CI: 1.02–1.09, p < 0.01), SLC31A1 (HR = 1.03, 95% CI: 1.01–1.05, p < 0.01), DBT (HR = 1.05, 95% CI: 0.94–1.18, p < 0.01), PDHA1 (HR = 1.01, 95% CI: 0.99–1.03, p < 0.05), ATP7A (HR = 1.00, 95% CI: 0.93–1.09, p < 0.05), and DLD (HR = 1.02, 95% CI: 0.99–1.04, p < 0.05) were “risk” factors with HR > 1.
Figure 1 The Flowchart of the Study Design.
Figure 2 The Landscape of Cuproptosis-Related Genes in BRCA. (A) The gene expression levels of CRGs in BRCA compared to normal tissue (Wilcoxon test; blue: normal; red: BRCA). (B) The correlation network of the 12 CRGs (red: positive correlation; blue: negative correlation). (C) The frequency of CNV variation in CRGs (green: CNV deletion; red: CNV amplification). (D) The location of the CNV alteration of the CRGs changes on 23 chromosomes. (E) The genetic alteration on a query of CRGs. CRGs, cuproptosis-related genes; BRCA, breast cancer; CNV, copy number variant. *p < 0.05, **p < 0.01, ***p < 0.001.
To explore their mutation landscape, the single nucleotide variation (SNV) and copy number variation (CNV) data were downloaded from the TCGA database. Moreover, CNVs were prevalent and mostly involved deletion, though DLD had a high frequency of amplification (Figure 2C). Figure 2D shows the locations of the CNV alterations in the CRGs on their respective chromosomes. Next, we analyzed the incidence of somatic mutations in these 12 CRGs, which showed that 29 (2.94%) of the 986 BRCA samples had mutations in the CRGs. Specifically, ATP7A had the highest mutation frequency (1%), followed by ATP7B, while others did not have any significant mutations (Figure 2E). Therefore, our findings on the landscape of CRGs in BRCA showed that they might play an essential role in the development and progression of BRCA.
Identification of CuproptosisClusters in BRCA
To further explore the expression pattern of CRGs implicated in BRCA tumorigenesis, we integrated patients from TCGA-BRCA (n = 1091) and GSE20685 (n = 327) as training cohort. First, the comprehensive landscape of CRG interactions and their prognostic value in the BRCA training cohort was demonstrated in a cuproptosis network (Figure 3A). We did the unsupervised clustering analysis using the “ConsensusClusterPlus” R package and picked k= 2 based on the empirical cumulative distribution function (CDF) plots, which suggested the highest intragroup correlations and the lowest intergroup correlations compared with others (Figures 3B, C, and Figure S1). Thus, two CRG-expression patterns were observed: CuproptosisCluster A and CuproptosisCluster B. Furthermore, the BRCA patients in the training cohort could be completely distinguished (Figure 3D). We also performed the Kaplan–Meier (K-M) survival analysis of the two clusters, suggesting a poor overall survival of patients in Cluster B (Figure 3E). Finally, we examined the clinical and pathological characteristics of the two BRCA clusters and the expression of the CRGs (Figure 3F).
Figure 3 Identification of CuproptosisClusters in BRCA. (A) The interaction among CRGs in BRCA (green: favorable factors for overall survival; purple: risk factors for overall survival). (B) The relative change in area under consensus CDF curve for k=2 to 9. (C) The consensus clustering of BRCA patients for k = 2. (D) The PCA analysis of the two CuproptosisClusters. (E) The OS of the two CuproptosisClusters (Log-rank test). (F) The heatmap for the connections between clinicopathologic features and the two CuproptosisClusters (blue: low expression; red: high expression). CDF, cumulative distribution function; PCA, Principal component analysis; OS, overall survival. *p < 0.05, **p < 0.01, ***p < 0.001.
Characteristics of the TME in CuproptosisClusters of BRCA
To thoroughly analyze the role of cuproptosis-related genes in the TME of BRCA, we conducted the GSVA enrichment analysis. As shown in Figure 4A, cluster A was significantly enriched in progesterone-mediated oocyte maturation, cell cycle, oocyte meiosis, cysteine and methionine metabolism, basal transcription factors, homologous recombination, DNA replication, mismatch repair, glycosphingolipid biosynthesis lacto and neolacto series, pathogenic escherichia coli infection, proteasome, pyruvate metabolism, pentose phosphate pathway, glycolysis gluconeogenesis, pyruvate metabolism, citrate cycle TCA cycle, terpenoid backbone biosynthesis, amino sugar and nucleotide sugar metabolism. Moreover, the results of the ssGSEA algorithm indicated that CuproptosisClusters A and B were rich in different innate immune cell infiltrations with significance (Figure 4B). We also examined the expression of 32 immune checkpoints in two clusters, which showed a higher expression of most immune checkpoints (BTLA, CDC20R1, CD244, CD27, CD274, CD28, CD40, CD40LG, CD48, CD80, CTLA4, HHLA2, ICOS, IDO1, IDO2, KIR3DL1, LAG3, LGALS9, PDCD1, TIGIT, TMIGD2, TNFRSF9, TNFSF14) in cluster B (Figure 4C). In addition, NRP1 and TNFRSF14 were highly expressed in cluster A. Using the “estimate” package, we evaluated the TME score in Figure 4D. Cluster A had a higher stromal score, while Cluster B had a higher immune score. According to the results above, we identified two clusters with distinct immunological and metabolic characteristics, suggesting that cuproptosis may affect the immune microenvironment and metabolic processes that lead to breast cancer progression.
Figure 4 Characteristics of the TME in CuproptosisClusters of BRCA. (A) The GSVA of biological pathways between the two CuproptosisClusters (Spearman test, blue: inhibited pathways; red: activated pathways). (B) The abundance of 23 infiltrating immune cell types in the two CuproptosisClusters (Spearman test). (C) The expression levels of 32 immune checkpoints in the two CuproptosisClusters (Wilcoxon test). (D) The TME score of the two CuproptosisClusters (Spearman test). GSVA, gene set variation analysis; TME, tumor microenvironment. **p < 0.01, ***p < 0.001.
Identification of gene subtypes based on CuproptosisClusters of BRCA
Next, we identified 968 CuproptosisCluster-related DEGs by performing the “limma” package further to explore the different biological behaviors of each cluster (Table S5). Firstly, functional enrichment and GO (Gene Ontology) analysis and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis were conducted among the CuproptosisCluster-related genes (Figure 5A, B, Table S6 and S7). Then a univariable Cox regression analysis was performed to screen out 25 genes with significant prognostic values (p < 0.001) for the subsequent investigation (Table S8). Moreover, we applied the consensus to divide patients into three gene subtypes based on the 25 prognostic genes (Figure S2). As shown in Figure 5C, patients of gene subtype A showed the best OS, while patients of gene cluster C showed the worst OS (p < 0.001). The comparison of the clinicopathological characteristics and the expression of DEGs between the three gene subtypes was shown in the heatmap (Figure 5D). Finally, we observed the different expressions of the cuproptosis-related genes in the three gene subtypes (Figure 5E).
Figure 5 Identification of Gene Subtypes based on CuproptosisClusters of BRCA. (A, B) The GO and KEGG enrichment analyses of DEGs among the two CuproptosisClusters. (C) The overall survival of the three gene subtypes (Log-rank test). (D) The heatmap for the connections between clinicopathologic features and the three gene subtypes (blue: low expression; red: high expression). (E) The differences in the expression of 12 cuproptosis-related genes among the three gene subtypes (Wilcoxon test). DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. *p < 0.05, **p < 0.01, ***p < 0.001.
Development of the cuproptosis-related gene signature and CRG_score in the BRCA training cohort
To further investigate the underlying mechanisms regulating breast cancer progression, we constructed a cuproptosis-related gene signature by performing the LASSO regression analysis based on the 25 prognostic subtype-related genes. And the Cuproptosis-Related Gene Score (CRG_score) was calculated as follows: CRG_score = (0.00523335734652904 * PGK1) + (0.0185186220360186 * SLC52A2) − (0.0261649532232335 * SEC14L2) + (0.0129284730002406 * RAD23B) − (0.0297061289632435 * SLC16A6) − (0.0142610013224485 * CCL5) + (0.0014557255170858 * MAL2). To better estimate the characteristics of patients with different levels of CRG_score, we divided the BRCA training cohort into high- and low-CRG_score groups depending on the median CRG_score. The alluvial diagram depicted the distribution of BRCA patients within the two cuproptosis clusters, three gene subtypes, and two CRG score groups (Figure 6A). As shown in Figure 6B, gene subtype C had significantly higher CRG_scores than the other two gene subtypes. Moreover, cuproptosis cluster B had higher CRG_scores than cluster A (Figure 6C). Figure 6D also showed the distribution plot of the survival of each BRCA patient from the training cohort, which indicated a higher death probability in the high-CRG_score group. The Kaplan-Meier curve consistently suggested a worse prognosis for the high-CRG_score group than the low-CRG_score group (p < 0.001; Figure 6E). Moreover, we performed the time-dependent receiver operating characteristic (ROC) analysis to calculate the AUC values of this cuproptosis-related gene signature (0.741 for 3-year, 0.707 for 5-year, and 0.716 for 10-year; Figure 6F). We also conducted a nomogram plot analysis, which suggested an excellent advantage for CRG_score in long-term survival prediction (Figure 6G). Finally, the calibration chart revealed an excellent performance of the CRG_score among the predicted and observed overall survival with a C index of 0.69 (Figure 6H).
Figure 6 Development of the Cuproptosis-Related Gene Signature and CRG_score in the BRCA Training cohort. (A)The alluvial diagram showing the connection between CuproptosisClusters, gene subtypes, and CRG_score. (B) The level of CRG_score in the three gene subtypes (Wilcoxon test). (C) The level of CRG_score in the two CuproptosisClusters (Wilcoxon test). (D) The ranked dot and scatter plots of CRG_score distribution and patient survival status in BRCA training cohort. (E) The overall survival of the high and low CRG_score groups in BRCA training cohort (Log-rank test). (F) The ROC curves for the predictive efficiency of the CRG_score in BRCA training cohort (green: 3 year; blue: 5 year; red: 10 year). (G) The Nomogram to predict 3-, 5- and 10-year OS in the BRCA training cohort. (H) The Calibration plots of the nomogram to predict OS at 3-, 5- and 10-year (green: 3 year; blue: 5 year; red: 10 year). CRG_score, Cuproptosis-Related Gene Score; ROC, receiver operating characteristic. *p < 0.05, **p < 0.01, ***p < 0.001.
Validation of the cuproptosis-related gene signature in BRCA testing cohorts
To validate the reliability and reproducibility of our cuproptosis-related gene signature, we calculated CRG_scores of three independent external validation BRCA groups, including GSE7390 with 198 BRCA patients, GSE58812 with 107 BRCA patients, and GSE42568 with 104 BRCA patients. Then the patients in each testing cohort were classified into high- and low-CRG_score groups based on the median CRG_score value. Patients with high CRG_scores in all three testing cohorts showed worse survival status (Figures 7A–C). Similarly, survival analysis revealed a significantly better overall survival of patients with low CRG_scores than those with high CRG_scores (p = 0.001, 0.038, 0.003, respectively; Figure 7D–F). Moreover, the high AUC values also suggested an excellent ability of the CRG_score to predict the long-term prognosis of BRCA patients in testing cohorts (Figures 7G–I). Thus, the results above showed a similar tendency to the training cohort, indicating the cuproptosis-related gene signature was stable and reliable.
Figure 7 Validation of the Cuproptosis-Related Gene Signature in BRCA Testing Cohorts. BRCA Testing Cohorts: GSE7390 (n = 198), GSE58812 (n = 107), and GSE42568 (n = 104). (A–C) The ranked dot and scatter plots of CRG_score distribution and patient survival status in BRCA testing cohorts (GSE7390, GSE58812, and GSE42568 sets, respectively). (D–F) The overall survival of the high and low CRG_score groups in BRCA testing cohorts (GSE7390, GSE58812, and GSE42568 sets, respectively). (G–I) The ROC curves for the predictive efficiency of the CRG_score in BRCA training cohort (red: 3 year; blue: 5 year; yellow: 8 year; green: 10 year).
Moreover, to further clarify the applicability of different molecular subtypes of breast cancer to our predictive gene signature, the PAM50 algorithm was conducted on the testing cohorts (GSE7390, GSE58812, and GSE42568; n = 409). The distribution of molecular subtypes according to the PAM50 signature was as follows: 13 normal-like (3%), 126 basal-like (31%), 111 luminal A (27%), 102 luminal B (25%), and 57 HER2+ (14%). We also calculated the CRG_scores of each PAM50 subtype and divided them into high- and low-CRG_score groups based on the median CRG_score value, respectively. The survival and ROC curve analyses were also performed for each PAM50 subtype (Figure S3). Patients of luminal A and HER2+ subtypes showed a better prognosis in the low CRG_score group (Figures S3C, E), while the survival of other subtypes was not statistically significant. In addition, the AUC values for all of the PAM50 subtypes were high, which indicated that our CRG_score could predict the long-term prognosis of BRCA patients with different molecular subtypes.
Relationship of clinical characteristics and TME characteristics with CRG_score in the BRCA training cohort
Further, we investigated the clinical characteristics of the training cohort’s high- and low-CRG_score groups. For ages, patients aged ≥65 had a higher CRG_score, and the high CRG_score of patients of different ages was positively related to poor prognosis (Figure 8A). Then a relatively higher CRG_score was observed in patients with T1-2 and a lower CRG_score in patients with T3-4. In addition, the low CRG_scores of patients with different T stages suggested a better overall survival than the high CRG_score group (Figure 8B). Similarly, patients with N1-3 were correlated with a relatively high CRG_score, while the low CRG_score of patients with different N stages both lived longer than those in the high-CRG_score group (Figure 8C). Therefore, patients with a high CRG_score of types ≥65 age, stage T1-2, and stage N1-3 had worse long-term survival.
Figure 8 Relationship of Clinical Characteristics and TME Characteristics with CRG_score in the BRCA Training Cohort. The relationship between age (A), T (B), N (C) and CRG_score. The TNM system of cancer staging reflects the extent of tumor growth, where primary tumor (T), and nodal status for metastasis (N) (Wilcoxon test for boxplot; Log-rank test for survival analysis). (D) Survival analysis of ImmuneScore, StromalScore, ESTIMATEScore, and TumorPurity in BRCA patients (Log-rank test). (E) The TME score of the two CuproptosisClusters in the high and low CRG_score groups (Spearman test). (F) The abundance of 23 infiltrating immune cell types in the high and low CRG_score groups (Spearman test). *p < 0.05, **p < 0.01, ***p < 0.001.
To analyze the relationship between TME and CRG_score, we calculated the TME score (ImmuneScore, StromalScore, ESTIMATEScore, and TumorPurity) using the “estimate” package (Figure 8D). We divided patients into low- and high-TME score groups depending on the median values. We found that low ImmuneScore and ESTIMATEScore were significantly correlated with poor overall survival (p < 0.001, p = 0.025, respectively), while high TumorPurity was significantly correlated with worse overall survival (p=0.025). Moreover, there was no significant correlation between the prognosis of BRCA patients and StromalScore. Besides, a low CRG_score was also significantly associated with ImmuneScore, StromalScore, and ESTIMATEScore (Figure 8E). Furthermore, the relationship between the cuproptosis-related gene signature and immune cell abundance was further estimated in box plots and scatter diagrams (Figure 8F). Here, the results above revealed a close relationship between the cuproptosis-related gene prognostic model and TME, which suggested the possibility of immunotherapy targeting cuproptosis-related genes to block BRCA progression.
Prediction of immunotherapy response in high- and low-CRG_score groups in BRCA
Since cuproptosis played an essential role in the TME, we further investigated its influence on immunotherapy for BRCA. First, the expression of 32 critical immune checkpoints in high- and low-CRG_score groups was examined (Figure 9A). We observed that most immune checkpoints were significantly overexpressed in the low-CRG_score group. Next, we investigated the application of CRG_score in the therapy of BRCA. The BRCA immunotherapy profile of patients from the TCIA database revealed that the low-CRG_score group had higher ICI scores than the high-CRG_score group and was more responsive to the immunotherapy than the high-CRG_score group (Figure 9B). In addition, we assessed the responses of the high- and low-CRG_score groups to conventional and novel chemotherapeutic agents. The high-CRG_score group was more sensitive to lapatinib, nilotinib, pazopanib, metformin, lenalidomide, camptothecin, cytarabine, bexarotene, midostaurin, shikonin, temsirolimus, and vorinostat. In contrast, the low-CRG_score group was more responsive to paclitaxel, imatinib, sorafenib, and rapamycin (Figure 9C). Therefore, BRCA patients with high-CRG_score were characterized by abundant immune infiltration, high expression of immune checkpoints, and better response to immunotherapy and chemotherapy.
Figure 9 Prediction of Immunotherapy Response in High- and Low-CRG_score groups in BRCA. (A) The expression levels of 32 immune checkpoints in the high and low CRG_score groups (Wilcoxon test). (B) The immunotherapy response between the high and low CRG_score groups (Wilcoxon test). (C) The relationships between CRG_score and chemotherapeutic sensitivity (Wilcoxon test). *p < 0.05, **p < 0.01, ***p < 0.001.
Expressions pattern of cuproptosis-related gene signature in BRCA
To further explore the important role of the cuproptosis-related gene signature in breast cancer, we analyzed their protein expression patterns in the normal and tumor samples in the Human Protein Atlas (HPA) database (Figure 10A). RAD23B protein was strongly expressed in breast cancer tumor tissues, SLC52A2 was moderately expressed, and SEC14L2 was weakly expressed. Additionally, PGK1 was not expressed in most breast cancers, some were weakly expressed, and few were moderately expressed. As for SLC16A6, CCL5, and MAL2, their protein expression was almost negative in breast cancer tissues based on the results of the HPA database. According to the previous study (8), we also examined the transcriptional expression of our cuproptosis-related gene signature in the cuproptosis cell model. As shown in Figure 10B, RAD23B was significantly down-regulated in SUM159 and MCF7 cells after cuproptosis induction; SLC16A6 was significantly down-regulated in MCF7 cells, while other genes showed no significant change.
Figure 10 Expressions of Cuproptosis-Related Gene Signature in BRCA. (A) Representative IHC images of PGK1, SLC52A2, SEC14L2, RAD23B, SLC16A6, CCL5, and MAL2 across clinical specimens of normal and tumor samples in the Human Protein Atlas (HPA) database. Bar charts represent IHC staining intensities of PGK1 (11 patients), SLC52A2 (10 patients), SEC14L2 (10 patients), RAD23B (12 patients), SLC16A6 (10 patients), CCL5 (11 patients), and MAL2 (11 patients). (B) The relative mRNA expression of PGK1, SLC52A2, SEC14L2, RAD23B, SLC16A6, CCL5, and MAL2 in the cuproptosis cell model (RT-PCR). *p < 0.05, **p < 0.01.
We also analyzed the mRNA expression of the cuproptosis-related gene signature in breast cancer tissues, suggesting a higher transcriptional expression of all 7 cuproptosis-related genes in BRCA tissues than in adjacent normal tissues (Figure S4A). Furthermore, overall survival analysis showed a worse prognosis for patients with high expression of PGK1, SLC52A2, RAD23B, and MAL2; while patients with up-regulated SEC14L2, SLC16A6, and CCL5 could live longer (Figure S4B).
The role of RAD23B in the regulation of breast cancer progression, chemotherapy and immunotherapy in vitro
The protein and mRNA expression levels of RAD23B were significantly higher in breast cancer tissues than in paraneoplastic tissues, and its expression was positively correlated with poor prognosis in BRCA patients. Moreover, the expression of RAD23B was decreased in the cuproptosis cell model, which indicated that it acted as a “risk” factor and antagonized cuproptosis in BRCA progression. Therefore, RAD23B was selected as a promising target for in-depth experimental validation. As shown in Figure 11A, the transcriptional expression of RAD23B was examined in breast epithelial cell lines (MCF10A), luminal breast cancer cells (MCF7), and triple-negative breast cancer cells (MDA-MB-231, SUM-159, and BT549). Due to the high expression of RAD23B in SUM-159 and MCF7 cells, two RAD23B siRNA (small interfering RNA) constructs were used to knockdown its expression in these two cell lines (Figure 11B). After the reduction of RAD23B, cell viability was significantly inhibited (Figure 11C). The MTT assay reduced the IC50 of Paclitaxel in the RAD23B down-regulated groups compared with that of the control group (Figure 11D). To further explore the expression pattern of RAD23B in breast cancer, we examined the mRNA expression level of RAD23B in 34 pairs of breast cancer tissues and adjacent normal breast tissues. It suggested that RAD23B was significantly overexpressed in breast cancer tissues than in adjacent normal breast tissues, and its expression was positively correlated with pathological grade (Figure 11E, F). Finally, we also tested the mRNA expression of PD1 and PDL1 in breast cancer tissues by performing real-time RT-PCR. Then the correlation between RAD23B and PD1/PDL1 was evaluated (Figure 11F; r = 0.774, p < 0.001; r = 0.577, p < 0.001, respectively).
Figure 11 The Role of RAD23B in the Regulation of the Breast Cancer Progression, Chemotherapy and Immunotherapy in Vitro. (A) The mRNA expression levels of RAD23B in different breast cancer cell lines (RT-PCR). (B) The relative mRNA expression of RAD23B in SUM-159 and MCF7 breast cancer cells (RT-PCR; One-way ANOVA). (C) The cell viability assay after RAD23B reduction in SUM-159 and MCF7 breast cancer cells (MTT; Two-way ANOVA). (D) The IC50 of Paclitaxel after RAD23B reduction in SUM-159 and MCF7 breast cancer cells (MTT). (E) The relative mRNA expression of RAD23B in 34 pairs of breast cancer tissues and adjacent normal breast tissues (RT-PCR; Student’s t-test). (F) The relative mRNA expression of RAD23B in different pathologic stages of breast cancer tissues (RT-PCR; Student’s t-test). (G) Correlation between RAD23B mRNA expression and PD1/PDL1 mRNA expression (RT-PCR; Spearman test). *p < 0.05, **p < 0.01, ***p < 0.001.
Discussion
Because of breast cancer’s high incidence and recurrence rates, its treatment has been an ongoing challenge for decades (4). Numerous studies have attempted to determine the significance of the immune microenvironment in breast cancer progression, and immunotherapy may be a viable treatment option for breast cancer patients (42–44). However, immunotherapy for breast cancer patients is still in its infancy, and more investigation is needed to benefit more people.
As an essential cofactor of key enzymes, copper must maintain a dynamic low concentration to maintain normal physiological activity. A few studies have noted the role of cooper in regulating cancer progression. It has been reported that the serum level of copper was significantly increased in the BRCA group compared to the control group, indicating its function in the early detection and monitoring of breast cancer (13, 24, 25). In triple-negative breast cancer (TNBC), inhibition of mitochondrial copper shifted tumor cells from respiration to glycolysis to reduce energy production, ultimately inhibiting tumor growth and improving prognosis (26–28). Recently reported, an excess load of intracellular copper could induce cell death termed cuproptosis (8). Independent of known cell death pathways, cuproptosis does not activate caspase-3 and cannot be blocked by apoptotic inhibitors (8). This study mainly discussed 12 cuproptosis-related genes, including FDX1, LIPT1, LIAS, DLD, DBT, GCSH, DLST, DLAT, and PDHA1, PDHB, SLC31A1, ATP7A, and ATP7B. These genes are primarily involved in processes such as glycolysis and the tricarboxylic acid (TCA) cycle (45, 46), steroids (47, 48), and vitamin D metabolism (49). SCL31A1 as a copper importer and ATP7A and ATP7B as copper exporters are essential for maintaining intracellular copper concentration (50, 51). As Peter et al. suggested, overexpression of SCL31A1 and deletion of ATP7B may increase susceptibility to cuproptosis (8). In addition, the knockout of nine genes (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, GCSH, and DBT) conferred resistance to cuproptosis (8). Although various studies have highlighted the significance of copper in breast cancer, there is a dearth of studies on the association between cuproptosis and breast cancer, mainly its function in the immune microenvironment and immunotherapy of breast cancer.
Our study first summarized the expression and mutation patterns of CRGs based on TCGA-BRCA and GSE20685. Although the frequency of global alterations was only 2.94%, all CRGs showed significant differences in expression and prognosis in BRCA samples compared to normal samples. Performing the unsupervised clustering algorithm, we classified breast cancer patients into two cuproptosis patterns (CuproptosisCluster A and CuproptosisCluster B). Compared to patients with Cluster A, those with Cluster B showed more advanced clinicopathological characteristics and a worse OS. A global examination of the TME for both clusters revealed that Cluster B presented enrichment in most immune cells and important immune checkpoints. In addition, Cluster A showed a higher stromal score with significance, while Cluster B had a higher immune score. Thus, the results above indicated that these two clusters were closely associated with TME in BRCA, implying a crucial role of CRGs in the immune regulation of breast cancer. Next, three cuproptosis-related gene subtypes were identified according to the DEGs of the two CuproptosisClusters. To further explore the role of cuproptosis in BRCA progression and TME, a cuproptosis-related gene signature (PGK1, SLC52A2, SEC14L2, RAD23B, SLC16A6, CCL5, and MAL2) and the CRG_score were constructed based on the training cohort and validated in the testing cohorts, as well as PAM50 subtyped testing cohorts. Furthermore, patients in high-CRG_score groups showed a worse prognosis under different clinicopathological features. And patients with a low-CRG_score exhibited enrichment in most immune cells and important immune checkpoints, and are more sensitive to immunotherapy. Also, a quantitative nomogram depending on CRG_score and tumor stage facilitated the prognostic stratification of breast cancer patients, further promoting the clinical application of CRG_score. In summary, the characteristics of TME differed significantly in the high- and low-CRG_score groups, suggesting CRGs could provide reasonable recommendations for personalized immunotherapy for breast cancer patients.
Regarding the cuproptosis-related gene signature, previous studies have identified the crucial roles of these seven genes in cell metabolism. PGK1, an ATP-generating enzyme, mediates mitochondrial metabolism and promotes tumorigenesis (52). SLC family genes such as SLC52A2 and SLC16A6 are important transporters in metabolic processes, and their dysregulation is associated with various diseases (53–56). SEC14L2 encodes lipid binding proteins and facilitates the uptake of Vitamin E (57). RAD23B is involved in nucleotide excision repair (NER) and is associated with cell apoptosis (58). CCL5, expressed and secreted by activated and normal T cells, could regulate the migration and chemotaxis of inflammatory cells (59–61). MAL2 has been reported to work as an essential component of the machinery for transcytosis in hepatoma HepG2 cells (62). However, there is still no research focus on the relationship between these seven genes and breast cancer cuproptosis. We then analyzed the protein expression of these seven genes in breast cancer tissues from the HPA database and their mRNA expression in the cuproptosis cell model. In particular, RAD23B was screened out for in vitro experimental validation. The results of RT-PCR, cell viability, and the IC50 assay illustrated that RAD23B expression was positively correlated with breast cancer progression, drug resistance, and poor prognosis in BRCA patients. More importantly, both PD1 and PDL1 were positively correlated with RAD23B, suggesting that patients with up-regulated RAD23B were more sensitive to immune checkpoint-blocking therapy targeting the PD-1/PD-L1 axis. Thus, our results confirmed the important role of cuproptosis in TME and immunotherapy for breast cancer, providing new ideas for immunotherapy with blocked immune checkpoints. Similarly, our cuproptosis-related scoring system was of great utility for clinical patient stratification, predicting the efficacy of adjuvant chemotherapy and patient prognosis. Also, our study laid an important research foundation for the role of cuproptosis in controlling the progression of breast cancer and made it easier to study its molecular mechanisms in more depth in the future.
Although we have performed a comprehensive analysis of cuproptosis in breast cancer and screened out potential targets to lay the foundation for future exploration of breast cancer progression, this study still has some limitations. Since our breast cancer samples are only obtained from retrospective studies from the TCGA and GEO databases, more cases from prospective research are required. In addition, experimental studies in vivo and in vitro are needed to validate our findings. Furthermore, additional research is necessary to identify the specific molecular mechanisms of cuproptosis regulating breast cancer progression.
Data availability statement
The data sets analyzed during the current study are available in the TCGA (https://portal.gdc.cancer.gov/), accession numbers TCGA-BRCA, BRCA-FPKM; GEO repository (https://www.ncbi.nlm.nih.gov/geo/), accession numbers GSE20685, GSE7390, GSE58812, and GSE42568.
Ethics statement
The studies involving human participants were reviewed and approved by The Ethics Committee on Human Research of the First Affiliated Hospital of Xi’an Jiaotong University. The patients/participants provided their written informed consent to participate in this study.
Author contributions
SS designed the study and wrote the manuscript. MZ and PX analyzed data and contributed to writing the manuscript. SW and YW revised the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This research was funded by the National Natural Science Foundation of China (No. 82173365), the Key projects of Natural Science Foundation of Shaanxi Province(2021JZ-36); Natural Science Basic Research Project of Shaanxi Province (2020JM-362); the Institutional foundation of The First Affiliated Hospital of Xi’an Jiaotong University (No.2020QN-34).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.978909/full#supplementary-material
Supplementary Figure 1 | PCA analysis for the construction of CuproptosisClusters.
Supplementary Figure 2 | PCA analysis for the construction of Gene Subtypes.
Supplementary Figure 3 | Validation of the cuproptosis-related gene signature in different BRCA subtypes of testing cohorts.
Supplementary Figure 4 | The survival analysis and expression of Cuproptosis-related gene signature in BRCA.
References
1. Guo F, Scholl M, Fuchs EL, Berenson AB, Kuo YF. BRCA testing in unaffected young women in the united states, 2006-2017. Cancer (2020) 126(2):337–43. doi: 10.1002/cncr.32536
2. Britt KL, Cuzick J, Phillips KA. Key steps for effective breast cancer prevention. Nat Rev Cancer (2020) 20(8):417–36. doi: 10.1038/s41568-020-0266-x
3. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA Cancer J Clin (2022) 72(1):7–33. doi: 10.3322/caac.21708
4. Harbeck N, Penault-Llorca F, Cortes J, Gnant M, Houssami N, Poortmans P, et al. Breast cancer. Nat Rev Dis Primers (2019) 5(1):66. doi: 10.1038/s41572-019-0111-2
5. Deepak KGK, Vempati R, Nagaraju GP, Dasari VR, Nagini S, Rao DN, et al. Tumor microenvironment: Challenges and opportunities in targeting metastasis of triple negative breast cancer. Pharmacol Res (2020) 153:104683. doi: 10.1016/j.phrs.2020.104683
6. Lim B, Woodward WA, Wang X, Reuben JM, Ueno NT. Inflammatory breast cancer biology: the tumour microenvironment is key. Nat Rev Cancer (2018) 18(8):485–99. doi: 10.1038/s41568-018-0010-y
7. Choi J, Gyamfi J, Jang H, Koo JS. The role of tumor-associated macrophage in breast cancer biology. Histol Histopathol (2018) 33(2):133–45. doi: 10.14670/HH11-916
8. Tsvetkov P, Coy S, Petrova B, Dreishpoon M, Verma A, Abdusamad M, et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science (2022) 375(6586):1254–61. doi: 10.1126/science.abf0529
9. Bock FJ, Tait SWG. Mitochondria as multifaceted regulators of cell death. Nat Rev Mol Cell Biol (2020) 21(2):85–100. doi: 10.1038/s41580-019-0173-8
10. Bandmann O, Weiss KH, Kaler SG. Wilson's disease and other neurological copper disorders. Lancet Neurol (2015) 14(1):103–13. doi: 10.1016/S1474-4422(14)70190-5
11. Bolea R, Hortells P, Martin-Burriel I, Vargas A, Ryffel B, Monzon M, et al. Consequences of dietary manganese and copper imbalance on neuronal apoptosis in a murine model of scrapie. Neuropathol Appl Neurobiol (2010) 36(4):300–11. doi: 10.1111/j.1365-2990.2010.01065.x
12. Steinbrueck A, Sedgwick AC, Brewster JT 2nd, Yan KC, Shang Y, Knoll DM, et al. Transition metal chelators, pro-chelators, and ionophores as small molecule cancer chemotherapeutic agents. Chem Soc Rev (2020) 49(12):3726–47. doi: 10.1039/C9CS00373H
13. Pavithra V, Sathisha TG, Kasturi K, Mallika DS, Amos SJ. Ragunatha s. serum levels of metal ions in female patients with breast cancer. J Clin Diagn Res (2015) 9(1):BC25–c7. doi: 10.7860/JCDR/2015/11627.5476
14. Diez M, Cerdan FJ, Arroyo M, Balibrea JL. Use of the copper/zinc ratio in the diagnosis of lung cancer. Cancer (1989) 63(4):726–30. doi: 10.1002/1097-0142(19890215)63:4<726::AID-CNCR2820630421>3.0.CO;2-P
15. Brady DC, Crowe MS, Turski ML, Hobbs GA, Yao X, Chaikuad A, et al. Copper is required for oncogenic BRAF signalling and tumorigenesis. Nature (2014) 509(7501):492–6. doi: 10.1038/nature13180
16. Gambling L, Kennedy C, McArdle HJ. Iron and copper in fetal development. Semin Cell Dev Biol (2011) 22(6):637–44. doi: 10.1016/j.semcdb.2011.08.011
17. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170
18. Wang S, Su W, Zhong C, Yang T, Chen W, Chen G, et al. An eight-CircRNA assessment model for predicting biochemical recurrence in prostate cancer. Front Cell Dev Biol (2020) 8:599494. doi: 10.3389/fcell.2020.599494
19. Zeng L, Fan X, Wang X, Deng H, Zhang K, Zhang X, et al. Bioinformatics analysis based on multiple databases identifies hub genes associated with hepatocellular carcinoma. Curr Genomics (2019) 20(5):349–61. doi: 10.2174/1389202920666191011092410
20. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7
21. Ferreira MR, Santos GA, Biagi CA, Silva Junior WA, Zambuzzi WF. GSVA score reveals molecular signatures from transcriptomes for biomaterials comparison. J BioMed Mater Res A (2021) 109(6):1004–14. doi: 10.1002/jbm.a.37090
22. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337
23. Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell (2015) 160(1-2):48–61. doi: 10.1016/j.cell.2014.12.033
24. Shen Y, Peng X, Shen C. Identification and validation of immune-related lncRNA prognostic signature for breast cancer. Genomics (2020) 112(3):2640–6. doi: 10.1016/j.ygeno.2020.02.015
25. Huang ZD, Lin LL, Liu ZZ, Hu C, Gu HY, Wei RX. m6A modification patterns with distinct immunity, metabolism, and stemness characteristics in soft tissue sarcoma. Front Immunol (2021) 12:765723. doi: 10.3389/fimmu.2021.765723
26. Wu J, Zhu Y, Luo M, Li L. Comprehensive analysis of pyroptosis-related genes and tumor microenvironment infiltration characterization in breast cancer. Front Immunol (2021) 12:748221. doi: 10.3389/fimmu.2021.748221
27. Hu FF, Liu CJ, Liu LL, Zhang Q, Guo AY. Expression profile of immune checkpoint genes and their roles in predicting immunotherapy response. Brief Bioinform (2021) 22(3):bbaa176. doi: 10.1093/bib/bbaa176
28. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007
29. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) (2021) 2(3):100141. doi: 10.1016/j.xinn.2021.100141
30. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an r package for comparing biological themes among gene clusters. OMICS (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118
31. Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst (2015) 1(6):417–25. doi: 10.1016/j.cels.2015.12.004
32. Liu J, Zhou S, Li S, Jiang Y, Wan Y, Ma X, et al. Eleven genes associated with progression and prognosis of endometrial cancer (EC) identified by comprehensive bioinformatics analysis. Cancer Cell Int (2019) 19:136. doi: 10.1186/s12935-019-0859-1
33. Tibshirani R. The lasso method for variable selection in the cox model. Stat Med (1997) 16(4):385–95. doi: 10.1002/(SICI)1097-0258(19970228)16:4<385::AID-SIM380>3.0.CO;2-3
34. Engebretsen S, Bohlin J. Statistical predictions with glmnet. Clin Epigenetics (2019) 11(1):123. doi: 10.1186/s13148-019-0730-1
35. Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med (2013) 32(30):5381–97. doi: 10.1002/sim.5958
36. Parker JS, Mullins M, Cheang MC, Leung S, Voduc D, Vickery T, et al. Supervised risk predictor of breast cancer based on intrinsic subtypes. J Clin Oncol (2009) 27(8):1160–7. doi: 10.1200/JCO.2008.18.1370
37. Gendoo DM, Ratanasirigulchai N, Schroder MS, Pare L, Parker JS, Prat A, et al. Genefu: an R/Bioconductor package for computation of gene expression-based signatures in breast cancer. Bioinformatics (2016) 32(7):1097–9. doi: 10.1093/bioinformatics/btv693
38. Uhlen M, Fagerberg L, Hallstrom BM, Lindskog C, Oksvold P, Mardinoglu A, et al. Proteomics. tissue-based map of the human proteome. Science (2015) 347(6220):1260419. doi: 10.1126/science.1260419
39. Thul PJ, Lindskog C. The human protein atlas: A spatial map of the human proteome. Protein Sci (2018) 27(1):233–44. doi: 10.1002/pro.3307
40. Song S, Liu J, Zhang M, Gao X, Sun W, Liu P, et al. Eukaryotic translation initiation factor 3 subunit b could serve as a potential prognostic predictor for breast cancer. Bioengineered (2022) 13(2):2762–76. doi: 10.1080/21655979.2021.2017567
41. Hochberg Y, Benjamini Y. More powerful procedures for multiple significance testing. Stat Med (1990) 9(7):811–8. doi: 10.1002/sim.4780090710
42. Xu Q, Chen S, Hu Y, Huang W. Landscape of immune microenvironment under immune cell infiltration pattern in breast cancer. Front Immunol (2021) 12:711433. doi: 10.3389/fimmu.2021.711433
43. Wang S, Xiong Y, Zhang Q, Su D, Yu C, Cao Y, et al. Clinical significance and immunogenomic landscape analyses of the immune cell signature based prognostic model for patients with breast cancer. Brief Bioinform (2021) 22(4):bbaa311. doi: 10.1093/bib/bbaa311
44. Mehta AK, Kadel S, Townsend MG, Oliwa M, Guerriero JL. Macrophage biology and mechanisms of immune suppression in breast cancer. Front Immunol (2021) 12:643771. doi: 10.3389/fimmu.2021.643771
45. Xu S, Liu CX, Xu W, Huang L, Zhao JY, Zhao SM. Butyrate induces apoptosis by activating PDC and inhibiting complex I through SIRT3 inactivation. Signal Transduct Target Ther (2017) 2:16035. doi: 10.1038/sigtrans.2016.35
46. He K, Chen Z, Ma Y, Pan Y. Identification of high-copper-responsive target pathways in Atp7b knockout mouse liver by GSEA on microarray data sets. Mamm Genome (2011) 22(11-12):703–13. doi: 10.1007/s00335-011-9359-x
47. Sheftel AD, Stehling O, Pierik AJ, Elsasser HP, Muhlenhoff U, Webert H, et al. Humans possess two mitochondrial ferredoxins, Fdx1 and Fdx2, with distinct roles in steroidogenesis, heme, and Fe/S cluster biosynthesis. Proc Natl Acad Sci U S A (2010) 107(26):11775–80. doi: 10.1073/pnas.1004250107
48. Strushkevich N, MacKenzie F, Cherkesova T, Grabovec I, Usanov S, Park HW. Structural basis for pregnenolone biosynthesis by the mitochondrial monooxygenase system. Proc Natl Acad Sci U S A (2011) 108(25):10139–43. doi: 10.1073/pnas.1019441108
49. Zalewski A, Ma NS, Legeza B, Renthal N, Fluck CE, Pandey AV. Vitamin d-dependent rickets type 1 caused by mutations in CYP27B1 affecting protein interactions with adrenodoxin. J Clin Endocrinol Metab (2016) 101(9):3409–18. doi: 10.1210/jc.2016-2124
50. Lutsenko S. Human copper homeostasis: a network of interconnected pathways. Curr Opin Chem Biol (2010) 14(2):211–7. doi: 10.1016/j.cbpa.2010.01.003
51. Moriya M, Ho YH, Grana A, Nguyen L, Alvarez A, Jamil R, et al. Copper is taken up efficiently from albumin and alpha2-macroglobulin by cultured human cells by more than one mechanism. Am J Physiol Cell Physiol (2008) 295(3):C708–21. doi: 10.1152/ajpcell.00029.2008
52. Lu Z, Hunter T. Metabolic kinases moonlighting as protein kinases. Trends Biochem Sci (2018) 43(4):301–10. doi: 10.1016/j.tibs.2018.01.006
53. Console L, Tolomeo M, Colella M, Barile M, Indiveri C. Reconstitution in proteoliposomes of the recombinant human riboflavin transporter 2 (SLC52A2) overexpressed in e. coli Int J Mol Sci (2019) 20(18):4416. doi: 10.3390/ijms20184416
54. Zhang L, Li M, Cui Z, Chai D, Guan Y, Chen C, et al. Systematic analysis of the role of SLC52A2 in multiple human cancers. Cancer Cell Int (2022) 22(1):8. doi: 10.1186/s12935-021-02432-7
55. Halestrap AP. The SLC16 gene family - structure, role and regulation in health and disease. Mol Aspects Med (2013) 34(2-3):337–49. doi: 10.1016/j.mam.2012.05.003
56. Higuchi K, Sugiyama K, Tomabechi R, Kishimoto H, Inoue K. Mammalian monocarboxylate transporter 7 (MCT7/Slc16a6) is a novel facilitative taurine transporter. J Biol Chem (2022) 298(4):101800. doi: 10.1016/j.jbc.2022.101800
57. Wright ME, Peters U, Gunter MJ, Moore SC, Lawson KA, Yeager M, et al. Association of variants in two vitamin e transport genes with circulating vitamin e concentrations and prostate cancer risk. Cancer Res (2009) 69(4):1429–38. doi: 10.1158/0008-5472.CAN-08-2343
58. Friedberg EC. How nucleotide excision repair protects against cancer. Nat Rev Cancer (2001) 1(1):22–33. doi: 10.1038/35094000
59. Schall TJ, Bacon K, Toy KJ, Goeddel DV. Selective attraction of monocytes and T lymphocytes of the memory phenotype by cytokine RANTES. Nat (1990) 347(6294):669–71. doi: 10.1038/347669a0
60. Tanaka T, Bai Z, Srinoulprasert Y, Yang BG, Hayasaka H, Miyasaka M. Chemokines in tumor progression and metastasis. Cancer Sci (2005) 96(6):317–22. doi: 10.1111/j.1349-7006.2005.00059.x
61. Romagnani S. Cytokines and chemoattractants in allergic inflammation. Mol Immunol (2002) 38(12-13):881–5. doi: 10.1016/S0161-5890(02)00013-5
Keywords: cuproptosis, breast cancer, cancer molecular subtypes, tumor microenvironment, immunotherapy
Citation: Song S, Zhang M, Xie P, Wang S and Wang Y (2022) Comprehensive analysis of cuproptosis-related genes and tumor microenvironment infiltration characterization in breast cancer. Front. Immunol. 13:978909. doi: 10.3389/fimmu.2022.978909
Received: 27 June 2022; Accepted: 06 October 2022;
Published: 20 October 2022.
Edited by:
Marta Truffi, Scientific Clinical Institute Maugeri (ICS Maugeri), ItalyReviewed by:
Jamshid Hadjati, Tehran University of Medical Sciences, IranLinlin Guo, The Ohio State University, United States
Copyright © 2022 Song, Zhang, Xie, Wang and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yaochun Wang, wangyc@xjtufh.edu.cn; Shuhong Wang, wangsh@xjtufh.edu.cn
†These authors have contributed equally to this work