Identification and study of cuproptosis-related genes in prognostic model of multiple myeloma

ABSTRACT

Background

Multiple myeloma (MM) is a highly heterogeneous disease. Cuproptosis is a novel mode of death that is closely associated with several diseases, such as hepatocellular carcinoma. However, its role in MM is unknown.

Methods

MM transcriptomic and clinical data were obtained from UCSC Xena and gene expression omnibus (GEO) databases. Following MM samples were divided into different subtypes based on the cuproptosis genes, the differentially expressed genes (DEGs) among different subtypes, namely, candidate cuproptosis related genes were analyzed by univariate Cox and least absolute shrinkage and selection operator (LASSO) regression to construct a cuproptosis-related risk model. After the independent prognostic analysis was performed, a nomogram was constructed. Finally, Functional enrichment analysis and immune infiltration analysis were performed in the high- and low-risk groups, potential therapeutic agents were then predicted.

Results

The 784 MM samples in UCSC Xena cohorts were divided into three different subtypes, and 4 out of 346 candidate cuproptosis related genes, namely CDKN2A, BCL3, KCNA3 and TTC14 were used to construct a risk model. Risk score was considered a reliable independent prognostic factor for MM patients. It was investigated that the pathway of cell cycle was significantly enriched in the high-risk group. In addition, immune score, ESTIMATE score and cytolytic activity were significantly different between different risk groups, as well as 13 immune cells such as memory B cells. Nine drugs were predicted in our study.

Conclusion

A cuproptosis-related prognostic model was constructed, which may have a potential guiding role in the treatment of MM.

KEYWORDS:

• Multiple myeloma
• cuproptosis
• prognosis
• immune infiltration
• consensus clustering

Abbreviations

CRGS	=	cuproptosis-related genes
GEO	=	gene expression omnibus
LASSO	=	least absolute shrinkage and selection operator
MDSCs	=	myeloid derived suppressor cells
MGUS	=	monoclonal gamma globulin disease
MM	=	multiple myeloma
OS	=	overall survival
SCT	=	stem cell transplantation
Tregs	=	tregulatory cells

1. Introduction

Multiple myeloma (MM) is a hematologic malignancy characterized by abnormal accumulation of myeloma plasma cells and secretion of monoclonal immunoglobulins, and causes osteolytic lesions and extramedullary organ invasion in the advanced stages of the disease [1]. In the past 20 years, application of new drugs such as proteasome inhibitors (e.g. bortezomib and carfifizomib), immunomodulatory drugs (e.g. thalidomide and lenalidomide), hematopoietic stem cell transplantation (SCT) have significantly improved complete remission rate and overall survival of MM patients [2]. Nevertheless, MM is still an incurable malignancy [3]. Response to treatment and survival of MM is heterogeneous and median overall survival (OS) ranges from 2 to 10 years. This apparent difference may be related to the heterogeneity of myeloma cell biology and multiple host factors [4]. Therefore, it is critical to identify disease-related biomarkers and use them to distinguish patients with different prognosis, and which will facilitate the development of individualized treatment to cope with tumor heterogeneity and improve the prognosis of patients.

Copper (Cu) is an essential mineral nutrient for all living organisms and it is fundamental in a huge number of biological processes, including mitochondrial respiration, iron uptake, antioxidant/detoxification processes [5]. In the field of cancer, several studies have reported higher levels of Cu in a variety of malignancies compared to normal tissues. It is evident that the dyshomeostasis of Cu plays a prominent role in cancer [6]. It has been reported that Cu may have anti-myeloma activity in vivo and not facilitate myeloma progress [7]. In recent decades, a growing body of evidence suggests that the Cu-complex might serve as a potential therapeutic target for treating cancer, as Cu has been shown to promote the death of cancer cells via apoptosis and/or an accumulation of free radicals [7,8]. Numerous studies have shown that the combination of DSF and copper (DSF/Cu) greatly enhances its antitumor efficacy [9]. Xu et al. demonstrated that the DSF/Cu complex induced apoptosis of MM cells and MM primary cells [10]. The findings showed that DSF/Cu has potent anti-myeloma activity in vitro and in vivo highlighting valuable clinical potential of DSF/Cu in MM treatment. Anshul Agarwal et al. showed that a combination of resveratrol-copper reduced the transplant related toxicities in patients with MM receiving high dose melphalan [11].

Recently, a novel way of copper-induced cell death that distinct from apoptosis, necroptosis, pyroptosis, and ferroptosis has been discovered and termed ‘cuproptosis’ by a study published on Science. The cuproptosis is a regulated cell death with copper-dependent and mitochondria-dependent, it occurs through direct binding of copper to lipid acylation component of the TCA cycle, resulting in abnormal aggregation of lipid acylation proteins and loss of iron-sulfur tuftin, and also causing cell death due to proteotoxic stress [12]. A lot of studies have suggested that cuproptosis directly related to the progression of several diseases, including Huntington’s disease (Hd), amyotrophic lateral sclerosis (ALS), and Alzheimer’s disease (Ad) [13]. These findings regarding cuproptosis provide new insights into potential applications for treating cancer. However, cuproptosis as a newly discovered form of regulating cell death, its potential in being the target of immunotherapy is far to understand. Besides, the role of cuproptosis in the onset, development and prognosis is still unknown.

In this study, we explored the correlation between cuproptosis-related genes (CRGs) and prognosis of MM through bioinformatics analysis, providing a new thought for the treatment and prognosis evaluation of MM. Moreover, our analysis highlights the importance of CRGs in MM development and lays a foundation for the therapeutic application of cuproptosis regulators in MM.

2. Materials and methods

2.1. Data source

In this study, 860 MM transcriptomic data and clinical data were downloaded from the UCSC Xena database (http://xena.ucsc.edu/), of which 784 MM samples had survival status and survival time. MM transcriptome data and clinical datasets were downloaded from the gene expression omnibus (GEO) databases (https://www.ncbi.nlm.nih.gov/geo/). In addition, GSE47552 contains 5 controls and 41 MM bone marrow donors, its chip detection platform is GPL6244 [HuGene 1 0 st] Affymetrix HumanGene 1.0 ST Array [transcript(gene)version]. GSE16558 contains 5 controls and 60 MM bone marrow donors. GSE4581 contains 414 MM samples with survival time and status, its chip detection platform is GPL570 [HG U133 Plus 2] Affymetrix HumanGenome U133 Plus 2.0 Array. Moreover, 10 cuproptosis genes were obtained [14], namely FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS and CDKN2A.

2.2. Identification of candidate cuproptosis related genes

The 784 MM samples from the UCSC Xena database source were classified different subtypes according to K values using the ConsensusClusterPlus package (version 1.54.0) [15], and subtypes were analyzed respectively for survival and clinical differences using the survival package and chi-square test. In addition, spearman correlation analysis was performed for 10 cuproptosis genes. Finally, we performed differential analysis of genes between subclasses to obtain candidate cuproptosis related gene.

2.3. Screening of candidate genes

DEGs were obtained by differential analysis of the GSE47552 and GSE16558 datasets using the limma package with the filter conditions of P < 0.05, | log2foldchange | ≥ 0.5[16]. The candidate cuproptosis related genes were intersected with the DEGs to obtain the intersected genes. Then, the candidate genes were screened by correlation analysis using Pearson setting screening criteria of P < 0.05, |Cor| > 0.5 for the intersecting genes and the 10 cuproptosis genes. Moreover, heat map of cuproptosis related genes was drawn using the pheatmap package.

2.4. Functional enrichment analysis

The study was carried out by enrichment analysis of candidate gene based on gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGg) databases using the clusterProfiler package, and the screening conditions were P < 0.05 and counts ≥ 2. And the pathway was displayed by drawing a bar graph in the enrichplot package (version 1.10.2) [17].

2.5. Construction of a risk model

The study extracted the expression data of candidate gene in 784 MM samples and combined them with clinical data such as overall survival (OS) to obtain the clinical expression data of MM samples. Then, the survival related cuproptosis genes were obtained by univariate Cox regression analysis using the survival package (P < 0.05) (version 3.2-3). Least absolute shrinkage and selection operator (LASSO) logistic regression was implemented by setting the parameter family to Cox using the glmnet package (version 4.1-1) [18], and the model genes were selected to construct the risk model.

2.6. Evaluation and verification of risk model

Firstly, the risk value of each MM patient was derived from the expression of model genes and LASSO regression, and the patients were divided into high-risk and low-risk groups by the median risk score, the risk score was calculated as follows: risk score = β1×X1+β2×X2+ … +βn × Xn. Survminer package (version0.4.8) was used to draw survival curve to compare the survival status of patients in high and low risk groups [19], and then survival was used to draw Kaplan–Meier (K–M) and receiver operating characteristic (ROC) curves to test the effectiveness of the model. Finally, the model is externally validated using the GSE4581 dataset.

2.7. Correlation analysis between risk score and clinical data

Among the 784 MM samples from UCSC Xena database, the correlation between age, gender, race, stage, sample type, survival status and risk score were analyzed. Chi-square test was used to analyze the differences in the number of patients with different clinical subtypes between high and low risk groups, and rank sum test was used to compare the differences in risk scores under different clinical information.

2.8. Independent prognostic analysis and construction of a nomogram

The age, gender and other clinicopathological factors and risk scores of 784 MM samples from UCSC Xena database were analyzed by univariate and multivariate Cox independent prognostic analysis. We used rms package to build a nomogram to show the prediction effect of risk model and used psm to plot the correction curve to further verify the accuracy [20].

2.9. Functional enrichment analysis of high and low risk groups

This study was performed by using the clusterProfiler (version 3.18.1) and org.Hs.eg.db (version 3.12.0) packages for GSEA in high and low risk groups. Differential analysis of genes was performed between high and low risk groups to obtain logFC of genes and rank them, and then gene set enrichment analysis (GSEA) was performed using the ranking results. The GSEA ridge chart was drawn by enrichplot (version 1.10.2) to visualize the TOP10 of the results of the high and low risk groups (sorted by the default corrected P value).

2.10. Immune infiltration analysis

ESTIMATE analysis was performed on 784 MM samples from the UCSC Xena database. Cell type identification by estimating relative subsets of RNA transcripts (CIBERSORT) algorithm (version 1.03) and LM22 gene set were used to calculate the proportion of 22 immune cells in MM samples, and samples with P > 0.05 were excluded. Moreover, the expression values of GZMA and PRF1 were extracted, and the FPKM value was converted to TPM, and the geometric mean of the two was calculated to obtain the cytolytic activity (CYT) of each sample. Single-sample gene set enrichment analysis (ssGSEA) was used to calculate the antigen presentation mechanism (APM), tumor infiltrating lymphoid cells (TILs) and T cell infiltration relative score (TIS) of each sample. In addition, ggplot2 package and ggpubr package were used to draw box line chart and stacked bar chart to visualize the results [21], and rank sum test was used to compare the results between high and low risk groups.

2.11. Drug sensitivity analysis

In this study, 138 chemotherapeutic drugs commonly used in oncology treatment were obtained from the genomics of drug sensitivity in cancer (GDSc) database (https://www.cancerrxgene.org/). The 50% inhibitory concentration (IC50) of each drug in each sample was first calculated based on MM gene expression data using the pRRophetic (version0.5) package [22], and then a correlation analysis between risk score and drug IC50 was performed using spearman to screen the drugs that met the criteria and compare IC50 between high and low risk groups to predict the potential drugs for treating MM.

3. Results

3.1. Classification of MM subtypes and their survival and clinical differences

The 784 MM samples were divided into three clusters (cluster 1 = 385, cluster 2 = 171, cluster 3 = 228) by consistent clustering analysis (Figure 1A,B). Survival analysis yielded significant differences in survival between subtypes, and cluster 3 patients had a worse survival prognosis (P < 0.0001) (Figure 1C). Analysis of clinical differences between subtypes using chi-square tests yielded significant differences in survival status, stage and race type between subtypes (Supplementary Figure 1). The heat map demonstrated that the differential expression of cuproptosis genes in different subtypes (Figure 1D). As seen by spearman correlation analysis, there was mainly a positive correlation between cuproptosis genes in MM (Figure 1E,F), such as LIAS with LIPT1 and PDHA1 were positively correlated.

Figure 1. Exploration for three cuproptosis-related MM subtypes. (A) Cumulative distribution function (CDF) curve, CDF delta area curve as well as the delta area curve of consensus clustering of samples from UCSC Xena exhibiting the relative alterations in the area under the CDF curve for each category number k compared with category number k-−1. (B) Sample clustering heatmap when consensus k = 3. (C) Survival analyzes for the three clusters (P < 0.0001). (D) Heatmap of cuproptosis-related genes expression in MM samples with different clinical data. (E) Correlation bubble diagram among cuproptosis genes. (F) The correlation network among cuproptosis genes.

3.2. Screening and functional enrichment analysis of candidate genes

Because of the significant difference in survival between cluster 3 and cluster 1 patients and the high number of patients in cluster 1, we performed a differential analysis between cluster 3 and cluster 1. The volcano plot (Figure 2A) showed there were a total of 346 DEGs between Cluster 3 and Cluster 1, of which 101 genes were up-regulated in cluster 3 and 245 genes were down-regulated in cluster 3. The DEGs between cluster 3 and cluster 1 were candidate cuproptosis related genes.

Figure 2. Identification and evaluation of candidate cuproptosis related genes and candidate genes. (A) Volcano plot of 346 DEGs between Cluster 3 and Cluster 1. (B) Volcano plot of 2931 DEGs between MM and normal samples in GSE47552. (C) Volcano plot of 2817 DEGs between MM and normal samples in GSE16558. For a–c, the screening criteria are set to P < 0.05, | log₂foldchange | ≥ 0.5. (D) Venn diagrams for up-regulated and down-regulated intersecting genes by overlapping DEGs and candidate cuproptosis related genes. (E) Heatmap for candidate genes expression between MM and normal samples in GSE47552. (F) Heatmap for 33 candidate genes expression between MM and normal samples in GSE16558. (G)The Gene Ontology (GO) analysis for candidate genes, including the enriched biological processes, cellular composition terms. (H) The top 4 enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) terms of candidate genes.

Differential analysis of the datasets GSE47552 and GSE16558 showed that a total of 2931 DEGs existed between MM and normal samples in GSE47552, of which 866 genes were up-regulated in MM samples and 2065 genes were down-regulated in MM samples (Figure 2B). A total of 2817 DEGs were present between MM and normal samples in GSE16558, of which 646 genes were up-regulated in MM samples and 2171 genes were down-regulated in MM samples (Figure 2C). Moreover, a total of 113 intersecting genes were obtained by taking the intersection of DEGs and candidate cuproptosis related genes, of which 23 were up-regulated and 90 were down-regulated (Figure 2D). The 33 genes screened by Pearson correlation analysis between the candidate cuproptosis related genes and the 10 cuproptosis genes were candidate genes (Table 1). The results of heat map plotting showed that candidate genes were differentially expressed in MM and normal samples (Figure 2E,F).

Table 1. Gene expression patterns of candidate genes.

	logFoldChange	AveExpr	P value	adjust P value	Type
ZBTB43	1.544748459	7.536689217	0.00002367361	0.0006497842641	up
ARID5A	1.042615	6.944444283	0.00000228484	0.0001151878029	up
PER1	0.83253721	7.145908435	0.01034641789	0.0547042450264	up
BCL3	0.714748337	7.4120375	0.00000801325	0.0002973500987	up
CDKN2A	0.559170063	6.717567109	0.00189882113	0.0169075257696	up
JRKL	−0.683364234	5.789102087	0.01181529817	0.0599163198804	down
KANK1	−0.707842112	5.994155848	0.01174341636	0.0596875870062	down
TTC14	−0.80318762	6.402041783	0.00220477968	0.0189275782815	down
XRRA1	−0.811652639	6.083252674	0.00766709658	0.0451257433726	down
TRMT13	−0.991852215	4.687170609	0.00000662633	0.0002594152620	down
MEF2C	−1.171066361	9.380368261	0.00088668751	0.0098097251531	down
ZNF780B	−1.210632273	5.728114391	0.00048748021	0.0063666977262	down
MLEC	−1.215633239	7.829310609	0.00003665635	0.0009159849807	down
PEX3	−1.232453766	4.830419	0.00000365124	0.0001648492703	down
PREB	−1.277475029	7.793607152	0.00001099724	0.0003736037440	down
BTN3A1	−1.369683356	7.80929187	0.00003133974	0.0008137038935	down
INPP4A	−1.381296985	7.788920391	0.00002760596	0.0007355350612	down
EAF2	−1.417448088	7.844076957	0.00023588483	0.0036257071673	down
AARS1	−1.460543112	7.997751174	0.00071162614	0.0083963589838	down
HSD17B4	−1.494252302	6.583870826	0.00090880375	0.0099588304279	down
DDX17	−1.496989312	8.872141717	0.00000616417	0.0002467378924	down
BTK	−1.519907346	6.768366174	0.00007457311	0.0015499017683	down
SEC23B	−1.590185946	8.401603152	0.00000079091	0.0000526015975	down
ATG4C	−1.602991488	5.699794717	0.00007536576	0.0015607482059	down
SLC35F5	−1.642497834	6.936458739	0.00000536976	0.0002223055541	down
LAX1	−1.719442605	7.60914313	0.00001178321	0.0003930059803	down
IRF2	−1.813843478	7.233650022	0.00000000055	0.0000002555443	down
SETDB2	−1.830299073	5.286186478	0.00000000676	0.0000014251953	down
PREPL	−1.888035932	6.551051565	0.00000992715	0.0003467935671	down
MANEA	−2.016639415	8.218766522	0.00001862183	0.0005384759578	down
KCNA3	−2.426247259	8.04320887	0.00000000240	0.0000007131565	down
TRIP11	−2.456738639	5.41834337	0.00000000098	0.0000003947739	down
CCR2	−3.205685102	6.63999813	0.00001908149	0.0005445188728	down

GO functional enrichment results in terms of biological processes showed that candidate genes were associated with the regulation of NMDA receptor activity, leukocyte proliferation, myeloid apoptotic process, monocyte proliferation, lymphocyte proliferation, positive regulation of leukocyte apoptotic process, regulation of myeloid apoptotic process, interferon-gamma production, leukocyte apoptotic process, inflammatory apoptotic process, regulation of skeletal muscle cell differentiation, regulation of leukocyte apoptotic process, germinal center formation, macrophage apoptotic process, and regulation of macrophage apoptotic process were significantly correlated (Figure 2G). In terms of cellular composition, candidate genes were significantly associated with intrinsic components of peroxisomal membranes, COPII vesicle coat, components of peroxisomal membranes, microsomes, peroxisomes, membrane regions, microsomal membranes, peroxisomal membranes, nuclear spots, and endoplasmic reticulum exit sites (Figure 2G).

KEGG functional enrichment analysis showed that candidate genes were significantly associated with biosynthesis of unsaturated fatty acids, protein processing in endoplasmic reticulum, primary bile acid biosynthesis, and peroxisome pathways, and bar graphs were drawn to demonstrate the pathways (Figure 2H).

3.3. Construction and assessment of a prognostic risk model

The expression data of 33 candidate genes were combined with clinical data such as OS to obtain the clinical expression data of MM samples, and the results of univariate Cox analysis showed that a total of 6 survival related cuproptosis genes were obtained (P < 0.05) (Figure 3A), as shown by the survival curves (Figure 3B), the median values of all genes except SLC35F5 divided the samples into two subtypes with significant differences in survival. The gene coefficients plotted by LASSO regression analysis and the error plot for 10-fold cross-validation are shown (Figure 3C) that when the lambda value was taken to be the smallest, the constructed model was the optimal model. When the lambdamin was taken as 0.0081, four genes namely CDKN2A, BCL3, KCNA3 and TTC14 were screened (Table 2).

Figure 3. Identification of prognostic genes (A) Forest map for univariate Cox analysis to screen six survival related cuproptosis genes. (B) Kaplan–Meier (K–M) survival curve of 784 MM samples with different expressed patterns of six survival related cuproptosis genes, including CDKN2A (P < 0.0001), BCL3 (P = 0.018), KCNA3 (P = 0.0019), TTC14 (P = 0.0053), XRRA1 (P = 0.00093), and SLC35F5 (P < 0.38). (C) Least absolute and selection operator (LASSO) regression analysis to screen prognostic genes. LASSO coefficient profiles (top) and cross-validation (bottom) to select the optimal tuning parameter.

Table 2. Four model genes screened using LASSO COX regression analysis.

	coef	exp(coef)	se(coef)	z	p
CDKN2A	0.54	1.71	0.14	3.8	0.0001
BCL3	0.13	1.14	0.11	1.2	0.2
KCNA3	−0.15	0.86	0.1	−1.5	0.1
TTC14	−0.04	0.96	0.14	−0.3	0.7

The median risk score (0.6124) was calculated from the expression of the four model genes and the risk coefficient obtained by LASSO regression, and the MM patients were divided into high-risk and low-risk groups. The risk profile plots showed differences in risk scores, survival time and expression of model genes between the high- and low-risk groups (Figure 4A). Survival curve plots for the high- and low-risk groups showed significant differences in patient survival between the high- and low-risk groups (P < 0.0001), and patients in the high-risk group had a lower survival rate (Figure 4B). The ROC curves were plotted based on the risk model with a survival time node of 1–5 years (Figure 4C), and the results yielded that the AUC of the ROC curves were all greater than 0.6, indicating a good efficacy of the risk model. In addition, the risk model was externally validated using the dataset GSE4581, and the same risk score was calculated for each MM patient, and the MedianRiskscore = 9.4291 was calculated as the boundary to classify patients into high and low risk groups. The risk curve, scatter plot, heat map (Figure 4D), survival curve for high- and low-risk groups (Figure 4E) and ROC curve (Figure 4F) showed that the survival prognosis was worse in the high-risk group in the validation set data (P = 0.005), and the AUCs were all greater than 0.6, which was consistent with the results in the training set.

Figure 4. Construction and assessment of the prognostic risk model. (A) Risk score curve, survival status and gene expression of four prognostic genes between high- and low-risk groups in training set. (B) K–M survival curve of high- and low-risk groups in training set (P < 0.0001). (C) Receiver operating characteristic (ROC) curves of training set. (D) Risk score curve, survival status and gene expression of four prognostic genes between high- and low-risk groups in validation set (GSE4581). (E) K–M survival curve of high- and low-risk groups in validation set. (F) ROC curves of validation set (P < 0.005).

3.4. Correlation analysis of risk model and clinical factors and nomogram construction

The number of patients with different clinical subtypes was compared between the high and low risk groups using chi-square test, and the trilinear table of clinical information showed that there were significant differences between patients in the high and low risk groups between the subtypes of survival status, sampletype, Stage and cluster (Table 3). We used the rank sum test to compare the risk scores under different clinical information, it is evident from the results that there is a significant difference in the risk scores of stage, sample type, cluster and vital subtypes (P < 0.05) (Supplementary Figure 2).

Table 3. Clinical correlation analysis between the risk score and clinical factors.

	Total	Risk		P-value
		high	low
age(year)
Mean (SD)	63.5 (±10.5)	63.7 (±9.8)	63.4 (±11.2)	0.94
gender
female	323 (41.2%)	169 (43.1%)	154 (39.3%)	0.31
male	461 (58.8%)	223 (56.9%)	238 (60.7%)
Vital
Alive	622 (79.3%)	287 (73.2%)	335 (85.5%)	< 0.001
Dead	162 (20.7%)	105 (26.8%)	57 (14.5%)
Race
asian	12 (1.5%)	7 (1.8%)	5 (1.3%)	< 0.001
Black or African American	104 (13.3%)	43 (11.0%)	61 (15.6%)
not allowed to collect	16 (2.0%)	7 (1.8%)	9 (2.3%)
not reported	135 (17.2%)	110 (28.1%)	25 (6.4%)
other	4 (0.5%)	0 (0.0%)	4 (1.0%)
white	513 (65.4%)	225 (57.4%)	288 (73.5%)
Sample_type
Primary	727 (92.7%)	353 (90.1%)	374 (95.4%)	0.006
Recurrent	57 (7.3%)	39 (9.9%)	18 (4.6%)
Stage
I	266 (34.9%)	110 (28.4%)	156 (41.5%)	< 0.001
II	276 (36.2%)	158 (40.8%)	118 (31.4%)
III	221 (29.0%)	119 (30.7%)	102 (27.1%)
Cluster
cluster1	375 (47.8%)	134 (34.2%)	241 (61.5%)	< 0.001
cluster2	167 (21.3%)	24 (6.1%)	143 (36.5%)
cluster3	242 (30.9%)	234 (59.7%)	8 (2.0%)

The results of univariate and multivariate Cox independent prognostic analysis of the risk model showed P < 0.05 for riskScore (Figure 5A,B), it indicated that the risk model constructed is a more reliable independent prognostic factor for MM patients. In addition, the nomogram was plotted and the C-index was calculated to be 0.7295 (Figure 5C), it indicated that the model was a better predictor. The calibration curve was plotted in this study (Figure 5D), from which it was known that the risk model has a high accuracy in predicting the survival rate of patients at 1–3 years, it indicated that the constructed prediction model could be used as a valid model.

Figure 5. Independent prognostic analysis and construction of the nomogram. (A) Forest map of univariate Cox independent prognostic analysis. (B) Forest map of multivariate Cox independent prognostic analysis. (C) nomogram was constructed based on the risk model and other independent prognostic factors. (D) Calibration curve of nomogram for survival prediction at 1-, 2-, 3-years. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001.

3.5. Differences in enrichment function between high and low risk groups

In terms of KEGG pathways, cell cycle, complement and coagulation cascades, neutrophil outer trap formation, nitrogen metabolism, p53 signaling pathway, ribosome, and spliceosome pathways were significantly enriched in the high-risk group (Figure 6A). In contrast, ascorbate and aldehyde metabolism, cofactor biosynthesis, fatty acid degradation, fatty acid metabolism, glycosaminoglycan biosynthesis, lysosomes, interconversion of pentose and glucuronide, peroxisomes, propionate metabolism, valine, leucine and isoleucine degradation were significantly enriched in the low-risk group (Figure 6B).

Figure 6. Gene set enrichment analysis (GSEA) for the high- and low-risk groups. (A) The activated KEGG terms in the high-risk group. (B) The activated KEGG terms in the low-risk group. (C) The enriched GO terms in the high-risk group. (D) The enriched GO terms in the low-risk group.

In terms of GO biological processes, antimicrobial humoral responses, antimicrobial peptide-mediated antimicrobial humoral immune responses, fungal defense responses, defense responses against Gram-negative bacteria, mucosal intrinsic immune responses, killing of cells in other organisms involved in symbiotic interactions, killing of cells in other organisms, mucosal immune responses, and responses against fungi were at high risk in the group was significantly enriched (Figure 6C). In contrast, cellular glucuronidation, glomerular vascular development, glucuronide metabolic processes, acetyl heparan sulfate proteoglycan biosynthetic processes, renal vascular development, posterior renal glomerular development, olefin compound biosynthetic processes, boosterol biosynthetic processes, renal system vasculature development, and glucuronide metabolic processes were significantly enriched in the low-risk group (Figure 6D).

3.6. Impact of risk model on immune heterogeneity

In this study, the results of ESTIMATE analysis of MM samples could be concluded that there were significant differences in immune scores and ESTIMATE scores between high and low risk groups (P < 0.05) (Figure 7A). In addition, the proportion of 22 immune cells in the samples was calculated using the CIBERSORT algorithm and the LM22 gene set to obtain the remaining samples high = 392, low = 392. We plotted box line plots and stacked bar graphs based on the proportion of immune cells in the remaining samples (Figure 7B,C). And meanwhile, the differences in immune cells infiltration were compared in the high and low risk groups using the rank sum test, which was plotted as a box line plot for visualization. As can be seen from Figure 7D, using the CIBERSORT algorithm, 13 immune cells such as memory B cells, plasma cells, and CD8 T cells differed between high and low risk groups (P < 0.05).

Figure 7. Immune related analyzes. (A) Box plot of stromal score, immune scores and ESTIMATE scores in high- and low-risk groups. (B) Histogram for 22 immune cells proportions in each MM patient. (C) Boxplot of the proportion of immune cells in the remaining samples. (D) Boxplot of 22 immune cells proportions in high- and low-risk groups. (E) Violin plot of cytolytic activity (CYT) levels in high- and low-risk groups. (F) Violin plot of antigen presentation mechanism (APM), tumor infiltrating lymphoid cells (TILs) and T cell infiltration relative score (TIS) in high- and low-risk groups. ns, not significant, * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001.

In addition, this study also used ssGSEA to calculate APM, TILs, and TIS for each sample and compared CYT, APM, TILs, and TIS between high- and low-risk groups using the rank sum test. The plotted box plots showed that CYT was significantly higher in the high-risk group than in the low-risk group (Figure 7E,F).

3.7. Potential chemotherapy drug forecast

In this study, 9 drugs were obtained through spearman correlation analysis, indicating the IC50 of BI.D1870, SL.0101.1 and KU.55933 was negatively correlated with risk score, the correlation coefficient was less than-0.3. While the IC50 of AP.24534, Embelin, NSC.87877, JNK.9L, Axitinib and MK.2206 was positively correlated with risk score, and the correlation coefficient was more than 0.3 (Figure 8A). In addition, the IC50 values of 9 drugs were compared between high and low risk score groups. The results showed that there were significant differences among these drugs (P < 0.05). It showed that these drugs may be some potential chemotherapy drugs (Figure 8B).

Figure 8. Potential drug prediction (A) Spearman correlation analysis between risk score and drugs from drug sensitivity in cancer (GDSc) database. (B) The 50% inhibitory concentration (IC50) of nine drugs in high- and low-risk groups. **** P < 0.0001.

4. Discussion

MM is a malignant hematological disease frequently developed in the elderly [23,24]. The aging of the population has witnessed increasing incidence of MM in China [25]. MM usually develops from monoclonal gamma globulin disease (MGUS) of unknown significance [26]. Although application of new drugs has proved great success in the clinical treatment of MM, they fail to prevent the later development of resistance and the recurrence of the disease [27,28]. Therefore, MM is still considered a fatal disease with no cure. A comprehensive understanding of the genetic backgrounds and tumor microenvironment is ultimate for the prevention, treatment, and prognosis evaluation of MM.

As a kind of PCD, the role of cuproptosis in tumorigenesis remains unknown, so we evaluated the prognostic value of CRGs in MM. In this study, MM patients were divided into three clusters, and it was found that there were significant differences in survival status, stages and ethnic types among different clusters, and the expression of copper death-related genes in different clusters was also different, especially in cluster 1 and cluster 3. After that, DEGs was selected as a candidate copper-related gene based on cluster 1 and cluster 3. After machine learning, four CRGs (CDKN2A, BCL3, KCNA3 and TTC14) were screened out. We explored the expression signature of 4 CRGs (CDKN2A, BCL3, KCNA3 and TTC14) in MM tissues and examined their relationships with OS. Clinical correlation analysis showed that there were significant differences in risk scores for stage, sample type, cluster and vital subtypes. We also verified the risk score was a reliable independent prognostic factor for MM patients. The CDKN2A expression had a function in the cell cycle control and was strongly correlated with the origin of variety tumors (including glioma, pancreatic cancer, clear-cell renal cell carcinoma, etc.) [29–32]. Furthermore, CDKN2A was found to be silenced by methylation in MM [33]. Bcl-3 is an established oncogene in hematologic malignancies, such as B-cell chronic lymphocytic leukemias [34]. In cells, Bcl-3 can regulate cell differentiation, control cell proliferation, promote apoptosis, and participate in the cell cycle and DNA damage response. In the human body, Bcl-3 is closely related to various immune and inflammatory pathologies, and also related to occurrence, development, invasion and prognosis of various blood and solid tumors [35]. Furthermore, Brocke-Heidrich K et al. have demonstrated the upregulation of BCL3 by IL-6 in INA-6 and other MM cell lines. KCNA3 and TTC14 as the newly CRGs we discovered [36], their prognostic mechanism for MM needs to be further studied.

Complex cancer–immune system interactions generate both pro- and anti-tumor effects, and whose balance can be altered in favor of immunosuppressive environment, promoting a tumor progression in MM. These include impairment of antigen processing and presentation, T cell, NK and NKT cell, local recruitment, expansion and activation of immune suppressor cells like T regulatory cells (Tregs) and myeloid derived suppressor cells (MDSCs), differentiation of the protumoral tumor-associated macrophages and Th17 cells [37–40]. In terms of KEGG pathways, cell cycle, complement and coagulation cascades, neutrophil outer trap formation, nitrogen metabolism, p53 signaling pathway, ribosome, and spliceosome pathways were significantly enriched in the high-risk group. In the GO biological processes, antimicrobial humoral responses and responses against fungi were significantly enriched at high-risk group. In addition, immune score, ESTIMATE score and CYT were significantly different between high and low risk groups, while 13 immune cells such as memory B cells, plasma cells, T cells CD8 were also different in the two groups. These results suggest that the prognostic risk model constituted by CRGs may play an important role via GO and KEGG pathways in MM. Through the relationship analysis between risk score and immune score of 4 CRGS, we found that the risk model may affect the immune process of MM.

A total of nine drugs were predicted in our study, namely BI.D1870, SL.0101.1, KU.55933, AP.24534, Embelin, NSC.87877, JNK.9L, Axitinib and MK.2206. Among them, BI.D1870 has been proved to be a potential drug for multiple myeloma [41,42]. In addition, plasma cell leukemia is secondary to multiple myeloma, which is one of the late manifestations of multiple myeloma. AP.24534 has been proved to be a potential therapeutic drug for plasma cell leukemia [43]. Embelin is an inhibitor of X-linked inhibitor of apoptosis protein, which has anti-inflammatory and anticancer activities, and can inhibit the activation of signal transducer and activator of transcription 3 in multiple myeloma, so it is considered as a potential therapeutic drug [44]. Axitinib is a common drug for the treatment of renal cell carcinoma, but its application in multiple myeloma has not been found [45–47]. MK.2206 is a novel allosteric pan-Akt inhibitor, which has good anti-myeloma activity in MM cell lines and patient cells [48,49]. Ataxia Telangiectasia Mutated (ATM) can repair DNA and activate cell cycle checkpoints. Incorrect DNA damage repair will lead to genome instability [50]. ATM inhibition will increase the sensitivity of cancer cells to radiotherapy [51], so ATM inhibitors are one of the development directions of anticancer drugs, and KU-55933 is one of them. Studies have shown that ATM can increase the glucose uptake of cancer cells, and KU-55933 can inhibit tumor growth and metastasis in mice breast tumors by inhibiting GLUT1 translocation and vimentin expression [52]. As for NSC.87877, an effective bispecific protein phosphatase 26 (DUSP26) inhibitor, it has been found that it can inhibit the growth of neuroblastoma cell line and induce apoptosis [53]. Unfortunately, the reports about the potential anticancer effects of JNK.9L and L.0101.1 have not been found.

We develop a prognostic model of CRGs in MM. Cuproptosis is a novel type of cell death, it may bring up new solutions for the treatment of cancer. Moreover, risk score was a reliable independent prognostic factor for MM patients, and the nomogram shows the prediction effect of risk model. Additionally, the prognostic score focusing on CRGs’ expression signature showed a good performance in the prediction of MM survivorship, and other significant factors (stage, sample type, cluster and vital subtypes) with predictive values were considered in this study. Thus, the CDKN2A, BCL3, KCNA3 and TTC14 showed a robust prognostic effect for MM.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Key discipline construction project of Shanxi 1331 Project approved by Shanxi Education Department and Shanxi Finance Department: [Grant Number 2107 #3]; the 2019 Annual Science and Technology Innovation Team Project of Changzhi Medical College: [Grant Number CX201901]; the 2023 Shanxi Provincial Health Commission scientific research project: [Grant Number 2023036].