A novel cuproptosis-associated immune risk model for prediction of prognosis and response to immunotherapy in triple-negative breast cancer

Cuproptosis regulated patterns in TNBC

The study design is presented in the flow chart (Figure 1). We observed differential expression of 12 cuproptosis-related genes between TNBC samples and paired normal tissues, indicating dysregulation of these genes in TNBC (Figure 2A). We further explored the relationship among these cuproptosis-related genes and found potential interactions between each gene (Figure 2B). Inspired by above results, we hypothesized that TNBC might present different regulated patterns according to cuproptosis-related gene expression, and then performed a consensus clustering analysis with TNBC samples. The TNBC patients were divided in two clusters (k value = 2), named cuproptosis cluster A and cluster B (Figure 2C, 2D). As expected, the differential expression of cuproptosis-related genes between the two clusters indicated distinct cuproptosis patterns in TNBC (Figure 2E). Importantly, the overall survival analysis revealed the difference prognosis in two cuproptosis-related clusters, suggesting that different cuproptosis patterns were associated with the clinical outcome of TNBC patients (Figure 2F).

Figure 1. Flow chart of our study.

Figure 2. Identification of different cuproptosis patterns in TNBC. (A) The expression of twelve cuproptosis-related genes in normal and TNBC tissues from TCGA data. (B) The correlation network between twelve cuproptosis genes in TNBC. (C) Cumulative distribution function (CDF) and relative change in the area under the CDF curve (CDF Delta area) of consensus clustering based on cuproptosis-related genes. (D) Unsupervised consensus clustering (k = 2) divides TNBC samples into cluster A and cluster B based on 12 cuproptosis-related genes. (E) The different expression of twelve cuproptosis genes between two clusters displayed with a heatmap. (F) KM curves of overall survival for the two cuproptosis clusters. (* P<0.05, **P< 0.01).

Identification of cuproptosis and immune-related genes

To further investigate the effect of cuproptosis in TNBC, we firstly identified the differentially expressed genes (DEGs) between cluster A and cluster B (Figure 3A and Supplementary Table 1). Subsequently, through a univariate Cox and intersection analysis, 207 prognostic DEGs were used to investigate the association between cuproptosis and biological characteristics of TNBC (Figure 3B). The Gene Ontology (GO) and HALLMARK analysis revealed that the genes were significantly enriched in HALLMARK_INTERFERON_GAMMA_RESPONSE and GOBP_IMMUNE_RESPONSE (Supplementary Table 2). The results above indicated the potential involvement of immune-related mechanisms in cuproptosis (Figure 3C).

Figure 3. Identification of cuproptosis and immune-related genes in TNBC. (A) Volcano plot exhibiting the differentially expressed genes (DEGs) between the two cuproptosis clusters. (B) Screening 207 prognosis DEGs through the intersection of DEGs and prognostic genes in TNBC. (C) Enrichment analysis with GO and HALLMARK pathways based on 207 prognosis DEGs. (D) Identification of 20 prognostic DEGs associated with GO_IMMUNE_RESPONSE and immune infiltrated cells. (E) The correlation heatmap of 20 immune genes with 12 cuproptosis genes. (* P<0.05, **P< 0.01).

Furthermore, we identified 20 prognostic DEGs which overlapped with the GO_BP_IMMUNE_RESPONSE gene set and correlated with immune cell infiltration in TNBC by the TIMER algorithm (Figure 3D). We also validated the significant correlation between these 20 prognostic DEGs and the 12 cuproptosis-related genes (Figure 3E). Therefore, we considered these 20 genes as cuproptosis and immune-related genes (CIRGs) which provided an opportunity for us to recognize the interaction between cuproptosis and tumor immunity.

Construction of prognostic cuproptosis and immune-related gene signature

The prognostic value of CIRGs were visualized by the forest plot (Figure 4A). Then, LASSO-Cox method successfully selected nine hub genes to construct a cuproptosis and immune related signature (CIS) and divided TNBC patients into high- and low- CIS groups according to the score of each TNBC patients (Figure 4B). The coefficient of each hub gene was displayed by the Lollipop plot (Figure 4C). The different mRNA expression profile of nine hub genes in high- and low-CIS groups were displayed by a heatmap (Figure 4D). Importantly, the overall survival of TNBC patients in high-CIS group was significantly better than low-CIS group in the TCGA cohort (Figure 5A). In addition, the area under curve (AUC) values of 1, 3, 5-year ROC curved were larger than 0.7, suggesting a satisfied predictive accuracy of CIS score for TNBC prognosis (Figure 5B). The distribution of patients with different survival status and CIS score also suggested that a high CIS score was associated with a worse outcome (Figure 5C). Consistently, the predictive ability of CIS was validated in GSE81540 and GSE58812 cohorts (Supplementary Figures 1A, 1B, 2A, 2B). The overall survival analysis showed poorer prognosis in high-CIS group, compared with low-CIS group, and the ROC curves indicated the high predictive accuracy of CIS for predicting TNBC overall survival (AUC value of 1-, 3-, 5-year OS was 0.705, 0.815 and 0.786 respectively) in GSE81540 cohort (Supplementary Figure 1A, 1B) and (AUC value of 1-, 3-, 5-year OS was 0.869, 0.752 and 0.763 respectively) in GSE58812 cohort (Supplementary Figure 2A, 2B). The distribution of CIS score and survival status of patients showed that higher CIS score was related to poorer prognosis in the validation set (Supplementary Figures 1C, 1D, 2C, 2D). In conclusion, these results revealed that CIS had a robust performance to predict prognosis of TNBC.

Figure 4. Construction of cuproptosis-related immune signature (CIS) based on the LASSO-Cox model. (A) Forest plot of hazard ratios manifesting the prognostic values of cuproptosis-related immune genes. (B) LASSO-Cox regression analysis for variable selection and construction of signature for TNBC. (C) The corresponding coefficient of nine selected genes. (D) The expression profile of nine selected genes in high- and low-risk group.

Figure 5. The clinical significance of cuproptosis-related immune signature in TNBC. (A) The Kaplan-Merier curves for overall survival of low- and high-risk groups. (B) The ROC curves of CIS for predicting 1-, 3-, and 5-year overall survival in the TGCA-TNBC cohort. (C) The distribution of TNBC patients according to CIS score and survival status. (D, E) Univariate Cox analysis (D) and, multivariate Cox analysis (E) of clinical factors and CIS score in TCGA cohort. (F) Construction of nomogram for predicting 1-, 3-, and 5-year overall survival of TNBC patients. (G) Calibration curves of the nomogram for predicting the probability of OS at 1-, 3-, and 5-years survival in the TCGA cohort.

CIS is an independent prognostic factor for TNBC

To identify the clinical factors associated with the prognosis of TNBC, we used univariate and multivariate Cox regression analysis to estimate the hazard ratio of CIS score and common clinicopathological features in TCGA-TNBC cohorts. The results of univariate Cox regression analysis showed that CIS was a strong risk factor for OS in TNBC patients (In TCGA cohort, Hazard Ratio (HR): 3.57, 95% Confidence Interval (CI): 2.26-5.66, P < 0.001, Figure 5D; in GSE81540 cohort, HR: 1.82, 95% CI:1.29-2.77, P = 0.004; Supplementary Figure 1E). The results of multivariate Cox regression analysis demonstrated that CIS score was an independent prognostic factor for TNBC patients after adjusting with clinical features (In TCGA, HR: 3.02, 95% CI: 1.81-5.04, P < 0.001 (Figure 5E); in GSE81540 cohort, HR: 1.02, 95% CI:1.32-1.87, P = 0.001; Supplementary Figure 1F). Thus, the above results indicated that CIS was an independent and effective prognostic factor for overall survival of TNBC patients. A nomogram was constructed based on CIS and prognostic clinical characteristics to predict 1-, 3-, and 5-year OS of individual TNBC patient to promote the application of CIS in clinical practice (Figure 5F). As expected, the high-nomogram total point was associated with a worse outcome in the TCGA-TNBC cohort (C-index = 0.829, P < 0.001). The calibration curves were highly concordance with the actual rates for 1-, 3-, and 5-year survival, indicating an excellent predictive capacity of the nomogram (Figure 5G).

Association between CIS and tumor immune phenotype

To gain insights into the immune status of individual TNBC patients, we firstly assessed immune cell infiltration and the expression of immune checkpoint genes in TNBC samples. The TIMER analysis revealed a significantly higher infiltration fraction of B cells, T cells CD8+ and T cells CD4+ in the high-CIS group (Figure 6A). Additionally, the CIS score exhibited a positive correlation with immune cell infiltration across various analyses (Figure 6B). Consistent with previous results, the CIBERSORT analysis showed that a high CIS score was correlated with high infiltration of B cells, T cells CD8+, T cells CD4+ and low infiltration of macrophages. Similar results were observed using MCPcounter, xCell, EPIC, and quantiseq algorithms (Figure 6B).

Figure 6. Association between CIS and tumor immune microenvironment. (A) The boxplot of the immune cell infiltration in low- and high-risk groups depending on TIMER estimation. (B) The significant correlation between immune cells and CIS score in CIBERSORT, MCPcounter, xCell, EPIC, and quantiseq analyses. (C) Comparison of the expression of genes from the signatures of Immune_check_point, HLA, B_cells_Rooney_et_al., B_cells_Danaher_et_al., T_cell_Danaher_et_al., and TIP_Cancer_antigen_presentation between low- and high -risk groups in heatmap and boxplot. (D) The correlation between CIS score and ESTIMATE algorithm evaluated scores in TNBC. (*P<0.05, **P< 0.01, ***P< 0.001).

In addition, we performed ESTIMATE analysis to further explore the immune status of individual TNBC patients. We found that CIS score was positively correlated to estimate_score and immune_score with statistical significance, and high CIS score was associated with high stromal score and low tumor purity in TNBC samples (Figure 6C). These findings suggested that CIS was significantly associated with the immune status of TNBC.

Moreover, we compared the expression of immune checkpoint genes between high- and low-CIS group. CTLA4, PDCD1, TIGIT were significantly higher expressed in high-CIS group. In the meanwhile, although no statistical significance, CD274 was slightly upregulation in high-CIS group (Figure 6D). Generally, dysregulated expression of HLA signaling genes has been reported to play important roles in the immune recognition, immune response and immunotherapy efficacy [20]. In this study, a majority of genes in HLA signature, including CR1/2, CD19, CD22, CCR5, MS4A1, were highly expressed in high-CIS group than that in low-CIS group, suggesting that the patients in high-CIS group might present induced immunological surveillance and boosted immune response that contribute to better response to immunotherapy (Figure 6C).

Expression level of genes in immune cell signatures were compared between high- and low-CIS groups to analysis the correlation between CIS and immune cell related genes (Figure 6C). Most of the genes in B_cells_Rooney_et_al. signature were upregulated in high-CIS group, including BTLA, FCRL3, BANK1, CD79A, BLK, FCRL1, BACH2. The expression of genes related to T_cells_Danaher_et_al. signature were investigated as well in TNCB samples, several genes such as CD3D/E/G and TRAT1 showed higher expression in high-CIS group than low-CIS group. Moreover, genes in TIP_Cancer_antigen_presentation signature were highly expressed in high-CIS group.

Above findings suggesting that CIS was significantly associated with anti-tumor immunity in TNBC, the patients in high-CIS group presented a higher rate of immune cell infiltration, higher expression of immune checkpoint genes and immune cell signature related genes, which might be more likely to benefit from immunotherapy, while the patients with low CIS score showed an immunosuppressive phenotype.

Predicition of response to immunotherapy and target drug

We further investigated whether the CIS score associated with the response to immunotherapy and predicted potential therapeutic drugs for patients in the low-CIS group. The TIDE algorithm was employed to predict the response of each TNBC patient to ICB treatments (Figure 7A). As anticipated, patients in the high-CIS group were more likely to exhibit favorable responses to immunotherapy (Fisher’s test, P < 0.05). The submap analysis further predicted that high CIS score might be significantly associated with anti-PD-1 response, while a low CIS score was related to resistance to anti-CTLA4 treatment (Figure 7B).

Figure 7. Predicting response to immunotherapy and potential therapeutic agents. (A) The distribution of the risk score and predicted immune therapy response based on the TIDE dataset (left). The numbers of patients who might potentially response or non-response to immunotherapy in two risk groups for Fisher’s test (right). (B) Submap analysis predicted potential response to immunotherapy in low- and high-risk groups. (C) Correlation between estimated AUC values of CTRP-derived compounds and CIS score. (D) Correlation between estimated AUC values of PRISM-derived compounds and CIS score. (E) Comparison of AUC values of CTRP-derived compounds between low- and high-risk groups. (F) Comparison of AUC values of PRISM-derived compounds between low- and high-risk groups. (*P<0.05, **P< 0.01, ***P< 0.001).

Since the low-CIS group might be insensitive to immunotherapies due to its immunosuppressive phenotype, we predicted the potential therapeutic small molecular agents for these patients. Four CTRP-derived compounds (KX2-391, gemcitabine, BI-2536, AZD8055) and four PRISM-derived compounds (SB743921, vandetanib, paclitaxel, BI-2536) were found to be sensitive in the low-CIS group (Figure 7C, 7D). These compounds had a negative correlation with the CIS score and lower estimated AUC values in the high-CIS group (Figure 7E, 7F).

Functional annotation of different CIS groups

To identify the biological processes associated with CIS in TNBC, GO and KEGG enrichment analysis were used for functional annotation for patients in different CIS groups. In GO analysis, high enrichment of the immune system process (GO: 0002376) was observed in the high-CIS group, while the cellular process (GO: 0009987) and metabolic process (GO: 0008152) were higher enriched in the low-CIS group (Figure 8A, 8B). The results of GSEA showed that KEGG pathways such as Steroid Hormone biosynthesis and Drug Metabolism Cytochrome P450 were highly enriched in the low-CIS group (Figure 8C). The immune-related pathways including B Cell Receptor Signaling Pathway and T cell Receptor Signaling Pathway were highly enriched in the high-CIS group (Figure 8D). Furthermore, since tumor metabolism was closely associated with tumor progression, we explored the difference of tumor metabolism regulation between two CIS groups by ssGSEA method to score 102 metabolism signatures based on gene expression profiles of TNBC samples [21]. We found that the scores of Glycogen Degradation, Drug Metabolism by Cytochrome P450, Steroid Biosynthesis and Tyrosine Metabolism were significantly higher in low-CIS group than high-CIS group (Figure 8E). Combining the results together, the dysregulated genes in high-CIS group mainly involved in regulation of tumor immunity, while genes in low-CIS group might participate in tumor metabolism.

Figure 8. Functional analyses for patients in two risk groups. (A) GO analysis based on differentially expressed genes in low-risk samples. (B) GO analysis based on differentially expressed genes in high-risk samples. (C) Identifying KEGG pathways significantly enriched in the low-risk samples by GSEA. (D) Identifying KEGG pathways significantly enriched in the high-risk samples by GSEA. (E) Comparing the scores of signatures relating to tumor metabolism between low- and high-risk groups. (*P<0.05, **P< 0.01, ***P< 0.001).

Validating the expression of gene expression

Previous studies have demonstrated that IDO1 and THY1 could serve as oncogenes in cancers. Similarly, our result predicted that IDO1 and THY1 were associated with high prognostic risk in TNBC patients. The mRNA expression retrieved from breast cancer and normal tissue in the TCGA cohort showed that IDO1 and THY1 were significantly higher in TNBC samples than in normal samples (Figure 9A, 9C). Then, the results of qRT-PCR assays detecting mRNA expression of IDO1 and THY1 in TNBC cell lines (MDA-MB-231, MDA-MB-468) and normal cell line (MCF10A) were consistent with TCGA samples (Figure 9B, 9D). Moreover, we detected the protein expression of IDO1 and THY1 in TNBC and adjacent normal tissues by IHC staining. These results confirmed that IDO1 and THY1 were upregulated in TNBC than normal cells and tissues (Figure 9E, 9F). In our analysis of TNBC scRNA-seq data, we observed elevated expression of IDO1 and THY1 in cluster 3 and cluster 2, respectively, compared to the other clusters (see Supplementary Figure 2E, 2F).

Figure 9. Upregulated expression of IDO1 and THY1 in TNBC cell lines and tissue samples. (A) Significantly different expression of the IDO1 among breast cancer tissues and normal tissues. (B) qRT-PCR examined the mRNA expression of IDO1 gene in MCF-10A, MDA-MB-231, and MDA-MB-468 cell lines. (C) Significantly different expression of the THY1 gene between breast cancer tissues and normal tissues. (D) qRT-PCR examined the mRNA expression of THY1 gene in MCF-10A, MDA-MB-231, and MDA-MB-468 cell lines. (E, F) IHC staining of IDO1 and THY1 protein in adjacent normal tissues and TNBC tissues. (*P<0.05, **P< 0.01, ***P< 0.001).

Data acquisition and processing

The clinicopathological information and mRNA expression data of TNBC were downloaded from The Cancer Genome Atlas TCGA (https://portal.gdc.cancer.gov/) and set as the training cohort. The data of GSE81540, GSE58812 and GSE188600 GSE58812 and GSE188600 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) using the “getGEO” function in the “GEOquery” R package. Log2 transformations were performed for all gene expression data. Somatic mutation and copy number variation (CNV) data were obtained from the TCGA database. A total of 156 patients from the TCGA-TNBC cohort were included in the training set, and 360 TNBC patients from the GSE81540 cohort were enrolled in the validation set. Cuproptosis-related genes were obtained from previous studies [21].

Consensus clustering

Cluster analysis was performed using the “ConsensusClusterPlus” R package, using a gglomerative pam clustering with a 1-pearson correlation distance and resampling 80% of the samples for 10 repetitions. The optimal number of clusters was determined using the empirical cumulative distribution function plot. In this study, the TNBC patients were divided into 2 clusters based on the gene expression level of 12 cuproptosis-related genes. In detail, we used the consensus clustering algorithm with 1,000 iterations by sampling 80% of the data in each iteration. Two cuproptosis clusters (namely, “cluster A” and “cluster B” groups) were distinguished to present different cuproptosis statuses.

The TNBC cells divided into 7 clusters through UMAP analysis by using “Seurat” R package.

Identification of cuproptosis and immune-related genes

The differentially expressed genes (DEGs) between cuproptosis cluster A and cluster B were screened out with limma R package with the criteria of p < 0.05, |Log2(Fold Change)| > 1.3. Then, the genes associated with the overall survival of TNBC patients in the TCGA cohort were identified by univariate Cox regression analysis with the criteria of P-value <0.05. The intersection of the DEGs and TNBC prognosis-related genes included 207 prognostic cuproptosis-related genes. Subsequently, enrichment analysis based on GO and HALLMARK pathways was performed to investigate the potential biological functions of the prognostic cuproptosis-related genes. Then the genes involved in the GO_BP_IMMUNE_RESPONSE pathway were selected as cuproptosis and immune-related genes (CIRGs).

Construction and validation of CIS scoring model

Firstly, we perform univariate Cox regression analysis to assess the prognostic value of 20 cuproptosis and immune-related genes (CIRGs) and visualized it by a forest plot. Then, LASSO-Cox regression was used to screen 9 hub genes and obtain corresponding coefficients for constructing the cuproptosis-associated immune signature CIS signature by using the R package “glmnet”. The calculation formula for the CIS risk score is shown below:

$$CIS={\displaystyle \sum _{i=1}^{n}coe{f}_i}\ast x_i$$

Where x_i is the mRNA expression level of each gene and Coefi is the corresponding coefficient. The CIS score was calculated for all patients according to the formula: Riskscore = (-0.3088)*CLEC4E + (-0.1055)*IDO1 + (-0.194)*GBP1 + (-0.1523)*APOBEC3C + (0.2277)*TAP1 + (-0.1583)*C1R + (-0.1891)*RAB2B + (0.6039)*THY1 + (-0.7136)*GPLD. We assessed the reliability of the risk score model by predicting the OS of TNBC through survival curves of both the train set and validation set. The cut-off value of the CIS score was selected as the median risk score of all the patients, and the patients were divided into two risk groups (low-CIS and high-CIS) according to the cut-off value.

Functional annotation and immune cell infiltration analysis

The Gene Ontology (GO), REACTOME, and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment in each TNBC sample were analyzed with the “ClusterProfiler” R package. Gene set enrichment analysis (GSEA) was applied to explore significantly correlated pathways in this study.

The CIBERSORT, MCPcounter, xCell, EPIC, and quantiseq algorithms were used to assess the immune cell infiltration of high- and low-CIS groups. The ESTIMATE (The Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression) method was used to calculate the estimate score, immune score, stromal score, and tumor purity of each TNBC sample. The genes from immune-related signatures summarized in “IOBR” R package published in previous studies were included in our study to evaluated the immune phenotype of TNBC patient. The signatures related to tumor metabolism were also measured by IOBR package via ssGSEA method [21].

Prediction of potential therapeutic agents and immunotherapy response

To predict the potential therapeutic drugs for TNBC patients in high- and low-CIS groups, two datasets (Cancer Therapeutics Response Portal (CTRP) 2.0 database (http://portals.broadinstitute.org/ctrp/) and the Profiling Relative Inhibition Simultaneously in Mixtures (PRISM) database (https://www.theprismlab.org/)) were used to obtain sensitive data. The AUC was used in the two datasets as a standard value for the evaluation of drug sensitivity, with a lower AUC value indicating better drug sensitivity. The drug sensitivity of each sample was predicted by the “pRophetic” package by calculating the AUC value of each compound in every TNBC sample. The correlation between the CIS score and the AUC of each drug was calculated, potentially sensitive agents were screened based on the relationship between AUC and CIS with the criterion of Spearman r > 0.3, and components with significantly different AUC values between the high- and low- CIS groups. Potential responses to immunotherapies (anti-PD1 and anti-CTLA4) of high- and low-risk groups were predicted by the Tumor Immune Dysfunction and Exclusion (TIDE) online Tool (http://tide.dfci.harvard.edu/) and Subclass Mapping method (https://cloud.genepattern.org/gp).

qRT-PCR assays

Total RNA was extracted from TNBC cells by RNAiso Plus (TaKaRa, Kyoto, Japan). The PrimeScript® RT Master Mix Perfect Real Time kit Q12 (TAKARA) was used to reverse transcription. PCR amplifications were performed with the SYBR PremixEx Taq II (TAKARA) and Bio-Rad system according to manufacturers’ protocols. All the primer pairs used in qRT-PCR were listed in Supplementary Table 3.

Immunohistochemistry staining

We collected 20 pairs of paraffin-embedded TNBC tissues and adjacent normal tissues from Peking University Cancer Hospital, and this study has been approved by the Ethics committees of Peking University Cancer Hospital. The TNBC tissues and adjacent normal breast tissues were fixed with 10% formalin, embedded by paraffin, and sectioned. After heating in tris-EDTA buffer, we blocked slides using 5% goat serum and incubated slides with primary antibody (anti-IDO1 [1:100 dilution, Cat No.: 13268-1-AP, Proteintech, China], anti-THY1/CD90 [1: 500 dilution, Cat No.: 66766-1-Ig, Proteintech, China]) at 4° C overnight. Then, the slides were incubated with a secondary antibody. The sections were scanned to obtain high-resolution digital images using a 3DHISTECH scanner (Pannoramic, Taipei). Protein expression levels were evaluated semi-quantitatively following the Allred scoring system guidelines and scored separately by two qualified pathologists [46].

Statistical analysis

In the present experiment, all statistical analysis was conducted by R 4.3. The Wilcoxon test and the Kruska-Wallis test were used for comparisons between two independent samples and comparisons among multiple samples for nonparametric data, respectively. The t-test was used for parametric data. P-value < 0.05 was considered statistically significant. Related R packages including “ggplot2”, “ggpubr”, “survival”, “survminer” and other related R packages were downloaded from Bioconductor packages or R packages.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.