The workflow of our analysis was shown in Figure 1. The gene expression RNAseq (HTSeq-FPKM), clinicopathological data, survival data, mutation data, breast cancer were downloaded from the online database TCGA(https://cancergenome.nih.gov). The data of TCGA breast cancer patients were randomly divided into training set and internal verification set at a ratio of 1:1. Chi square test verified the fairness of grouping.

According to gencode (https://www.gencodegenes.org), lncRNAs was extracted. Seven cuproptosis-related genes including FDX1, LIPT1, LIAS, DLD, DLAT, PDHA1 and PDHB (Tsvetkov et al., 2022). Pearson method was used to analyze the correlation between lncRNAs and cuproptosis-related genes in breast cancer. 153 LncRNAs with correlation coefficient R2 > .3 and p < .05 were considered as cuproptosis-related lncRNAs. The expression data of these 153 lncRNAs were extracted from the breast cancer data of TCGA. Then the data were equally divided into training set (n = 555) and internal verification set (n = 554). All RNA sequencing data were normalized by log2 conversion.

Univariate Cox regression analysis was used to obtain lncRNAs associated with prognosis from the training dataset. LncRNAs with p < .05 were then included in multivariate Cox risk model analysis to generate cuproptosis-related lncRNA model. To avoid the model being too complicated, Akaike information criterion (AIC) was used to evaluate the complexity of the model, and model with lowest AIC value (AIC = 695.15) and no reduction in prediction efficiency was evaluated as best. The risk score of each breast cancer patient was calculated according to the following formula: BCCuS (risk score) = γ1 × expression (lncRNA1)+ γ2 × expression (lncRNA2) + … + γn × expression (lncRNAn). γn is the regression coefficient of the corresponding lncRNA, expression (lncRNAn) is the expression level of lncRNA, and the unit is FPKM. According to the median of risk score, 555 breast cancer patients in the training data set were divided into two groups: high-risk group and low-risk group. To compare OS differences between high-risk and low-risk groups, Kaplan Meier plotter analysis was implemented. In addition, the predictive ability of the model was evaluated by performing ROC curves using the survivalROC R package. The diagnostic value of BCCuS was further validated by Risk Plot analysis by using pheatmap package. Finally, the expression and prognostic data of 554 patients in the test set were used to test the efficacy of the model.

Nomogram and C-index was used to prove that BCCuS has strong correlation with clinical outcomes and can tutor decision-making. Then patients were splitting into Stage I-II and Stage III-IV; Age < 60 and Age > 60, BCCuS still maintain inspection efficiency. TMB (tumor mutation burden) data and base mutation data of breast cancer samples were downloaded from TCGA database. Breast cancer samples were divided into high and low groups according to the median risk scores, and the frequency of gene mutations and TMB in the two groups was calculated by using maftools package. GSVA and limma package was used to analyze the differences in immune function between BCCuS-high and BCCuS-low risk groups. We also evaluated the relationship between BCCuS and tumor immune dysfunction and exclusion (TIDE) score (https://tide.dfci.harvard.edu/).

By using pRRophetic package and machine learning, 36 anticancer drugs whose sensitivity (IC50) was significantly related to the model were screened out. The threshold of p-value was .001, and the anticancer drug sensitivity database used was CPG 2016, which was included in pRRophetic package.

We used the primary breast cancer tissues and normal tissues adjacent to the cancer from 10 patients with breast cancer diagnosed in the breast center of Jiangsu Province, China. The deadline for follow-up was April 2022. All patients provided written informed consent. This was done in accordance with the declaration of Helsinki. All samples were obtained with the approval of the hospital ethics committee.

Total RNA was isolated from tissues and cells using Trizol reagent (Invitrogen, United States) according to the manufacturer’s protocol. CDNA was synthesized using hiscript II (vazyme, China). Then, qRT-PCR of mRNA and lncRNA was performed on a stepone plus real-time PCR system (Applied Biosystems, United States). U6 and β- Actin was used as a standard control for lncRNA and mRNA detection, respectively. The gene expression in PCR was obtained by logarithmic conversion of CT value. All PCR primers (lifetech, China) are listed in Table 1.

All statistical analysis and visualization were performed using R software (version 4.1.3). The R packages “ggplot2”, “ggalluvial” and “limma” were used to analyze cuproptosis related lncRNAs in TCGA breast cancer patients. Use Kaplan Meier plotter to evaluate the OS difference between the BCCuS-high and BCCuS-low risk groups. Univariate and multivariate Cox regression were constructed to assess whether the prognostic characteristics were independent of Age, Sex, Clinical Stage, and TNM stage. The above analysis is implemented using the R package “Survival”. The lasso (Least absolute shrinkage and selection operator) analysis was operated using R package “caret” and “glmnet”. The Nomogram is implemented using the “rms” R package to predict the OS. All p-value<.05 were considered statistically significant. The qPCR data were analyzed and mapped with GraphPad Prism8 software.

We obtained the expression of 1,109 breast cancer samples and 103 adjacent samples from TCGA dataset. A total of 13,493 lncRNAs were screened out. The lncRNAs, which differentially expressed in breast cancer and adjacent samples obtained by Wilcox test were included in the follow-up analysis (p < .05). We determined the cuproptosis related lncRNAs by correlation test. Pearson correlation analysis showed that 153 lncRNAs were associated with cuproptosis related genes in breast cancer (Cor <.3, p-value <.05). Sankey map was used to show the corresponding relationship between Cuproptosis related seven genes and 153 Cuproptosis related lncRNAs (Figure 2A).

The data of TCGA breast cancer patients were randomly divided into training set and internal verification set at a ratio of 1:1. Chi square test verified the fairness of grouping. The p values of chi-square test for all indicators between Test and Train were >.05, indicating the fairness of grouping.

By using univariate cox regression analysis (p < .05), we identified five lncRNAs highly correlated with patients’ OS from 153 cuproptosis related lncRNAs in the training data set (Figure 2B). Then Lasso analysis screened out two genes (USP2-AS1, NIFK-AS1) with the largest relative weight for subsequent model construction (Figure 2C). Partial Likelihood deviance of Lasso analysis was also provided (Figure 2D). Next, we used the multivariate COX regression to fit cuproptosis related lncRNAs risk model (BCCuS) of breast cancer patients. The formula of prognostic risk score: BCCuS = (1.117,558 × USP2-AS1)+(−1.05486 × NIFK-AS1) (Table 2). According to the median BCCuS score, 1,109 breast cancer patients in the training data set were divided into two groups: high- BCCuS group and low- BCCuS group. In Figures 2B, C and Table 2, we can know that lncRNA USP2-AS1 is a favorable prognostic factor and lncRNA NIFK-AS1 is an unfavorable prognostic factor for breast cancer patients.

According to the expression data of TCGA samples, the correlation heatmap shows the expression correlation between two lncRNA used to build the model and seven cuproptosis genes (Figure 3A). USP2-AS1 was positively correlated with four genes (DLAT, PDHA1, FDX1, DLD) and negatively correlated with two genes (LIAS, PDHB). NIFK-AS1 was positively correlated with three genes (LIAS, LIPT1, PDHB) and negatively correlated with two genes (DLD, DLAT) (Figure 3A). To further evaluate the predicting ability of BCCuS in breast cancer, Kaplan Meier survival analysis was conducted in the training cohort and internal validation set from TCGA data set, based on the patient OS (overall survival) data obtained from TCGA website and the patient PFS (prognosis free survival) data obtained from XENA website (https://xena.ucsc.edu/). In train set, the results showed that the OS of breast cancer patients in the high-risk group was worse than that in the low-risk group (p = .026) (Figure 3B), and the PFS was also worse than that in the low-risk group (p < .001) (Figure 3E). In internal validation set, the OS of breast cancer patients in the high-risk group was worse than that in the low-risk group (p = .033) (Figure 3C), as well as PFS (p = .012) (Figure 3F). In all TCGA-BRCA 1109 samples, the OS of breast cancer patients in the high-risk group was worse than that in the low-risk group (p = .001) (Figure 3D), as well as PFS (p = .024) (Figure 3G). Univariate and multivariate Cox regression were constructed to assess whether the prognostic characteristics were independent of age, sex, clinical Stage, and TNM stage (Figures 3H,I). Results show that BCCuS might be an independent prognostic factor and could partially eliminate the interference of clinical factors. Risk plot was mapped using R package “heatmap” and in both sets, USP2-AS1 was appeared to express higher in BCCuS-high group, NIFK-AS1 express higher in BCCuS-low group, which shows the reliability of the model gene (Figures 4A–C). Further validation of risk model was conducted using c to check the predicting efficiency, with AUC (Area Under Curve): 1-year: .730, 3-year: .677, 5-year: .610 (Figure 5A). ROC curve contains BCCuS, and other clinical factors shows BCCuS was a better predicting signature than clinical stage (BCCuS AUC = .730, stage AUC = .717) (Figure 5B). Concordance index of BCCuS and other clinical information indicated that BCCuS could be a reliable clinical reference index (Figure 5C). The nomogram was constructed based on the score for clinical convenient predicting of single-patient prognosis (Figure 5D). Further, KM-plot was mapped after grouping of overall samples according to clinical stages (stage 1-2: p = .024; stage 3-4: p = .003), and according to age (≤60: p = .005; >60: p = .011) (Figures 5E–H).

First, we used limma package to analyze the differentially expressed genes between BCCuS-high and BCCuS-low group (p < .001, |log2Foldchange|>1) (Figure 6A). 26 differential genes are screened out and displayed using volcano plot. Among these 26 genes, 18 genes down-regulated, and eight genes up-regulated, including calcium binding proteins S100A8 and S100A9. S100A8/A9 expresses in both the tumor microenvironment and the external environment and can be classified as a signal during the tumorigenesis (Moon et al., 2008; Woo et al., 2021). GO enrichment analysis reveals that differentially expressed genes enriched in Metabolic metabolism related molecular function and immune response (Ashburner et al., 2000), such as sequestering of metal ion and leukocyte migration (Figure 6B). KEGG results shows that BCCuS-high and BCCuS-low group enriched differentially in pathways related to immune response and energy metabolism, for example, serotonergic synapse and viral protein interaction with cytokine receptor (Figure 6C) (Kanehisa et al., 2017).

We use R package “GSVA” to calculate the significant difference enrichment of immunological function between high and low risk groups. Only type II IFN response was significantly inhibited in BCCuS-high group. HLA, T cell co-inhibition, check point, APC co-stimulation, CCR, APC co-inhibition, Para inflammation, MHC class I and Type I IFN response all significantly promoted in BCCuS-high group (p < .05) (Figure 7A). Tumor immune dysfunction score (TIDE) is consistent with tumor immune escape characteristics and can predict the effect of immunosuppression therapy (Jiang et al., 2018). BCCuS-high group shows significantly lower TIDE score, which may reveal the immune escape of tumor cells and may lead to poor prognosis (Figure 7B).

1,109 Breast cancer samples were divided into high and low groups according to the median risk scores. Frequency of mutations and TMB in both groups were calculated by using “maftools” R package. In BCCuS-low group, the mutation frequency of PIK3CA, CDH1, MAP3K1 up-regulated, these genes have been reported as tumor suppressor in breast cancer (Figure 8A) (Jiang et al., 2014; Couch et al., 2017; Ge et al., 2017; Zacksenhaus et al., 2017; Xue et al., 2018; Wijshake et al., 2021). While TP53 and TTN, were upregulated in BCCuS-high group, functioning as a promoter in breast cancer (Figure 8B) (Zhang et al., 2018a; Mao et al., 2018; Badve and Gökmen-Polar, 2019; Zheng et al., 2021). TMB (tumor mutation burden) data and base mutation data of breast cancer samples were downloaded from TCGA database. KM plot was mapped between TMB-high and TMB-low groups and TMB-high shows lower survival probability (p = .023) (Figure 8C). Then according to BCCuS and TMB, TCGA samples were divided into four groups, and survival curve was fitted. BCCuS-high and TMB-high group shows the highest risk (p = .005) (Figure 8D). Violin plot reveals that high-risk group has close ties with higher TMB (Figure 8E).

By using pRRophetic package and machine learning (Geeleher et al., 2014), 36 anticancer drugs whose sensitivity (IC50) was significantly related to the model were screened out. The threshold of p-value was set to .001. It is worth noting that several kinds of tumor chemotherapy drugs and targeted drugs have higher sensitivity in high-risk group, for example, Cisplatin, Masitinib, Gefitinib, and Tivozanib (Figure 9). Among 36 anticancer drugs, 31 had higher sensitivity in high BCCuS group, while five had lower sensitivity in high BCCuS group.

qPCR of primary breast cancer tissues and normal tissues adjacent to the cancer from 10 patients with breast cancer were operated. Results illustrated that two model lncRNAs both differentially expressed in cancer and Para cancerous samples, USP2-AS1 and NIFK-AS1 both had lower expression in cancer samples. FDX1, PDHA1 and DLAT also significantly down regulated in breast cancer. Correlation heatmap of seven cuproptosis-related genes and two model lncRNAs was constructed using experimental expression. The result shows the experimental expression of USP2-AS1 was positively correlated with PDHA1, NIFK-AS1 was positively correlated with LIPT1, this is consistent with the previous bioinformatic prediction.