Skip to main content

The prognostic significance of ubiquitination-related genes in multiple myeloma by bioinformatics analysis

Abstract

Background

Immunoregulatory drugs regulate the ubiquitin-proteasome system, which is the main treatment for multiple myeloma (MM) at present. In this study, bioinformatics analysis was used to construct the risk model and evaluate the prognostic value of ubiquitination-related genes in MM.

Methods and results

The data on ubiquitination-related genes and MM samples were downloaded from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. The consistent cluster analysis and ESTIMATE algorithm were used to create distinct clusters. The MM prognostic risk model was constructed through single-factor and multiple-factor analysis. The ROC curve was plotted to compare the survival difference between high- and low-risk groups. The nomogram was used to validate the predictive capability of the risk model. A total of 87 ubiquitination-related genes were obtained, with 47 genes showing high expression in the MM group. According to the consistent cluster analysis, 4 clusters were determined. The immune infiltration, survival, and prognosis differed significantly among the 4 clusters. The tumor purity was higher in clusters 1 and 3 than in clusters 2 and 4, while the immune score and stromal score were lower in clusters 1 and 3. The proportion of B cells memory, plasma cells, and T cells CD4 naïve was the lowest in cluster 4. The model genes KLHL24, HERC6, USP3, TNIP1, and CISH were highly expressed in the high-risk group. AICAr and BMS.754,807 exhibited higher drug sensitivity in the low-risk group, whereas Bleomycin showed higher drug sensitivity in the high-risk group. The nomogram of the risk model demonstrated good efficacy in predicting the survival of MM patients using TCGA and GEO datasets.

Conclusions

The risk model constructed by ubiquitination-related genes can be effectively used to predict the prognosis of MM patients. KLHL24, HERC6, USP3, TNIP1, and CISH genes in MM warrant further investigation as therapeutic targets and to combat drug resistance.

Peer Review reports

Introduction

Multiple myeloma (MM) is a prevalent hematological malignancy characterized by the malignant proliferation of plasma cells in the bone marrow. This results in the secretion of monoclonal immunoglobulin or its fragments, leading to bone marrow failure, multiple bone destruction, and damage to various organs or tissues [1]. The number of MM patients is increasing worldwide [2]. The incidence is about 1.03% in China, and there are more male patients than female patients [3]. However, the etiology of MM is still unclear. Most researchers agree that abnormal plasmacytes originate from memory B lymphocytes or proplasmacytes with C-myc gene recombination and high expression levels of certain N-ras genes.This leads to unrestricted plasmacyte proliferation and abnormal increase of IL-6 in the bone marrow. The clinical manifestations include anemia, kidney impairment, hypercalcemia, and other symptoms, while severe cases may lead to death [4]. Therefore, exploring the genes associated with MM will help provide a theoretical basis for early diagnosis and treatment targets.

The ubiquitin-proteasome system (UPS) plays a critical role in the therapy of MM [5]. The inhibition of UPS has been considered a new strategy against MM. For example, the proteasome inhibitor (PI) can reduce the expression of TNF-α and NF-κB, induce apoptosis, and suppress drug efflux [6]. Immunomodulatory drugs are mainly applied to the treatment of MM, targeting the UPS. Immunomodulatory drugs induce selective ubiquitination and degradation of MM-associated lymphatic transcription factors IKZF1 and IKZF3 by binding to E3 ubiquitin ligase substrate receptor protein targets (CRBN) [7, 8]. Neuronally expressed developmentally downregulated 4 − 1 (NEDD4-1) inhibits the AKT signaling pathway by inducing proteasome degradation through ubiquitinating phosphorylated AKT-Ser473 and promoting apoptosis of MM cells. The PI Bortezomib can lead to the accumulation of polyubiquitylated proteins. Bortezomib promotes apoptosis of MM cells by inhibiting the 26 S proteasome [9]. Anti-MM drugs show their activity mainly by inhibiting ubiquitin-related enzymes, and target ubiquitin pathway promotes MM cell death. But drug resistance often develops [10]. The evidence showed that hyperactive small ubiquitin-like modifier is closely related to the progression of MM, while the higher the level of hyperactive small ubiquitin-like modifier, the poorer the survival [11, 12]. The research by Du et al. indicated that the ubiquitin receptor PSMD4/Rpn10 is effective in regulating cytotoxicity in MM, making it a potential therapeutic target [13]. According to the findings, it can be seen that ubiquitination is closely related to the development and treatment of MM. Although immunomodulatory drugs can prolong survival, there is still the possibility of relapse and drug resistance [14]. Moreover, PI showed serious side effects, causing cardiovascular damage [15]. Hence, there is an urgent need to develop new targets to treat MM and reduce drug resistance.

Methods

Download expression data and ubiquitination-related genes

We employed the TCGA database to download gene expression data and clinical information of MM as a training dataset of the risk model under the project ID MMRF-COMMPASS by using the R statistical software (version 4.2.1, https://mirrors.tuna.tsinghua.edu.cn/CRAN/bin/windows/base/R-4.2.1-win.exe) based on R package GDCRNATools (v1.16.2, https://github.com/Jialab-UCR/GDCRNATools, access time: 2023-6-20) [16]. A total of 859 samples were included in the TCGA MMRF-COMMPASS dataset. Among them, 764 MM samples were selected as a case group by PrimaryBloodDerivedCancer-BoneMarrow screening. Then, we used BoneMarrowNormal with project ID TARGET-AML downloaded from the leukemia dataset to screen normal bone marrow tissue samples due to the absence of normal samples in the MMRF-COMMPASS dataset. As a result, we got 20 normal samples for the control group. The expression data from 784 integrated samples were used for subsequent analysis. In the meanwhile, we obtained a gene expression dataset of MM samples with prognostic information GSE2658 [17] (549 tissue samples of MM at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE2658, Chip platform GPLGPL570), probe data and survival information necessary to verify the risk model. We searched for the keyword “ubiquitin” to acquire the gene set from the Signatures Database [18] v7.5.1 (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). A total of 197 C2 curated sets were obtained, and 1282 related genes were selected after merging data and de-duplicating for subsequent analysis.

Analysis of MM differentially expressed genes (DEGs) and differentially curated ubiquitination-related genes and identification of ubiquitination-related subtypes

The analysis of DEGs was performed based on MM expression data and normal bone marrow expression data by using the limma method [19] of the built-in gdcDEAnalysis function in the GDCRNATools package. In addition, we conducted an intersection between DEGs and curated ubiquitination-related genes. In order to conduct a more accurate model, individual differences in ubiquitination gene expression were excluded. We applied R package ConsensusClusterPlus (http://bioconductor.org/packages/release/bioc/html/ConsensusClusterPlus.html) based on the expression quantity of differentially curated ubiquitination genes to conduct consistent cluster analysis on the samples of the training dataset and set parameters: maxK = 10 (maximum number to evaluate), reps = 100 (number of subsamples), while the most appropriate cluster number was selected, the samples were divided into different subtypes, and the classification reliability was verified by PCA.

K-M survival analysis among different ubiquitination subtypes, analyzing tumor immune infiltration analysis, and comparison of clinical information differences

The survival differences between different ubiquitination subtypes were analyzed. The Estimate algorithm [20] was employed to perform immune scores for samples of each subtype and to compare the differences in immune scores. Moreover, CIBERSORT [21] was used to evaluate the proportion of immune cell infiltration and compare the difference in samples of each subtype.

Differences in checkpoint genes and HLA family genes between subtypes and representative gene collections

The list of immune checkpoint genes was obtained based on published article [22], and the list of HLA family genes was obtained by searching the HLA nomenclature database (http://hla.alleles.org/genes/index.html). Wilcoxon test was used to compare the expression of immune checkpoint genes and HLA family genes among subtypes.

For obtaining representative genes of subtypes, we used the limma method from the R package to conduct difference analysis between samples of a specific subtype and samples of other subtypes. Only genes highly expressed in the subtype (logFC > 1, P < 0.05) were selected as representative genes of the subtype.

Prognostic risk models were constructed based on single-factor and multiple-factor analysis [23] and scored after picking the union set de-duplicating representative genes for each subtype and collecting the survival information of each subtype.

The predictive power evaluation of the risk model based on the TCGA dataset

According to the median risk score, the samples were divided into a high-risk group and a low-risk group, and the survival difference between the high- and low-risk groups was compared respectively. Finally, the ROC curve was plotted.

The independent predictive power of the risk model was evaluated based on a nomogram and a calibration graph combined with the risk model score and clinical phenotype.

Risk model validation based on the GSE2658 dataset

According to the risk model, risk scores were assigned to the samples in the validation dataset, and the survival difference between high-risk and low-risk samples was compared using the median as the standard. The predictive power of the validation dataset is verified by drawing the nomogram and a calibration graph simultaneously.

Prediction of drug sensitivity and subtype characterization correlation analysis

We utilized the Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org/) database to estimate the sensitivity of each patient to chemotherapy drugs. The half maximal inhibitory concentration (IC50, where a lower IC50 value indicates greater cell sensitivity to the drug) was quantified using the pRRophetic package in the R language. The Wilcoxon test was used to compare drug sensitivity between high- and low-risk groups. The classification samples of different subtypes within the high- and low-risk groups were collected and analyzed.

External validation of model genes

The Kaplan-Meier Plotter was used to verify the correlation between model genes and MM survival. The P-value < 0.05 was considered statistically significant.

Results

DEGs analysis

A total of 784 samples were obtained after screening for subsequent analysis. In the meanwhile, we obtained 197 ubiquitination-related genes from the Signatures Database. A total of 6533 DEGs were identified using the limma method with |log2 (FC)| > 0.58 and P < 0.05 as the criterion through the built-in gdcDEAnalysis function of the GDCRNATools package (Fig. 1A). Furthermore, we got 87 differential ubiquitination-related genes. Among them, 47 genes were highly expressed in the tumor, and 40 genes were highly expressed in the control group, suggesting a close relationship between ubiquitination-related genes and MM (Supplementary Table 1, Fig. 1B). The expression levels of ubiquitination-related genes are presented in Supplementary Table 2. The results of the Wilcoxon test were statistically significant (Fig. 1C). The correlation analysis of ubiquitination-related DEGs was shown in Fig. 1D. The results suggest a strong correlation among highly expressed ubiquitination-related genes.

Fig. 1
figure 1

A DEGs volcano map of MM and normal bone marrow samples; B Ubiquitination-related DEGs volcano of MM and normal bone marrow samples; C There were significant differences in ubiquitination-related genes; D Correlation analysis of ubiquitination-related DEGs. *: 0.01 < P < 0.05, **: 0.001 < P < 0.01, ***: 0.0001 < P < 0.001, ****: P < 0.0001

Identification of molecular subtypes

According to the method, we performed a consistent cluster analysis of 764 TCGA MM samples based on the 84 ubiquitination-related genes obtained earlier. We selected the Delta area to show the inflection point (Fig.2A), and the clustering heat-map was well-structured (Fig. 2B), while the CDF curve was leveling off (Fig.2C). The corresponding k value was used as the number of clusters, so we finally determined that the number of clusters is 4, that is, the samples are divided into four subtypes. Principal Component Analysis (PCA) was conducted using the expression levels of ubiquitination genes, followed by ANOSIM analysis to compare thevariations among subtypes (Fig. 2D-E).

Fig. 2
figure 2

A Change of area of CDF curve under the number of clusters; B Clustering hierarchy diagram; C CDF curve for each cluster number; D PCA of differentially ubiquitination-related genes in 3D; E PCA Anosim in 2D, R > 0 refers there is difference between clusters, P < 0.05 refers the difference is statistically significant

K-M survival analysis between different immune subtypes, analyze immune cell infiltration, and correlate with clinical data

Based on the survival information of samples in each subtype, we utilized the R packages survival and survminer to compare the survival conditions of different subtypes.Our analysis revealed significant differences in survival between subtypes (Fig. 3A). Among them, cluster 1 had a poor prognosis, while the rest had a better prognosis. At the same time, the results of enrichment scores showed that the score of cluster 1 was significantly lower than that of cluster 2 and cluster 4 (Fig. 3B), which illustrated that ubiquitination-related genes were related to the prognosis of MM patients. In addition, the immune scores were performed by using the ESTIMATE algorithm, and we compared the differences in immune scores among subtypes (Fig. 3C). As a result, the differences in immune scores between the subtypes were statistically significant. The tumor purity of cluster 1 and cluster 3 was significantly higher than those the other two clusters. The immune score and stromal score in cluster 2 and cluster 4 were higher than in cluster 1 and cluster 3. The GSVA score of each cluster was consistent with the ESTIMATE score, which confirmed the correlation between ubiquitination-related genes and MM pathogenic progression. The proportion of immune cells in tumor tissues, evaluated by the CIBERSORT algorithm showed that the differences in CIBERSORT fraction in clusters were statistically significant (Fig. 3D). In terms of immune cells, B cells memory in cluster 1 accounted for the highest proportion, while the proportion of plasma cells in cluster 1 was less than in cluster 2. Additionally, T cells CD8 accounted for the highest proportion in cluster 4.

Fig. 3
figure 3

A The differences in survival analysis in ubiquitination-related gene clusters; B The differences in enrichment scores in ubiquitination-related gene clusters; C The differences of immune infiltration scores in ubiquitination-related gene clusters; D The differences in the proportion of immune cells in ubiquitination-related gene clusters. Wilcoxon Test *: 0.01 < P < 0.05, **: 0.001 < P < 0.01, ***: 0.0001 < P < 0.001, ****: P < 0.0001

The collection of representative genes in each cluster

The expression levels of checkpoint genes and HLA family genes are presented in Supplementary Table 3. The results revealed significant differences in checkpoint genes and HLA family genes among various molecular subtypes, indicating notable distinctions significant differences between the subtypes. For example, the expression level of BTN3A1 in cluster 4 was lower than in other clusters, and the expression level of HLA-A, B, and C in clusters 1 and 4 was higher than in clusters 2 and 3 (Fig. 4A). Moreover, we utilized the limma package to analyze the differences between the samples of one subtype and other subtypes, and selected the highly expressed genes as representative genes of each subtype (Fig. 4B). The intersection of representative genes in subtypes was limited, meaning that the differences in representative genes of each subtype were significant. In the meanwhile, we selected all sets of ubiquitination-related genes in the Signatures Database to intersect with all subtypes, and we got 72 ubiquitination representative genes for subsequent survival analysis (Fig. 4C). Then, we constructed a risk model combined with the survival information of the samples. The reason LASSO regression only considers survival and death time, and does not include survival information, is that we need to select genes that are significantly related to survival time and survival state in the single-factor analysis to ensure the reliability of the model. Based on single-factor analysis, the LASSO regression was performed on genes with P < 0.005 (Fig. 4D).

Fig. 4
figure 4

A The differences between checkpoint genes and HLA family genes in subtypes. (Wilcoxon Test *: 0.01 < P < 0.05, **: 0.001 < P < 0.01, ***: 0.0001 < P < 0.001, ****: P < 0.0001); B The representative genes of clusters (The red dots represent the representative genes of each cluster, log2FC > 1 and P < 0.05); The veen diagram of representative genes in each cluster; LASSO regression

The risk model construction and validation based on the TCGA dataset

The risk score model was constructed according to the results of LASSO model: Riskscore = 0.0046833357739203*SOCS3 + 0.0296897767259966*UBE2T + 0.0185158927170784*TOP2A + (-0.0189779895672147)*HERC6 + (-0.0261853093234728)*USP3 + (-0.0124611715684952)*TNIP1 + (-0.0668058378550144)*KLHL24 + 0.0114569251299551*SDE2 + 0.0439772194879231*VCPIP1 + 0.0270466928975837*YOD1 + (-0.00796472303046347)*CISH”.

After scoring the TCGA MM samples, the samples were divided into a high-risk group and a low-risk group according to the median risk score. The differences in gene expression (Fig. 5A) and survival probability (Fig. 5B) between high- and low-risk groups were statistically significant, while the ROC curve had a high interpretive ability (Fig. 5C). The model genes UBE2T, TOP2A, VCPIP1, SDE2, YOD1, and SOCS3 were highly expressed in the low-risk group, while KLHL24, HERC6, USP3, TNIP1, and CISH were highly expressed in the high-risk group. Based on the prognostic model and TCGA dataset, we conducted univariate and multivariate analyses with independent prognostic clinical information and risk score, and found that age and risk score were significantly correlated with prognostic survival. According to the age and risk score, the risk model nomogram (Fig. 5D) and calibration graphs of 1-, 2-, and 3-year survival rates were conducted (Fig. 5E).

Fig. 5
figure 5

A The heat map of gene expression in low- and high-risk groups; The survival probability of low- and high-risk groups; C The ROC curve of the risk model in the TCGA dataset; D Nomogram of risk model; E The calibration graph of 1-, 2- and 3-year survival rate

The validation of the risk model based on the GSE2658 dataset

Based on the expression of model genes in the validation dataset GSE2658, we scored the risk of MM samples (Supplementary Table 4), and divided them into high- and low-risk groups according to the median risk score. According to the expression of model genes in the high- and low-risk groups (Fig. 6A), the model genes UBE2T, TOP2A, VCPIP1, SDE2, YOD1, and SOCS3 were highly expressed in the low-risk group, while KLHL24, HERC6, USP3, TNIP1, and CISH were highly expressed in the high-risk group. In addition, survival analysis was conducted for the samples in high- and low-risk groups (Fig. 6B). The survival probability was higher in the low-risk group than in the high-risk group, and the difference was statistically significant. The ROC curve indicated that the risk model effectively predicted the prognosis of MM (Fig. 6C).

The results of single-factor analysis indicated that the risk score was significantly correlated with survival time (Supplementary Table 5). The higher the risk score, the shorter the survival time. The nomogram (Fig. 6D) and calibration graphs of 1-, 2- and 3-year survival rates (Fig. 6E) were plotted.

Fig. 6
figure 6

The expression of model genes in high- and low-risk groups based on the GEO validation dataset; B The survival analysis of model genes in high- and low-risk groups based on the GEO validation dataset; C The ROC of risk model prediction in high- and low-risk groups based on validation dataset; The nomogram of validation dataset; E The calibration graphs of 1-, 2- and 3-year survival model based on GEO dataset

The prediction of drug sensitivity and subtype characterization correlation analysis

Based on the TCGA MM sample expression dataset and GDSC database to estimate the sensitivity to chemotherapy drugs for each patient, and compare the drug sensitivity between high- and low-risk groups. Most drugs had significantly different sensitivity in high- and low-risk groups. For instance, 5-aminoimidazole-4-carboxamide-1-β-riboside (AICAr) and BMS.754,807 had higher drug sensitivity in a low-risk group than in a high-risk group, while Bleomycin had higher drug sensitivity in the high-risk group (Fig. 7A). We analyzed the enrichment of high- and low-risk samples in various clusters in the TCGA dataset using the fisher test, and the results of the enrichment are displayed in (Fig. 7B). High-risk samples were significantly enriched in cluster 1, while low-risk samples were significantly enriched in cluster 2 and cluster 4.

Fig. 7
figure 7

A The prediction of drug sensitivity in the TCGA dataset and the difference in high- and low-risk groups; B The enrichment analysis of high- and low-risk groups in the TCGA dataset

External validation of model genes

The validation of these 11 model genes on the Kaplan-Meier Plotter database indicated that the expression of the model genes was significantly correlated with the survival of MM patients (P < 0.05, Fig. 8).

Fig. 8
figure 8

Correlation between the expression of model genes and survival of MM patients

Discussion

Protein modification plays a crucial role in the regulating of cellular processes, in which protein ubiquitination can alter tumor metabolism and immune regulation [24]. Ubiquitination is related to the stability of oncoproteins in MM [25]. Nowadays, the proteasome inhibitor has been widely used in treatment of MM [11, 26]. However, drug resistance is one of the major challenges we face. And the correlation between ubiquitination-related genes and MM has not been thoroughly examined.

According to the DEGs analysis, a total of 47 ubiquitination-related genes were highly expressed in the MM group, while 40 genes were highly expressed in the control group, indicating a connection between ubiquitination-related genes and MM. Four clusters exhibited significant differences according to the consistent cluster analysis. At the same time, the survival information was significantly among different clusters. The prognosis of cluster 1 was poor, and the enrichment score was significantly lower than that of cluster 2 and cluster 4. It is suggested that the genes related to to ubiquitination are associated with the prognosis of MM patients. Patients with poor survival may experience early death. Developing targeted therapies that act on genes associated with aggressive disease or drug resistance can be applied to improve outcomes for MM patients with high-risk [27]. The tumor purity in clusters 1 and 3 was higher than that in clusters 2 and 4 consistent with the trend of GSAV score.

In this study, the analysis of gene expression levels of using TCGA and GEO datasets revealed that the that model genes UBE2T, TOP2A, VCPIP1, SDE2, YOD1, and SOCS3 were highly expressed in the low-risk group, whereas KLHL24, HERC6, USP3, TNIP1, and CISH were highly expressed in the high-risk group. Regarding genes highly expressed in high-risk groups, KLHL24 is the main substrate adaptor protein of Cullin3-RING ligase (CRL3) which is one of the most common E3 ubiquitin ligases. CRL3 can affect the stability of functional proteins by mediating substrate ubiquitination modifications. Dysregulation of CRL3 can lead to the development of multiple diseases. KLHL24 inactivation is associated with hypertrophic cardiomyopathy [28, 29]. KLHL24 inhibits the activation of fibroblasts, leading to skin fibrosis, and hinders skin wound healing [30]. People who carry a mutation in the KLHL24 gene are at risk of developing epidermolysis bullosa [31]. HERC6 is also an E3 ubiquitin ligase, primarily expressed in the testicles and fetal brain, with rare expression in the heart and skeletal muscle. HERC6 participates in a variety of cellular activities, including cell proliferation, cell migration, and neurodevelopment. HERC6 can be cancer-promoting factors and tumor suppressor genes, depending on the type of cancer [32]. At the same time, the up-regulation of HERC6 is present in systemic lupus erythematosus (SLE) and promotes inflammation [33]. Ubiquitin-specific protease 3 (USP3) is highly expressed in multiple cancers. It plays an important role in tumor proliferation and invasion. The overexpression of USP3 leads to an unfavorable prognosis in breast cancer patients and stomach cancer metastasis [34]. Furthermore, the activation of USP3 induces neuroblastoma [35]. The tumor necrosis factor α-induced protein 3-interacting protein 1 (TNIP1) plays a part in mitophagy [36]. It has an impact on the development of autoimmune diseases [37], such as lupus nephritis [38]. The overexpression of cytokine-inducible SH2-containing protein (CISH) promotes inflammation in the elderly [39]. The knockdown of CISH may regulate the metabolic activity of NK cells to exert an anti-tumor effect [40, 41]. In addition, a previous study has shown that the inhibiting of CISH could improve the outcomes of immune checkpoint blockade therapy [42]. The existing study has confirmed that HERC4 is involved in inhibiting the proliferation of myeloma cells [43]. More important, model genes were verified to be significantly correlated with the survival of MM patients through external validation. The role of KLHL24, HERC6, USP3, TNIP1, and CISH genes in MM has not been reported in the literature, so it is worth further investigation in vitro and in vivo.

Effector lymphocyte dysfunction and malignant plasma cells are associated with MM which is accompanied by immunosuppression [44]. Immune escape and loss of antigenicity appear in MM. Tumor cells in patients with MM exhibit increased expression of the immune checkpoint receptor programmed death receptor ligand that promotes immune escape [45]. The proportions of B cells memory, plasma cells, T cells CD4 naïve, and T cells CD8 showed significant differences between the clusters in this study. In addition, the differences in ICGs and HLA genes between the clusters were statistically significant, except for BTN2A2, HLA-DQA1, and HLA-DRB5. After the intersection, a total of 72 ubiquitination representative genes were screened and divided into high- and low-risk groups for survival analysis.

In the present study, AICAr exhibited higher drug sensitivity in a low-risk group. It has been reported the effects in accelerating apoptosis [46] and inhibits the growth of MM cells [47]. Insulin-like growth factor 1 receptor/insulin receptor family kinases inhibitor BMS.754,807 also showed higher drug sensitivity in the low-risk group in this study. The evidence showed that it is effective in inhibiting MM cells [48]. Bleomycin demonstrated higher drug sensitivity in the high-risk group in the present study. It could be used in drug-resistant MM patients, which plays a role in the treatment of drug-resistant MM [49]. In Fig. 2E, cluser4 is clearly distinguished by PCA, but cluster1 \ 2 \ 3 has significant overlap. Exploring additional biological features related to gene clusters of these different clustering patterns, in order to identify patient sub clusters and develop therapies in clinical, will be a highly promising work in the future. Although this study identified ubiquitination genes associated with multiple myeloma, But differences in the transcriptome of hematopoietic system tumors may result in variations in drug efficacy among different samples and developmental stages [50]. Additionally, this study lacks analysis based on ubiquitination-related genes, risk scores, and drug sensitivity across samples at different developmental stages, which is a potential research point.

The risk model had exhibited satisfactory predictive power. It was elucidated that there was poor survival in the high-risk group. Furthermore, genes related to ubiquitination were initially explored in connection with MM in our study. Ubiquitination has a significant impact on celluar protein interactions and is linked to cancer progression. Due to the wide range of viral and non viral etiologies in MM [51, 52], clustering based on ubiquitination can better distinguish the differences in etiology of multiple myeloma, identify different interaction patterns between these clusters, and provide personalized treatment for patients based on these etiologies.

Whereas, there were some limitations in the present study. First, we did not find normal bone marrow samples in the same dataset but screened control samples in the leukemia dataset, and there is a lack of ethnic details in clinical data, which may have led to a decrease in reliability. Furthermore, we lacked survival information other than survival time and death time for LASSO regression analysis to enhance the effectiveness of the risk model. Due to the fact that only a very small number of enriched immune cells show insignificant differences in immune scores across different clusters, validation of all immune cells is necessary. However, due to limitations in research conditions, we are unable to conduct cell experiments for validation. This is a limitation of the study and we hope to verify it in future research.

Conclusion

The risk model is constructed to predict the prognosis of MM patients. Overexpression of ubiquitination-related genes such as KLHL24, HERC6, USP3, TNIP1, and CISH may indicated a poor prognosis and lower survival rate in MM patients. These genes have the potential as therapeutic targets for MM. The mechanisms of these genes in drug sensitivity and MM pathogenesis deserve further study.

Availability of data and materials

The data in this study are available with permission from the corresponding author.

References

  1. Cowan AJ, et al. Diagnosis and management of multiple myeloma: a review. JAMA. 2022;327(5):464–77.

    Article  CAS  PubMed  Google Scholar 

  2. Padala SA, et al. Epidemiology, staging, and management of multiple Myeloma. Med Sci (Basel). 2021;9(1):3.

    CAS  PubMed  Google Scholar 

  3. Liu J, et al. Union for China Lymphoma Investigators of the Chinese Society of Clinical Oncology; Union for China Leukemia Investigators of the Chinese Society of Clinical Oncology. Incidence and mortality of multiple myeloma in China, 2006-2016: an analysis of the Global Burden of Disease Study 2016. J Hematol Oncol. 2019;12(1):136.

  4. van de Donk N, Pawlyn C, Yong KL. Multiple myeloma. Lancet. 2021;397(10272):410–27.

    Article  PubMed  Google Scholar 

  5. Eichner R, et al. Cross Talk networks of mammalian target of Rapamycin Signaling with the Ubiquitin Proteasome System and their clinical implications in multiple myeloma. Int Rev Cell Mol Biol. 2019;343:219–97.

    Article  CAS  PubMed  Google Scholar 

  6. Narayanan S, et al. Targeting the ubiquitin-proteasome pathway to overcome anti-cancer drug resistance. Drug Resist Updat. 2020;48:100663.

    Article  PubMed  Google Scholar 

  7. Krönke J, et al. Lenalidomide causes selective degradation of IKZF1 and IKZF3 in multiple myeloma cells. Science. 2014;343(6168):301–5.

    Article  PubMed  Google Scholar 

  8. Holstein SA, McCarthy PL. Immunomodulatory drugs in multiple myeloma: mechanisms of action and clinical experience. Drugs. 2017;77(5):505–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Park J, Cho J, Song EJ. Ubiquitin-proteasome system (UPS) as a target for anticancer treatment. Arch Pharm Res. 2020;43(11):1144–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Beider K, et al. Blocking of transient receptor potential vanilloid 1 (TRPV1) promotes terminal mitophagy in multiple myeloma, disturbing calcium homeostasis and targeting ubiquitin pathway and bortezomib-induced unfolded protein response. J Hematol Oncol. 2020;13(1):158.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Heynen G, et al. SUMOylation inhibition overcomes proteasome inhibitor resistance in multiple myeloma. Blood Adv. 2023;7(4):469–81.

    Article  CAS  PubMed  Google Scholar 

  12. Du L, Liu W, Rosen ST. Targeting SUMOylation in cancer. Curr Opin Oncol. 2021;33(5):520–5.

    Article  CAS  PubMed  Google Scholar 

  13. Du T, et al. Ubiquitin receptor PSMD4/Rpn10 is a novel therapeutic target in multiple myeloma. Blood. 2023;141(21):2599–614.

    CAS  PubMed  PubMed Central  Google Scholar 

  14. Lub S, et al. Novel strategies to target the ubiquitin proteasome system in multiple myeloma. Oncotarget. 2016;7(6):6521–37.

    Article  PubMed  Google Scholar 

  15. Zheng Y, et al. Cardiovascular toxicity of Proteasome inhibitors in multiple myeloma therapy. Curr Probl Cardiol. 2023;48(3):101536.

    Article  PubMed  Google Scholar 

  16. Li R, et al. GDCRNATools: an R/Bioconductor package for integrative analysis of lncRNA, miRNA and mRNA data in GDC. Bioinformatics. 2018;34(14):2515–7.

    Article  PubMed  Google Scholar 

  17. Hanamura I, et al. Prognostic value of cyclin D2 mRNA expression in newly diagnosed multiple myeloma treated with high-dose chemotherapy and tandem autologous stem cell transplantations. Leukemia. 2006;20(7):1288–90.

    Article  CAS  PubMed  Google Scholar 

  18. Liberzon A, et al. The Molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1(6):417–25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Ritchie ME, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Yoshihara K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.

    Article  PubMed  Google Scholar 

  21. Chen B, et al. Profiling Tumor infiltrating Immune cells with CIBERSORT. Methods Mol Biol. 2018;1711:243–59.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Yu Y, et al. Immune Checkpoint Gene expression profiling identifies programmed cell death Ligand-1 centered immunologic subtypes of oral and squamous cell Carcinoma with favorable survival. Front Med (Lausanne). 2021;8:759605.

    Article  PubMed  Google Scholar 

  23. Liu Q, et al. An Individualized Prognostic signature for Clinically Predicting the survival of patients with bladder Cancer. Front Genet. 2022;13:837301.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Pellegrino NE, et al. The Next Frontier: Translational Development of Ubiquitination, SUMOylation, and NEDDylation in Cancer. Int J Mol Sci. 2022;23(7):3480.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Paulmann C, et al. The OTUD6B-LIN28B-MYC axis determines the proliferative state in multiple myeloma. EMBO J. 2022;41(20):e110871.

  26. Accardi F, et al. The Proteasome and Myeloma-Associated Bone Disease. Calcif Tissue Int. 2018;102(2):210–26.

    Article  CAS  PubMed  Google Scholar 

  27. Hanamura I. Multiple myeloma with high-risk cytogenetics and its treatment approach. Int J Hematol. 2022;115(6):762–77.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Hedberg-Oldfors C, et al. Cardiomyopathy with lethal arrhythmias associated with inactivation of KLHL24. Hum Mol Genet. 2019;28(11):1919–29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Walsh R, et al. Minor hypertrophic cardiomyopathy genes, major insights into the genetics of cardiomyopathies. Nat Rev Cardiol. 2022;19(3):151–67.

    Article  PubMed  Google Scholar 

  30. Liu Y, et al. Excess KLHL24 impairs skin Wound Healing through the degradation of Vimentin. J Invest Dermatol. 2023;143(7):1289–e129815.

    Article  CAS  PubMed  Google Scholar 

  31. Lin Z, et al. Stabilizing mutations of KLHL24 ubiquitin ligase cause loss of keratin 14 and human skin fragility. Nat Genet. 2016;48(12):1508–16.

    Article  PubMed  Google Scholar 

  32. Sala-Gaston J, et al. HERC ubiquitin ligases in cancer. Cancers (Basel). 2020;12(6):1653.

    Article  CAS  PubMed  Google Scholar 

  33. Cao L, et al. HERC6 is upregulated in peripheral blood mononuclear cells of patients with systemic lupus erythematosus and promotes the disease progression. Autoimmunity. 2022;55(8):506–14.

    Article  CAS  PubMed  Google Scholar 

  34. Wu X, et al. USP3 promotes gastric cancer progression and metastasis by deubiquitination-dependent COL9A3/COL6A5 stabilisation. Cell Death Dis. 2021;13(1):10.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Nagy Z, et al. An ALYREF-MYCN coactivator complex drives neuroblastoma tumorigenesis through effects on USP3 and MYCN stability. Nat Commun. 2021;12(1):1881.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Le Guerroué F, et al. TNIP1 inhibits selective autophagy via bipartite interaction with LC3/GABARAP and TAX1BP1. Mol Cell. 2023;83(6):927–e9418.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Shamilov R, Aneskievich BJ. TNIP1 in autoimmune Diseases: regulation of toll-like receptor signaling. J Immunol Res. 2018:3491269.

  38. Brady MP, et al. TNIP1/ABIN1 and lupus nephritis: review. Lupus Sci Med. 2020;7(1):e000437.

  39. Jin J, et al. CISH impairs lysosomal function in activated T cells resulting in mitochondrial DNA release and inflammaging. Nat Aging. 2023;3(5):600–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Zhu H, et al. Metabolic reprograming via deletion of CISH in human iPSC-Derived NK cells promotes in vivo persistence and enhances anti-tumor activity. Cell Stem Cell. 2020;27(2):224–e2376.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Bernard PL, et al. Targeting CISH enhances natural cytotoxicity receptor signaling and reduces NK cell exhaustion to improve solid tumor immunity. J Immunother Cancer. 2022;10(5):e004244.

  42. Kumar S, et al. Epitranscriptomic Approach: to improve the efficacy of ICB Therapy by Co-targeting Intracell Checkp CISH. Cells. 2021;10(9):2250.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Mao X, et al. The emerging roles of the HERC Ubiquitin ligases in Cancer. Curr Pharm Des. 2018;24(15):1676–81.

    Article  CAS  PubMed  Google Scholar 

  44. Nakamura K, Smyth MJ, Martinet L. Cancer immunoediting and immune dysregulation in multiple myeloma. Blood. 2020;136(24):2731–40.

    Article  CAS  PubMed  Google Scholar 

  45. Swamydas M, et al. Deciphering mechanisms of immune escape to inform immunotherapeutic strategies in multiple myeloma. J Hematol Oncol. 2022;15(1):17.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Bardeleben C, et al. Metabolomics identifies pyrimidine starvation as the mechanism of 5-aminoimidazole-4-carboxamide-1-β-riboside-induced apoptosis in multiple myeloma cells. Mol Cancer Ther. 2013;12(7):1310–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Baumann P, et al. Activation of adenosine monophosphate activated protein kinase inhibits growth of multiple myeloma cells. Exp Cell Res. 2007;313(16):3592–603.

    Article  CAS  PubMed  Google Scholar 

  48. Carboni JM, et al. BMS-754807, a small molecule inhibitor of insulin-like growth factor-1R/IR. Mol Cancer Ther. 2009;8(12):3341–9.

    Article  CAS  PubMed  Google Scholar 

  49. Bennett JM, et al. Phase II study of adriamycin and bleomycin in patients with multiple myeloma. Cancer Treat Rep. 1978;62(9):1367–9.

    CAS  PubMed  Google Scholar 

  50. Mukherjee S, et al. Silico Integration of Transcriptome and Interactome predicts an ETP-ALL-Specific transcriptional footprint that decodes its Developmental Propensity. Front Cell Dev Biol. 2022;10:899752.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Kar A, et al. The HBV web: an insight into molecular interactomes between the hepatitis B virus and its host en route to hepatocellular carcinoma. J Med Virol. 2023;95(1):e28436.

    Article  CAS  PubMed  Google Scholar 

  52. Marsh DJ, Ma Y, Dickson KA. Histone monoubiquitination in chromatin remodelling: focus on the histone H2B interactome and Cancer. Cancers (Basel). 2020;12(11):3462.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This study was funded by the Yunnan Science and Technology Planning Project, the Kunming Medical University Joint Special Project (202201AY070001-201); the Kunming Health Science and Technology Personnel Training Project and the “Ten Hundred Thousand” Project (2020-SW-09).

Author information

Authors and Affiliations

Authors

Contributions

F.Z. and X.C. designed the study concept, F.Z., X.C., and Q.P. wrote the main manuscript text, H.W. and T.G. collected data, and J.Y. and Z.J. analyzed data. All authors reviewed the manuscript.

Corresponding author

Correspondence to Feng zhang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

zhang, F., Chen, XL., Wang, HF. et al. The prognostic significance of ubiquitination-related genes in multiple myeloma by bioinformatics analysis. BMC Med Genomics 17, 164 (2024). https://doi.org/10.1186/s12920-024-01937-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12920-024-01937-0

Keywords