Skip to main content

ORIGINAL RESEARCH article

Front. Pharmacol., 17 August 2022
Sec. Pharmacology of Anti-Cancer Drugs
This article is part of the Research Topic Anti-cancer Effects of Botanical Drugs Targeting the Tumor Microenvironment View all 15 articles

A ubiquitin-related gene signature for predicting prognosis and constructing molecular subtypes in osteosarcoma

Nan Wei,&#x;Nan Wei1,2Gong Chao-yang,&#x;Gong Chao-yang1,2Zhou Wen-ming,&#x;Zhou Wen-ming1,2Lei Ze-yuan,&#x;Lei Ze-yuan1,2Shi Yong-qiang,Shi Yong-qiang1,2Zhang Shun-bai,Zhang Shun-bai1,2Zhang Kai,Zhang Kai1,2Ma Yan-chao,
Ma Yan-chao1,2*Zhang Hai-hong,
Zhang Hai-hong1,2*
  • 1Orthopaedics Key Laboratory of Gansu Province, Lanzhou, China
  • 2Lanzhou University Second Hospital, Lanzhou, China

Background: Ubiquitination is medicated by three classes of enzymes and has been proven to involve in multiple cancer biological processes. Moreover, dysregulation of ubiquitination has received a growing body of attention in osteosarcoma (OS) tumorigenesis and treatment. Therefore, our study aimed to identify a ubiquitin-related gene signature for predicting prognosis and immune landscape and constructing OS molecular subtypes.

Methods: Therapeutically Applicable Research to Generate Effective Treatments (TARGET) was regarded as the training set through univariate Cox regression, Lasso Cox regression, and multivariate Cox regression. The GSE21257 and GSE39055 served as the validation set to verify the predictive value of the signature. CIBERSORT was performed to show immune infiltration and the immune microenvironment. The NMF algorithm was used to construct OS molecular subtypes.

Results: In this study, we developed a ubiquitin-related gene signature including seven genes (UBE2L3, CORO6, DCAF8, DNAI1, FBXL5, UHRF2, and WDR53), and the gene signature had a good performance in predicting prognosis for OS patients (AUC values at 1/3/5 years were 0.957, 0.890, and 0.919). Multivariate Cox regression indicated that the risk score model and prognosis stage were also independent prognostic prediction factors. Moreover, analyses of immune cells and immune-related functions showed a significant difference in different risk score groups and the three clusters. The drug sensitivity suggested that IC50 of proteasome inhibitor (MG-132) showed a notable significance between the risk score groups (p < 0.05). Through the NMF algorithm, we obtained the three clusters, and cluster 3 showed better survival outcomes. The expression of ubiquitin-related genes (CORO6, UBE2L3, FBXL5, DNAI1, and DCAF8) showed an obvious significance in normal and osteosarcoma tissues.

Conclusion: We developed a novel ubiquitin-related gene signature which showed better predictive prognostic ability for OS and provided additional information on chemotherapy and immunotherapy. The OS molecular subtypes would also give a useful guide for individualized therapy.

Introduction

Osteosarcoma (OS) serves as the most common primary solid malignancy of bone in children and adolescents. Its characteristic is higher malignant potential and occurrence of metastases (Botter et al., 2014; Kansara et al., 2014). At present, the combination of surgery and adjuvant chemotherapy is considered a potential therapeutic approach and has also improved the overall survival of OS (Harrison et al., 2018). However, gradually increasing metastatic and relapsed OS has been a new threat to the prognosis of OS (Meazza and Scanagatta, 2016). The study reported that the long-term survival rate of metastatic OS was <25% (Anwar et al., 2020). Hence, it is crucial and necessary to develop more effective novel agents and novel therapeutic targets.

Post-translational modification (PTM) is an important step to regulate various kinds of cellular processes, including phosphorylation, ubiquitination, acetylation, and other modifications (Chen et al., 2020). Amongst these PTMs, dysregulation of ubiquitination has been involved in multiple cancers (Nakayama and Nakayama, 2006). Generally, ubiquitination is medicated by three classes of enzymes, such as ubiquitin-activating enzymes (E1), ubiquitin-conjugating enzymes (E2), and ubiquitin ligases (E3) (Hershko and Ciechanover, 1998). Meanwhile, dysregulation of ubiquitination has received a growing body of attention in OS tumorigenesis and treatment. For example, Zhou et al. (2020) reported that the E3 ubiquitin ligase tripartite motif 7 (TRIM7) regulates metastasis and chemoresistance in OS through ubiquitination of breast cancer metastasis suppressor 1 (BRMS1). In addition, Zhang et al. (2017) verified that ubiquitin-conjugating enzyme E2 variant 1A (UEV1A) played a suppressive role in OS differentiation by promoting Smurf1-mediated Smad1 ubiquitination. Taken together, increasing evidence has shown that ubiquitination is closely associated with OS initiation and progression. However, the studies on ubiquitin-related gene signatures and ubiquitin-related molecular subtypes for OS still remain limited.

In this study, we constructed a ubiquitin-related gene signature including seven ubiquitin-related genes through the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) and Gene Expression Omnibus (GEO) database. Meanwhile, we confirmed that the risk score model has good performance in predicting the overall survival and was an independent prognostic factor. Based on the risk score, we evaluated the immune infiltration and chemotherapy drug sensitivity between high- and low-risk groups. Importantly, we also developed novel molecular subtypes for OS by prognostically significant ubiquitin-related genes.

Materials and methods

Data processing

The study included three datasets (Table 1). RNA sequence data and clinical follow-up information of the training set were downloaded from public databases: the Therapeutically Applicable Research to Generate Effective Treatments (TARGET; https://ocg.cancer.gov/programs/target). Then, FPKM values were transformed into TPM values through the R package “dplyr” (Mangiola et al., 2021). Gene expression matrix and clinical follow-up information of validation set (GSE21257 and GSE39055) were downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/). The workflow of the study is shown in Figure 1. The code of analysis was showed in Supplementary Data Sheet 1.

TABLE 1
www.frontiersin.org

TABLE 1. Clinical characteristics for all patients.

FIGURE 1
www.frontiersin.org

FIGURE 1. The flow chart for this study.

Construction and evaluation of a ubiquitin-related gene signature

We acquired 911 ubiquitin-related genes from the integrated annotations for Ubiquitin and Ubiquitin-like Conjugation Database (iUUCD, http://iuucd.biocuckoo.org/browse.php). Finally, we obtained the expression matrix of 895 ubiquitin-related genes from the TARGET dataset. Then, univariate Cox regression and Lasso Cox regression were used to select prognosis-related ubiquitin-related genes by using R packages “survival” (Lorent et al., 2014) and “glmnet” (Engebretsen and Bohlin, 2019). Multivariate Cox regression was performed to obtain the risk score and risk classification of each OS case. The algorithm formula was as follows: Risk score = (Exp of gene1 x Coefficient of gene 1)+(Exp of gene 2 x Coefficient of gene 2)+…+(Exp of gene N x Coefficient of gene N).

The receiver-operating characteristic (ROC) curves were used to evaluate the diagnostic efficacy of the risk score model through the “timeROC” (Blanche et al., 2013) software package. Based on a ubiquitin-related gene signature, the principal component analysis (PCA) was applied to show distribution patterns of OS case by using the “ggplot2” software package (Gómez-Rubio, 2017). The Kaplan–Meier curve analysis was performed to evaluate overall survival in the high- and low-risk group by using R package “survival” package. Moreover, multivariate Cox regression was applied to evaluate the independent prognostic value through R package “forestplot”.

Construction nomogram

In order to predict overall survival probability at 1, 3, and 5 years for osteosarcoma patients in the Target cohort, we constructed a nomogram to evaluate the probable overall survival of the OS patients through the R package “rms” (Bandos et al., 2009).

Immune cell infiltration and immune microenvironment

The CIBERSORT (Newman et al., 2019) was applied to assess the difference in the expression matrix of immune cells in the high- and low-risk group. Meanwhile, we compared the immune microenvironment (stromal score, immune score, and ESTIMATE score) by using the R package “estimate” (Yoshihara et al., 2013) in Cluster 1, Cluster 2, and Cluster 3.

Gene set enrichment analysis, Gene Ontology, and Kyoto Encyclopedia of Genes and Genomes

GSEA was performed to evaluate the significant biological functions between high- and low-risk groups through the R package “clusterProfiler” (Yu et al., 2012), and the significant enrichment was selected with | NES | > 1 and p value < 0.05. In addition, the annotated gene sets (“c5. go.v7.4. symbols.gmt” and “c2. cp.kegg.v7.4. symbols.gmt”) served as the reference set. The prognosis-related ubiquitin-related genes were used to conduct GO and KEGG enrichment analyses by using the R package “clusterProfiler.” The significant criterion was p value < 0.05.

Prediction of drug sensitivity and nonnegative matrix factorization clustering

In order to further explore the clinical reliability of the risk score model, we predicted the half-maximal inhibitory concentration (IC50) of chemotherapy drugs through R package “pRRophetic” (Geeleher et al., 2014). Moreover, the “NMF” method was applied to identify OS subtypes through the R package “NMF” (Lee and Seung, 1999).

Quantitative real-time PCR

A total of six osteosarcoma and six normal tissues were obtained from 12 patients at the Department of Orthopaedics, Lanzhou University Second Hospital in June 2022, and this research was supported by the Independent Ethics Committee (IEC). Total RNA from the tissues was extracted by the TRIzol (AG21101, Hunan, China). Reverse transcription was carried out by using the Evo M-MLV RT Kit (AG11705, Hunan, China). cDNA amplification was carried out by using SYBR® Green Premix Pro Taq HS qPCR Kit (AG11718, Hunan, China). GAPDH was used as an internal control. All gene primers were listed in Supplementary Table 2.

Statistical analysis

All statistical analyses and figure visualization were used by R software (version 4.0.2) in this study. Perl (version 5.8.3) was applied to combine gene expression data with clinical information. The p value < 0.05 was considered statistical significance.

Results

Construction of a ubiquitin-related gene signature in osteosarcoma

In order to construct a ubiquitin-related gene signature, we firstly identified 42 ubiquitin-related genes that were significantly associated with the overall survival of OS patients through univariate Cox regression (Figure 2A). Meanwhile, we also found that the expression of 26 ubiquitin-related genes was positively correlated with overall survival, and others were negative. Then, lasso Cox regression was applied to further screen candidate genes of ubiquitin-related gene signature. The results showed that 12 ubiquitin-related genes were regarded as good candidates (Figures 2B,C). Eventually, we screened the seven ubiquitin-related genes (UBE2L3, CORO6, DCAF8, DNAI1, FBXL5, UHRF2, and WDR53) to construct prognostic signature through multivariate Cox regression. The coefficient of each gene was shown in Figure 2D. Based on the results of PCA, we can find that all OS cases were well divided into different distribution patterns (Figure 2E).

FIGURE 2
www.frontiersin.org

FIGURE 2. Construction of a ubiquitin-related gene signature. (A) Identifying 42 ubiquitin-related genes involved with the overall survival through Univariate Cox regression analysis. (B–C) Twelve ubiquitin-related genes were screened through Lasso Cox regression analysis. (D) Coefficients of the seven ubiquitin-related genes through multivariate Cox analyses. (E) Principal component analysis (PCA) based on the signature.

Evaluation and validation of a ubiquitin-related gene signature

After constructing a ubiquitin-related gene signature, all samples were divided into high-risk and low-risk groups. The Figure 3E showed the differential expression of these seven genes between high and low-risk groups. The Kaplan–Meier survival curve suggested that there was a significantly difference in overall survival between the high- and low-risk groups in the Target set (Figure 3A) and validation (GSE21257) set (Figure 4B). Further, the ROC curves presented good prediction efficiency for OS overall survival in the Target set and validation set (Figure 3B and Figures 4C,D). And in the Target set, the respective AUC values at 1/3/5 years were 0.957, 0.890, and 0.919. Moreover, risk score distribution and survival status of OS patients were shown in (Figure 3C and Figures 4E–H), and we constructed the heatmap of seven ubiquitin-related genes by combining clinical information and risk score groups in the Target set (Figure 3D). To further explore the function of seven ubiquitin-related genes, we performed a co-expression analysis of these seven genes in the data matrix containing all genes and obtained the 226 co-expression genes with seven ubiquitin-related genes in the Supplementary Table 3 (R > 0.5, p < 0.001). Functional annotation of co-expression genes showed in Figure 3F through the Metascape. Taken together, the above results proved that the signature had a good performance in predicting prognostic efficiency.

FIGURE 3
www.frontiersin.org

FIGURE 3. Evaluating a ubiquitin-related gene signature. (A) The Kaplan–Meier survival curve in the Target set. (B) The ROC curves in the Target set. (C) Risk score distribution and survival status of OS patients in the Target set. (D) Heatmap of seven ubiquitin-related genes with combining clinical information in the Target set. (E) Differential expression of these seven genes between high and low-risk groups. (F) Functional annotation of co-expression genes through the metascape. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

FIGURE 4
www.frontiersin.org

FIGURE 4. Validation of a ubiquitin-related gene signature in two independent cohorts. (A,B) The Kaplan–Meier survival curve in two validation sets. (C,D) The ROC curves in the validation sets. (E-H) Risk score distribution and survival status of OS patients in the validation sets.

Assessment of the prognostic value of the risk score model and establishment of the nomogram

As above described, the risk score model had an advantage in prediction efficiency. Then, we compared the risk score model with other clinicopathological features. The ROC curves showed that the prediction ability of the risk score model was better than other clinicopathological features, including age (AUC = 0.427), gender (AUC = 0.421), and prognosis stage (AUC = 0.701) (Figure 5A). Moreover, multivariate Cox regression indicated that the risk score model and prognosis stage were also independent prognostic prediction factors (Figure 5B). Additionally, a nomogram with a risk score model and other clinicopathological features was constructed to better predict the survival rate of OS patients (Figure 5C). Based on the nomogram, we can obtain the survival rate of 1/3/5 years by calculating the total points of every patient. Furthermore, the 1/3/5 years calibration curve was shown in Figure 5D. Overall, data showed that the risk score model had a better advantage in prognostic prediction than other clinicopathological factors.

FIGURE 5
www.frontiersin.org

FIGURE 5. Assessment prognostic value of risk score model and the nomogram. (A) The prediction ability of the risk score model and other clinicopathological features. (B) Analyses of prognostic prediction factors through multivariate Cox analyses. (C) Nomogram for predicting 1/3/5-years survival rates of OS patients. (D) 1/3/5-year calibration curve of the nomogram.

Correlation between immune and risk score model

In order to further explore the correlation between the immune and risk score model, we carried out the analyses of immune cells and immune-related functions in the matrix of seven genes after screening and the reference dataset of CIBERSORT. As shown in Figure 6A, we described the 22 immune cells landscape between high- and low-risk groups. The results showed that the percentage of T cells CD4 memory resting, Macrophages M0, and Macrophages M2 were higher than others in OS samples. Furthermore, the fraction of T cells CD4 memory activated (p = 0.035) and NK cells resting (p = 0.027) had a significant difference in the high- and low-risk groups (Figures 6B–D). Moreover, according to the result of immune-related functions, we found that inflammation-promoting, T cell co-inhibition, and parainflammation were significantly different in the risk score groups.

FIGURE 6
www.frontiersin.org

FIGURE 6. Correlation between immune and risk score model. (A) The 22 immune cells landscape between high- and low-risk groups. (B) Comparison with the percentage of immune cells in the high- and low-risk groups. (C–D) A significant percentage of T cells CD4 memory activated and NK cells resting. (E) The significant immune-related functions in the different risk score groups (p < 0.05).

GSEA, GO, and KEGG

GSEA was applied to explore significant biological functions and pathways in the risk score groups.

As shown in Figure 7A, the top five enrichment of signaling pathways, such as cytosolic DNA sensing pathway, regulation of autophagy, RIG I like receptor signaling pathway were significantly enriched in the high-risk groups. Furthermore, according to the result of enrichment analysis, significant biological functions including B cell receptor signaling pathway, complement activation, FC epsilon receptor signaling pathway, FC epsilon receptor-mediated stimulatory signaling pathway, humoral immune response mediated circulating immune signaling pathway were significantly enriched in the low-risk groups (Figure 7B). Taken together, we can find that most of the GSEA were related to immune between high- and low-risk groups. Moreover, GO and KEGG were used to explore the biological functions and pathways of the prognosis-related ubiquitin-related genes. The top five enrichment of GO and KEGG were shown in Figures 7C,D.

FIGURE 7
www.frontiersin.org

FIGURE 7. GSEA, GO, and KEGG. (A) The GSEA showed significant enrichment in the low-risk groups. (B) The GSEA showed significant enrichment in the low-risk groups. (C) The significant GO terms are based on the prognosis-related ubiquitin-related genes. (D) The significant KEGG pathways are based on the prognosis-related ubiquitin-related genes.

Chemotherapy drug sensitivity

Standard chemotherapy for OS patients included doxorubicin, cisplatin, methotrexate, and paclitaxel (Patel et al., 1996). In order to further explore the clinical significance of the risk score model, we especially analyzed the drug sensitivity of the above chemotherapy drug and a proteasome inhibitor (MG-132) between high- and low-risk groups (Figures 8A–E). According to the results of chemotherapy drug sensitivity, half-maximal inhibitory concentration (IC50) of doxorubicin, cisplatin, methotrexate, and paclitaxel were not significant in different risk groups. Interestingly, IC50 of MG-132 was notable significance between high- and low-risk groups. Based on the result, we proposed a hypothesis that MG-132 may be a novel candidate for OS chemotherapy.

FIGURE 8
www.frontiersin.org

FIGURE 8. Chemotherapy drug sensitivity between the high- and low-risk group. (A) IC50 of Cisplatin. (B) IC50 of Methotrexate. (C) IC50 of MG-132. (D) IC50 of Paclitaxel. (E) IC50 of Doxorubicin.

Different immune microenvironment and immune-related functions in the osteosarcoma subtypes

Furthermore, in order further to explore ubiquitin-related genes, we developed novel OS subtypes. First, 42 prognosis-related ubiquitin-related genes were selected to cluster. Based on NMF algorithm, when the cophenetic was 3, we eventually obtained molecular subtypes of three cluster (Figure 9A). In addition, the cluster map was shown in Figure 9B. Meanwhile, we also made a comparison of prognosis and immune with three clusters. The Kaplan–Meier survival curve showed that cluster 2 had the worst survival outcomes than cluster 1 and cluster 3 (Figure 9C). In addition, we also analyzed the difference in the immune microenvironment (stromal score, immune score, ESTIMATE score), and the results indicated that the stromal score, immune score, ESTIMATE score of cluster 3 was higher than cluster 2 (Figures 9D–F). Additionally, differences in immune cell infiltration and immune-related functions were shown in Figures 9G,H. The results confirmed that scores of macrophages and T helper cells, cytolytic activity and type II IFN response were significant differences in the three clusters. From the above results, we can find that cluster 3 showed better survival outcomes, and we speculated the reason that cluster 3 could have stronger immunoreactivity.

FIGURE 9
www.frontiersin.org

FIGURE 9. Different immune microenvironment and immune-related functions in the osteosarcoma subtypes. (A) The value of the nonnegative matrix factorization (NMF). (B) The NMF clustering-based value of best cophenetic. (C) Kaplan–Meier curves of OS in the three clusters. (D–F) The ESTIMATE score, stromal score, and immune score of OS subtypes. (G) The percentage of immune cells in the three clusters. (H) The significant immune-related functions in the three clusters.

Validating expression of the seven genes in normal and osteosarcoma tissues

Based on better prediction efficiency, we further obtained differential expression of these seven genes in six normal and six osteosarcoma sample tissues through the qRT-PCR. The results showed that the mRNA expression level of four ubiquitin-related genes (CORO6, UBE2L3, FBXL5, DNAI1) was significantly increased in osteosarcoma tissues as compared with adjacent normal tissues. Moreover, the expression of DCAF8 was markedly decreased in osteosarcoma tissues than in adjacent normal tissues, and the expression of UHRF2 and WDR53 was not significant between osteosarcoma and normal tissues (Figure 10A). Based on the differential expression of these seven genes in clinical samples, it is necessary to explore the biological functions and mechanisms for OS treatment.

FIGURE 10
www.frontiersin.org

FIGURE 10. Validating expression of the seven genes in normal and osteosarcoma tissues. (A) The mRNA expression level of four ubiquitin-related genes (UBE2L3, CORO6, DCAF8, DNAI1, FBXL5, UHRF2, WDR53) in six normal and six osteosarcoma samples tissues. *p < 0.05, **p < 0.01, ***p < 0.001.

Discussion

In recent years, a growing number of evidence has suggested that PTM was closely related to cell survival, cell cycle, differentiation, and innate and adaptive immunity for multiple cancer (Wu et al., 2019; Chen et al., 2020; Yang et al., 2020). As one of the most important PTM, ubiquitination has also attracted much attention in OS (Chen et al., 2014; Ying et al., 2016; Zhao et al., 2019). With the gradually increasing studies, regulating ubiquitination has been regarded as a novel therapeutic for OS (Jiang et al., 2020; Jiang et al., 2021; Mullard et al., 2021). Based on the above evidence, we further explored the significance of ubiquitin-related genes in OS prognosis through bioinformatic analysis.

In this study, we developed a ubiquitin-related gene signature including seven genes (UBE2L3, CORO6, DCAF8, DNAI1, FBXL5, UHRF2, WDR53) through univariate Cox regression, Lasso Cox regression, and multivariate Cox regression. Meanwhile, we confirmed that the gene signature had a good performance in predicting prognosis for OS patients. For the seven ubiquitin-related genes, there were a number of studies about their functions in human disease. For example, Tao et al. (2020) reported that overexpression of UBE2L3 could promote tumor progression through decreasing protein stability of GSK3β. Moreover, the small molecule PSSM0332 has been proved to inhibit inflammatory response by decreasing the DCAF8-medicated ubiquitination in the myocardial dysfunction model (Peng et al., 2020). Djakow et al. (2012) suggested that mutations DNAH5 and DNAI1 had an important significance in the diagnosis of primary ciliary dyskinesia. In addition, the research confirmed that FBXL5 could regulate HIF-1α activity by promoting the ubiquitination of CITED2 (Machado-Oliveira et al., 2015). Taken together, although there were a series of studies on seven ubiquitin-related genes, the reports associated with OS were poor. Hence, it will be necessary to explore their molecular biological functions in OS initiation and progression. In addition, many kinds of gene signature for OS were gradually increasing, such as immune-related gene signature (Xiao et al., 2020), hypoxia-associated prognostic signature (Fu et al., 2021), ferroptosis-related gene signature (Lei et al., 2021), and pyroptosis-related signature (Zhang et al., 2021). Among the above prognostic signature, our signature showed better predictive efficacy than others. Notably, there was still a need to verify the efficacy of the signature in the prospective study.

Currently, many studies suggested that ubiquitination had an important effect on the immune response. For example, increasing evidence demonstrated that targeting ubiquitination of programmed death 1 (PD-1)/programmed death-ligand 1 (PD-L1) could be a promising therapeutic approach for cancer immunotherapy (Hu et al., 2021). Hence, we explored the immune cells and immune-related functions between the high- and low-risk groups. Moreover, inflammation-promoting, T cell co-inhibition, and parainflammation exhibited a significant difference in the risk score model. Then, the results of GSEA showed some enrichment with immune. For instance, B cell receptor signaling pathway, complement activation, FC epsilon receptor signaling pathway, FC epsilon receptor-mediated stimulatory signaling pathway. In summary, the ubiquitin-related gene signature provided a good prediction for immune response.

All the time, doxorubicin, cisplatin, methotrexate, and paclitaxel are among the most chemotherapeutic agents for OS treatment (Patel et al., 1996; Bielack et al., 2015). However, a growing number of evidence indicated that resistance to standard chemotherapy is still a common problem for OS treatment. In our study, the drug sensitivity of the standard chemotherapy had no significance in the different risk models. Interestingly, IC50 of MG-132 showed a significant difference between high- and low-risk groups. Actually, emerging research have suggested that MG-132 played an important role in cancer chemotherapy (Pajonk et al., 2005; Li et al., 2007; Zhang et al., 2008). Moreover, Saini et al. (2021) have reported that MG-132 could be regarded as an effective therapeutic approach for OS treatment by inducing regulation of autophagy and protein homeostasis. MG-132 was also proved to modulate OS differentiation by regulating the expression of IRS-1 (Contaldo et al., 2014). Although some studies have verified that MG-132 has a pharmacological effect on OS treatment, it is necessary to explore the function of MG-132 in the clinical cohort.

In recent years, molecular subtypes have been proved to have instructive roles in preclinical and clinical therapy in many cancer types (Collisson et al., 2019). And in our study, we developed the three ubiquitination patterns through NMF algorithm. The results confirmed that cluster 3 showed better survival outcomes. Moreover, cluster 3 had the higher score in the stromal score, immune score, ESTIMATE score. Meanwhile, based on the relationship between ubiquitination and immune response, we compared the immune cell infiltration and immune-related functions in the three clusters. Overall, the above data indicated the probable conclusion that a better prognosis in cluster 3 was associated with stronger immunoreactivity.

Although we constructed a novel ubiquitin-related gene signature for predicting OS prognosis, there are still some limitations in this study. First, the ubiquitin-related gene signature only showed the relationship to OS prognosis, so it is necessary to explore their biological function for OS tumorigenesis. Second, due to the limitations of clinical factors, such as tumor grade and stage, we further demonstrated that the risk score model is an independent prognostic factor correlated with OS prognosis in large prospective clinical trials. Moreover, the clinical significance of OS subtypes should also be verified in the future.

In conclusion, we developed a novel ubiquitin-related gene signature which showed better predictive prognostic ability for OS and provided additional information on chemotherapy and immunotherapy. The OS molecular subtypes would also give a useful guide for individualized therapy.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm.nih.gov/geo/, GSE21257 https://www.ncbi.nlm.nih.gov/geo/, GSE39055.

Author contributions

ZH-h and MY-c worked on designing this study. NW, GC-y, ZW-m, and LZ-y were involved in data analysis and drafting the manuscript. SY-q, ZS-b, and ZK participated in collecting the data and revising English grammar. All authors consistently approved the final manuscript.

Funding

This study was supported by the Natural Science Foundation of China (31960175); Gansu Youth Science and Technology Fund (20JR10RA752); Gansu Key R&D Program (21YF5FA125); Lanzhou talent innovation and entrepreneurship project (2021-RC-99); and Cuiying Scientific and Technological Innovation Program (CY2021-MS-B17) of Lanzhou University Second Hospital.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2022.904448/full#supplementary-material

References

Anwar, M. A., El-Baba, C., Elnaggar, M. H., Elkholy, Y. O., Mottawea, M., Johar, D., et al. (2020). Novel therapeutic strategies for spinal osteosarcomas. Semin. Cancer Biol. 64, 83–92. doi:10.1016/j.semcancer.2019.05.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Bandos, A. I., Rockette, H. E., Song, T., and Gur, D. (2009). Area under the free-response ROC curve (FROC) and a related summary index. Biometrics 65, 247–256. doi:10.1111/j.1541-0420.2008.01049.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bielack, S. S., Smeland, S., Whelan, J. S., Marina, N., Jovic, G., Hook, J. M., et al. (2015). Methotrexate, doxorubicin, and cisplatin (MAP) plus maintenance pegylated interferon alfa-2b versus MAP alone in patients with resectable high-grade osteosarcoma and good histologic response to preoperative MAP: First results of the EURAMOS-1 good response randomized controlled trial. J. Clin. Oncol. 33, 2279–2287. doi:10.1200/JCO.2014.60.0734

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanche, P., Dartigues, J. F., and Jacqmin-Gadda, H. (2013). Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat. Med. 32, 5381–5397. doi:10.1002/sim.5958

PubMed Abstract | CrossRef Full Text | Google Scholar

Botter, S. M., Neri, D., and Fuchs, B. (2014). Recent advances in osteosarcoma. Curr. Opin. Pharmacol. 16, 15–23. doi:10.1016/j.coph.2014.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., Liu, S., and Tao, Y. (2020). Regulating tumor suppressor genes: Post-translational modifications. Signal Transduct. Target. Ther. 5, 90. doi:10.1038/s41392-020-0196-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Shen, J., Li, X., Wang, X., Long, M., Lin, F., et al. (2014). Rlim, an E3 ubiquitin ligase, influences the stability of Stathmin protein in human osteosarcoma cells. Cell. Signal. 26, 1532–1538. doi:10.1016/j.cellsig.2014.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Collisson, E. A., Bailey, P., Chang, D. K., and Biankin, A. V. (2019). Molecular subtypes of pancreatic cancer. Nat. Rev. Gastroenterol. Hepatol. 16, 207–220. doi:10.1038/s41575-019-0109-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Contaldo, C., Myers, T. J., Zucchini, C., Manara, M. C., Chiodoni, C., Colombo, M. P., et al. (2014). Expression levels of insulin receptor substrate-1 modulate the osteoblastic differentiation of mesenchymal stem cells and osteosarcoma cells. Growth factors (Chur, Switz. 32, 41–52. doi:10.3109/08977194.2013.870168

CrossRef Full Text | Google Scholar

Djakow, J., Svobodová, T., Hrach, K., Uhlík, J., Cinek, O., and Pohunek, P. (2012). Effectiveness of sequencing selected exons of DNAH5 and DNAI1 in diagnosis of primary ciliary dyskinesia. Pediatr. Pulmonol. 47, 864–875. doi:10.1002/ppul.22520

PubMed Abstract | CrossRef Full Text | Google Scholar

Engebretsen, S., and Bohlin, J. (2019). Statistical predictions with glmnet. Clin. Epigenetics 11, 123. doi:10.1186/s13148-019-0730-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, Y., Bao, Q., Liu, Z., He, G., Wen, J., Liu, Q., et al. (2021). Development and validation of a hypoxia-associated prognostic signature related to osteosarcoma metastasis and immune infiltration. Front. Cell Dev. Biol. 9, 633607. doi:10.3389/fcell.2021.633607

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS one 9, e107468. doi:10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Gómez-Rubio, V. (2017). ggplot2 - elegant graphics for data analysis. 2nd Edition Foundation for Open Access Statistics. Journal of Statal Software, Book Reviews 77, 1–3.

Google Scholar

Harrison, D. J., Geller, D. S., Gill, J. D., Lewis, V. O., and Gorlick, R. (2018). Current and future therapeutic approaches for osteosarcoma. Expert Rev. Anticancer Ther. 18, 39–50. doi:10.1080/14737140.2018.1413939

PubMed Abstract | CrossRef Full Text | Google Scholar

Hershko, A., and Ciechanover, A. (1998). The ubiquitin system. Annu. Rev. Biochem. 67, 425–479. doi:10.1146/annurev.biochem.67.1.425

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, X., Wang, J., Chu, M., Liu, Y., Wang, Z. W., and Zhu, X. (2021). Emerging role of ubiquitination in the regulation of PD-1/PD-L1 in cancer immunotherapy. Mol. Ther. 29, 908–919. doi:10.1016/j.ymthe.2020.12.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, L., Liao, J., Liu, J., Wei, Q., and Wang, Y. (2021). Geranylgeranylacetone promotes human osteosarcoma cell apoptosis by inducing the degradation of PRMT1 through the E3 ubiquitin ligase CHIP. J. Cell. Mol. Med. 25, 7961–7972. doi:10.1111/jcmm.16725

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, W., Cai, X., Xu, T., Liu, K., Yang, D., Fan, L., et al. (2020). Tripartite motif-containing 46 promotes viability and inhibits apoptosis of osteosarcoma cells by activating NF-B signaling through ubiquitination of PPAR. Oncol. Res. 28, 409–421. doi:10.3727/096504020X15868639303417

PubMed Abstract | CrossRef Full Text | Google Scholar

Kansara, M., Teng, M. W., Smyth, M. J., and Thomas, D. M. (2014). Translational biology of osteosarcoma. Nat. Rev. Cancer 14, 722–735. doi:10.1038/nrc3838

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, D. D., and Seung, H. S. (1999). Learning the parts of objects by non-negative matrix factorization. Nature 401, 788–791. doi:10.1038/44565

PubMed Abstract | CrossRef Full Text | Google Scholar

Lei, T., Qian, H., Lei, P., and Hu, Y. (2021). Ferroptosis-related gene signature associates with immunity and predicts prognosis accurately in patients with osteosarcoma. Cancer Sci. 112, 4785–4798. doi:10.1111/cas.15131

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W., Zhang, X., and Olumi, A. F. (2007). MG-132 sensitizes TRAIL-resistant prostate cancer cells by activating c-Fos/c-Jun heterodimers and repressing c-FLIP(L). Cancer Res. 67, 2247–2255. doi:10.1158/0008-5472.CAN-06-3793

PubMed Abstract | CrossRef Full Text | Google Scholar

Lorent, M., Giral, M., and Foucher, Y. (2014). Net time-dependent ROC curves: A solution for evaluating the accuracy of a marker to predict disease-related mortality. Stat. Med. 33, 2379–2389. doi:10.1002/sim.6079

PubMed Abstract | CrossRef Full Text | Google Scholar

Machado-Oliveira, G., Guerreiro, E., Matias, A. C., Facucho-Oliveira, J., Pacheco-Leyva, I., and Bragança, J. (2015). FBXL5 modulates HIF-1α transcriptional activity by degradation of CITED2. Arch. Biochem. Biophys. 576, 61–72. doi:10.1016/j.abb.2015.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Mangiola, S., Doyle, M. A., and Papenfuss, A. T. (2021). Interfacing Seurat with the R tidy universe. Bioinformatics 37 (2021), 4100–4107. doi:10.1093/bioinformatics/btab404

CrossRef Full Text | Google Scholar

Meazza, C., and Scanagatta, P. (2016). Metastatic osteosarcoma: A challenging multidisciplinary treatment. Expert Rev. Anticancer Ther. 16, 543–556. doi:10.1586/14737140.2016.1168697

PubMed Abstract | CrossRef Full Text | Google Scholar

Mullard, M., Lavaud, M., Regnier, L., Tesfaye, R., Ory, B., Rédini, F., et al. (2021). Ubiquitin-specific proteases as therapeutic targets in paediatric primary bone tumours? Biochem. Pharmacol. 194, 114797. doi:10.1016/j.bcp.2021.114797

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakayama, K. I., and Nakayama, K. (2006). Ubiquitin ligases: Cell-cycle control and cancer. Nat. Rev. Cancer 6, 369–381. doi:10.1038/nrc1881

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Steen, C. B., Liu, C. L., Gentles, A. J., Chaudhuri, A. A., Scherer, F., et al. (2019). Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat. Biotechnol. 37, 773–782. doi:10.1038/s41587-019-0114-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Pajonk, F., van Ophoven, A., Weissenberger, C., and McBride, W. H. (2005). The proteasome inhibitor MG-132 sensitizes PC-3 prostate cancer cells to ionizing radiation by a DNA-PK-independent mechanism. BMC cancer 5, 76. doi:10.1186/1471-2407-5-76

PubMed Abstract | CrossRef Full Text | Google Scholar

Patel, S. R., Papadopoulos, N. E., Plager, C., Linke, K. A., Moseley, S. H., Spirindonidis, C. H., et al. (1996). Phase II study of paclitaxel in patients with previously treated osteosarcoma and its variants. Cancer 78, 741–744. doi:10.1002/(SICI)1097-0142(19960815)78:4<741:AID-CNCR8>3.0.CO;2-H

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, Q., Xu, H., Xiao, M., and Wang, L. (2020). The small molecule PSSM0332 disassociates the CRL4A(DCAF8) E3 ligase complex to decrease the ubiquitination of NcoR1 and inhibit the inflammatory response in a mouse sepsis-induced myocardial dysfunction model. Int. J. Biol. Sci. 16, 2974–2988. doi:10.7150/ijbs.50186

PubMed Abstract | CrossRef Full Text | Google Scholar

Saini, H., Sharma, H., Mukherjee, S., Chowdhury, S., and Chowdhury, R. (2021). Verteporfin disrupts multiple steps of autophagy and regulates p53 to sensitize osteosarcoma cells. Cancer Cell Int. 21, 52. doi:10.1186/s12935-020-01720-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, N. N., Zhang, Z. Z., Ren, J. H., Zhang, J., Zhou, Y. J., Wai Wong, V. K., et al. (2020). Overexpression of ubiquitin-conjugating enzyme E2 L3 in hepatocellular carcinoma potentiates apoptosis evasion by inhibiting the GSK3β/p65 pathway. Cancer Lett. 481, 1–14. doi:10.1016/j.canlet.2020.03.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Z., Huang, R., and Yuan, L. (2019). Crosstalk of intracellular post-translational modifications in cancer. Arch. Biochem. Biophys. 676, 108138. doi:10.1016/j.abb.2019.108138

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, B., Liu, L., Li, A., Xiang, C., Wang, P., Li, H., et al. (2020). Identification and verification of immune-related gene prognostic signature based on ssGSEA for osteosarcoma. Front. Oncol. 10, 607622. doi:10.3389/fonc.2020.607622

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Zhou, W., and Lin, H. (2020). Posttranslational modifications of smurfs: Emerging regulation in cancer. Front. Oncol. 10, 610663. doi:10.3389/fonc.2020.610663

PubMed Abstract | CrossRef Full Text | Google Scholar

Ying, M., Zhang, L., Zhou, Q., Shao, X., Cao, J., Zhang, N., et al. (2016). The E3 ubiquitin protein ligase MDM2 dictates all-trans retinoic acid-induced osteoblastic differentiation of osteosarcoma cells by modulating the degradation of RARα. Oncogene 35, 4358–4367. doi:10.1038/onc.2015.503

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics a J. Integr. Biol. 16, 284–287. doi:10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

Zhang, N. H., Song, L. B., Wu, X. J., Li, R. P., Zeng, M. S., Zhu, X. F., et al. (2008). Proteasome inhibitor MG-132 modifies coxsackie and adenovirus receptor expression in colon cancer cell line lovo. Cell cycleGeorget. Tex.) 7, 925–933. doi:10.4161/cc.7.7.5621

CrossRef Full Text | Google Scholar

Zhang, W., Zhuang, Y., Zhang, Y., Yang, X., Zhang, H., Wang, G., et al. (2017). Uev1A facilitates osteosarcoma differentiation by promoting Smurf1-mediated Smad1 ubiquitination and degradation. Cell Death Dis. 8, e2974. doi:10.1038/cddis.2017.366

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., He, R., Lei, X., Mao, L., Jiang, P., Ni, C., et al. (2021). A novel pyroptosis-related signature for predicting prognosis and indicating immune microenvironment features in osteosarcoma. Front. Genet. 12, 780780. doi:10.3389/fgene.2021.780780

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, W., Zhang, D., Qin, P., Zhang, J., Cui, X., Gao, J., et al. (2019). Long non-coding RNA EPIC1 inhibits viability and invasion of osteosarcoma cells by promoting MEF2D ubiquitylation. Int. J. Biol. Macromol. 128, 566–573. doi:10.1016/j.ijbiomac.2019.01.156

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, C., Zhang, Z., Zhu, X., Qian, G., Zhou, Y., Sun, Y., et al. (2020). N6-Methyladenosine modification of the TRIM7 positively regulates tumorigenesis and chemoresistance in osteosarcoma through ubiquitination of BRMS1. EBioMedicine 59, 102955. doi:10.1016/j.ebiom.2020.102955

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: osteosarcoma, ubiquitin, ubiquitination, immune, molecular subtypes

Citation: Wei N, Chao-yang G, Wen-ming Z, Ze-yuan L, Yong-qiang S, Shun-bai Z, Kai Z, Yan-chao M and Hai-hong Z (2022) A ubiquitin-related gene signature for predicting prognosis and constructing molecular subtypes in osteosarcoma. Front. Pharmacol. 13:904448. doi: 10.3389/fphar.2022.904448

Received: 25 March 2022; Accepted: 19 July 2022;
Published: 17 August 2022.

Edited by:

Chien-Shan Cheng, Shanghai Jiao Tong University, China

Reviewed by:

Yun Liu, Guangxi Medical University, China
Zhuochao Liu, Shanghai Jiao Tong University, China
Hong Sun, Affiliated Hospital of Guizhou Medical University, China

Copyright © 2022 Wei, Chao-yang, Wen-ming, Ze-yuan, Yong-qiang, Shun-bai, Kai, Yan-chao and Hai-hong. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Ma Yan-chao, bWF5Y2xkZXlAMTYzLmNvbQ==; Zhang Hai-hong, ZXJ5X3poYW5naGhlckBsenUuZWR1LmNu

These authors have contributed equally to this work

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.