Skip to main content
Download PDF

System analysis based on Anoikis-related genes identifies MAPK1 as a novel therapy target for osteosarcoma with neoadjuvant chemotherapy

BMC Musculoskeletal Disorders volume 25, Article number: 437 (2024) Cite this article

Abstract

Background

Osteosarcoma (OS) is the most common bone malignant tumor in children, and its prognosis is often poor. Anoikis is a unique mode of cell death.However, the effects of Anoikis in OS remain unexplored.

Method

Differential analysis of Anoikis-related genes was performed based on the metastatic and non-metastatic groups. Then LASSO logistic regression and SVM-RFE algorithms were applied to screen out the characteristic genes. Later, Univariate and multivariate Cox regression was conducted to identify prognostic genes and further develop the Anoikis-based risk score. In addition, correlation analysis was performed to analyze the relationship between tumor microenvironment, drug sensitivity, and prognostic models.

Results

We established novel Anoikis-related subgroups and developed a prognostic model based on three Anoikis-related genes (MAPK1, MYC, and EDIL3). The survival and ROC analysis results showed that the prognostic model was reliable. Besides, the results of single-cell sequencing analysis suggested that the three prognostic genes were closely related to immune cell infiltration. Subsequently, aberrant expression of two prognostic genes was identified in osteosarcoma cells. Nilotinib can promote the apoptosis of osteosarcoma cells and down-regulate the expression of MAPK1.

Conclusions

We developed a novel Anoikis-related risk score model, which can assist clinicians in evaluating the prognosis of osteosarcoma patients in clinical practice. Analysis of the tumor immune microenvironment and chemotherapeutic drug sensitivity can provide necessary insights into subsequent mechanisms. MAPK1 may be a valuable therapeutic target for neoadjuvant chemotherapy in osteosarcoma.

Peer Review reports

Introduction

Osteosarcoma is a threatening bone malignancy that is common in adolescents and children [1]. Current treatments include surgical resection and pre-and postoperative chemotherapy, with five-year survival rates of up to 70% for non-metastatic osteosarcoma patients [2]. Osteosarcoma, a type of bone cancer, is known for its resistance to conventional chemotherapy treatments, which can have a significant impact on the patient’s prognosis [3]. Our aim is to identify the genes that hinder patients’ response to treatment and find suitable chemotherapy drugs. This would help us discover more effective targeted therapies that can overcome treatment resistance and serve as prognostic indicators for patients with osteosarcoma [4].

Anoikis is a type of programmed cell death that occurs when cells separate from the correct extracellular matrix, thereby disrupting integrin ligation. It is a critical mechanism that prevents dysplastic cells from growing abnormally or attaching to inappropriate substrates [5]. Anoikis resistance was found to be an important mechanism of cancer growth, invasion, and metastasis [6]. Numerous studies have shown that Anoikis plays a significant role in the occurrence and development of cancer. CPT1A promotes colorectal cancer cell metastasis by inhibiting Anoikis [7]. In addition, the CamKK2-AMPK signaling pathway leads to the metastasis and invasion of LKB1-deficient lung cancer [8]. However, the correlation between Anoikis and the prognosis of osteosarcoma still needs to be elucidated.

In the present study, we constructed a new model of risk score based on the Anoikis-related genes in osteosarcoma. We further explored the predictive value of these genes and investigated the association between Anoikis and the tumor immune microenvironment. We also focused on the analysis of related genes and chemotherapeutic susceptibility. Our findings will provide a novel perspective for predicting individualized survival and better treatment of osteosarcoma patients.

Materials and methods

Data collection

The gene sets information, and corresponding clinical data of 84 osteosarcoma patients(63 nonmetastasis and 21 metastasis) were obtained from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database. Clinical information included gender, age, diagnosis, metastasis, survival time, status, and follow-up. We downloaded another dataset of 53 osteosarcoma patients from the Gene Expression Omnibus (GEO) database (GSE21257). We retrieved the 639 Anoikis-related gene lists from the Genecards database (https://www.genecards.org/) and Harmonizome portals (https://maayanlab.cloud/Harmonizome/) [9] (Shown in Table S1 and Table S2). The expression data were normalized to fragment per kilobase million (FPKM) values. The calculation of FPKM for gene i uses the following formula [10]:

$$FPK{M_i} = {{{q_i}} \over {{{{l_i}} \over {{{10}^3}}}*{{\sum\nolimits_j {{q_j}} } \over {{{10}^6}}}}} = {{{q_i}} \over {{l_i}*\sum\nolimits_j {{q_j}} }}*{10^9}$$

qi are raw read or fragment counts, \(l_{i}\) is feature length, \(\sum\limits_j {{q_j}}\) and corresponds to the total number of mapped reads or fragments.

Identification of differentially expressed Anoikis-related genes

We applied the “limma” R package to find different expressions of Anoikis-related genes. Then we used the “heatmap” R package to accomplish the heatmaps of differentially expressed genes (DEGs).

Machine learning algorithm

LASSO logistic regression and support vector machine recursive feature elimination (SVM-RFE) algorithms were used to screen out characteristic genes. LASSO logistic regression was performed based on the “glmnet” R package [11]. SVM analysis, a machine learning method that depends on the “e1071” R package, can find the best variables by deleting the feature vectors generated by SVM [12].

Development and validation of the Anoikis-related genes (ARGS) prognostic model

Univariable Cox regression analysis was used to assess the prognostic value of each Anoikis-related DEG in the TARGET cohort. Those genes with p < 0.05 were chosen for further investigation, and multivariate Cox regression analysis was applied to shrink the potential genes and build the prognostic prediction model. The formula calculating for risk score was “Riskscore = Gene A∗Coef A + Gene B ∗ Coef B + … + Gene N ∗ Coef N” [13]. High- and low-risk score groups were divided according to the median Anoikis-related risk score of the training cohort. And the Kaplan-Meier analysis was conducted to compare survival possibility and overall survival time between the high- and low-risk groups. The area under the curve (AUC) was calculated to assess the sensitivity and specificity of the risk score system. Then, we used the univariate and multivariate Cox regression analysis to validate the independence of our prognosis model. Finally, we underwent external validation for the Anoikis-related risk score system in the GEO cohort.

Development of a nomogram

We applied the univariate and multivariable Cox regression analysis (“survival” R package) to assess the risk score combined with clinical information (age, gender, and metastatic status) in the TARGET cohort. Based on the results of Cox regression analysis, a nomogram was presented to predict the prognosis of osteosarcoma patients. Besides, the calibration curve was shown to evaluate the nomogram’s performance.

Gene set enrichment analysis (GSEA)

Gene set enrichment analyses (GSEA) software was used to identify the enriched pathways between the two cluster groups based on the curated genesets (go.v7.4.symbols.gmt and kegg.v7.4.symbols.gmt) [14].

Immune Cell Infiltration and Immune function analysis

We used the “GSVA” R package to quantify immune cell infiltration and immune function in the TARGET cohort [15]. Subsequently, we analyzed the correlation between the two clusters and different enriched pathways.

Drug susceptibility analysis

The “OncoPredict” R package was applied to predict in vivo drug responses in cancer patients [16]. OncoPredict fits the gene expression profile of tissues to the half-maximal inhibitory concentration (IC50) of the cancer cell lines to drugs downloaded from Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org/) and the gene expression profile of cancer lines from the Broad Institute Cancer Cell Line Encyclopedia (CCLE; https://portals.broadinstitute.org/ccle_legacy/home). 198 drugs were calculated in total, and the sensitivity of the drugs was analyzed using unpaired t-tests. p < 0.05 was considered statistically significant.

In addition, we downloaded the NCI-60 human cancer cell lines from the CellMiner database (https://discover.nci.nih.gov/cellminer) [17]. Pearson correlation analysis was applied to assess the relationship between MAPK1 and chemosensitivity.

Single-cell sequencing analysis

The Tumor Immune Single Cell Hub (TISCH, http://tisch.comp-genomics.org) database is a comprehensive website that can realize the visualization of the tumor immune microenvironment [18]. Gene expression data were downloaded from the GEO database (GSE162454). All data were uniformly processed with a standardized cell-type annotation and differential expression analysis.

Cell lines and cultures

One human osteoblast cell line (hFOB1.19) and two human osteosarcoma cell lines (U-2OS and MG-63) were obtained from the National Collection of Authenticated Cell Cultures (Shanghai, China). Dulbecco’s modified Eagle’s medium (DMEM, Gibco) contains 1% penicillin/streptomycin (Thermo Fisher Scientific, United States) and 10% fetal bovine serum (FBS, Gibco). We cultured human osteoblast cells in the medium at 34℃ with 5% CO2, and the osteosarcoma cells at 37℃ with 5% CO2.

Cell viability assay

Following the manufacturer’s protocol, U-2OS and MG-63 cell survival was assessed via the CCK-8 kit. About 5 × 10³ cells were extracted from U-2OS and MG-63 cell suspensions, respectively, and incubated at 37℃ and 5% CO2 for 24 h in each well of the 96-well plate. The cells were then treated with different concentrations of Nilotinib (0, 10, 20, and 30 µM). Cells were washed using phosphate-buffered saline (PBS) after incubation. 100 µl of DMEM containing 10 µl CCK-8 solution was added to each well, then the mixture was incubated for 2–4 h. The absorbance of the wells was measured using a microplate reader at 450 nm.

Western blotting

After cell treatments, whole-cell proteins from hFOB, U-2OS, and MG-63 cells were extracted using commercial kits (Beyotime) according to the manufacturer’s instructions. The protein quantification was determined using the BCA protein assay kit. Then 20 ng of protein from each group was resolved via sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS PAGE) and transferred to a PVDF membrane. Membranes were blocked with 5% non-fat milk. Then the blots were cut prior to incubated with primary antibodies against MYC (1:1000), MAPK1 (1:1000), GAPDH (1:1000), Bcl-2 (1:1000), Bax (1:1000), and CASP3 (1:1000) overnight at 4℃. Membranes were washed and further incubated with the respective secondary antibodies. Electrochemiluminescence plus reagent (Invitrogen) was used to detect the bands. Blots were imaged and quantified using Image Lab 3.0 software.

TUNEL staining

To measure U-2OS apoptosis under various 24 h treatments, TUNEL labeling was used. After being fixed for 15 min at room temperature in 4% paraformaldehyde (PFA), U-2OS were rinsed with PBS and permeabilized for 3 min on ice using 0.1% Triton X-100 buffer. Following observation under a confocal microscope, apoptotic U-2OS were stained with the TUNEL staining kit reagents, and the nuclei were counterstained with DAPI for 10 min. Apoptotic U-2OS were then counted and analyzed.

RNA extraction and real-time PCR analysis (RT-PCR)

Following the manufacturer’s instructions, total RNA was isolated from the hFOB, U-2OS, and MG-63 using TRIzol (Invitrogen), and then reverse-transcribed into cDNA (MBI Fermentas, Germany), and the RT-PCR reaction was performed in the RT-PCR system (Bio-Rad Laboratories, CA, USA), according to the operating instruction. The expression levels of relative genes were calculated via a comparative quantification method (2 − ΔΔCt formula) and were normalized to internal control, GAPDH. MYC and MAPK1 primers were designed with the NCBI Primer-Blast Tool, and they are presented below:

MYC

  1. a.

    (F) 5’ GGCTCCTGGCAAAAGGTCA 3’,

  2. b.

    (R) 5’ CTGCGTAGTTGTGCTGATGT 3’.

MAPK1

  1. c.

    (F) 5’ TACACCAACCTCTCGTACATCG 3’,

  2. d.

    (R) 5’ CATGTCTGAAGCGCAGTAAGATT 3’.

Statistical analysis

All statistical analysis was conducted by using R software (Version:3.6.1) and GraphPad Prism (Version:7.00). Comparisons between two independent groups were applied using a two-tailed, unpaired t-test. Two-way analysis of variance (ANOVA) with Tukey’s multiple comparisons test was applied to analyze differences among three or more groups when the data were normally distributed. What’s more, nonparametric Mann-Whitney U tests were used for groups if the data were not normally distributed. P value < 0.05 was considered statistically significant.

Results

Identification of characteristic genes

The 639 Anoikis-related gene expression levels were compared in the TARGET database from metastatic and non-metastatic tissues. We identified 28 DEGs (all P < 0.05). The RNA levels of these Anoikis-related genes are shown in Fig. 1A (red: high expression level; blue: low expression level). LASSO logistic regression and SVM-RFE algorithms were applied to screen out the characteristic genes related to metastasis (Fig. 1B, C). Finally, we ended up with 12 genes which were selected between LASSO and SVM- RFE algorithms for further research., including MAPK1, ERBB2, CALR, PIK3CG, MYC, SIRT3, PIN1, CASP6, FBLIM1, PIP5K1C, DNMT1, and CXCL14 (Fig. 1D).

Fig. 1
figure 1

Identification of anoikis-related gene signature. (A) A heatmap (blue: low expression level; red: high expression level) of genes between the metaststic and the non-metaststic tumor tissues. (B) LASSO coefficient profiles. (C) SVM- RFE algorithm. (D) Twelve characteristic genes were selected between LASSO and SVM- RFE algorithms

Full size image

Cancer classification based on the DEGs

We further investigated the association between the expression level of 12 Anoikis-related DEGs and osteosarcoma. Hence, we performed a consensus clustering analysis with 84 osteosarcoma patients. The intragroup correlations decreased as the clustering variable (k) increased from 2 to 10. when k = 2, the 84 osteosarcoma patients could be well divided into two clusters (Fig. 2A). We compared survival rates among the two clusters and found significant differences among these clusters’ patients (p < 0.01, Fig. 2B). Then the principal components analysis (PCA) results indicated that patients separated into A and B groups had significantly different discrimination (Fig. 2C). The relationship between the expression of these genes and the clinical characteristics, including gender, age (< 15 or > 15 years), metastasis status (metastatic or non-metastatic), and primary tumor site (Arm/hand Leg/Foot and other), is displayed in a heatmap (Fig. 2D).

Fig. 2
figure 2

Tumor classification based on DEGs. (A) 84 osteosarcoma patients were grouped into two clusters via the consensus clustering matrix (k = 2). (B) Kaplan–Meier overall survival curves for the two clusters. (C) PCA plot. (D) A heatmap (blue: low expression level; red: high expression level)

Full size image

Immune status and tumor microenvironment (TME)

To study the TME differences between the two clusters, we further used ssGSEA to evaluate the enrichment scores of 23 types of immune cells. Several immune cell infiltrates were significantly enriched in group A, such as Macrophages, CD8 T cells, activated NK cells, and others (Fig. 3A). Gene Set Variation Analysis (GSVA) revealed that the following pathways were significantly activated in group A: JAK/STAT signaling pathway, chemokine signaling pathway, and others. Interestingly, Arginine and Proline Metabolism pathways were much more active in group B (Fig. 3B, C).

Fig. 3
figure 3

Immune status and tumor microenvironment. (A) The ssGSEA analysis for immune cells between two clusters. (B) The GSVA (Gene Set Variation Analysis) for KEGG pathways between two clusters. (C) GSEA (Gene set enrichment analysis) in B cluster (low survival group)

Full size image

Development of a prognostic model in the TARGET cohort

A total of 84 osteosarcoma specimens were matched to corresponding patients with complete survival information. Univariate Cox regression analysis was used to initially screen genes associated with survival. The 4 genes (MAPK1, PIK3CG, MYC, EDIL3) that met the criteria of P < 0.05 were further analyzed. Then multivariate Cox regression analysis was performed, and a three-gene signature was constructed. The risk score was calculated as follows: risk score = (-0.884* MAPK1 exp.) + (0.656* MYC exp.) + (-0.772* EDIL3 exp.). According to the median risk score, we divided 84 osteosarcoma patients into high- and low-risk subgroups (Fig. 4A). High-risk patients had more deaths and shorter survival times than low-risk patients (Fig. 4B). The Kaplan-Meier curve showed that overall survival time and possibility were significantly lower in the high-risk group (Fig. 4C, P < 0.001). The value of the area under the curve was 0.948 for 1-year, 0.788 for 3-year, and 0.783 for 5-year survival prediction (Fig. 4D).

Fig. 4
figure 4

Construction of risk signature in the TARGET cohort. (A) Distribution of patients based on the risk score. (B) The survival status for each patient (left side of the dotted line: low-risk population; right side of the dotted line: high-risk population). (C) Kaplan–Meier curves for the overall survival of patients between the high- and low-risk groups. (D) ROC curves

Full size image

External validation of risk score in a GEO cohort

53 osteosarcoma patients from a GEO cohort (GSE21257) were extracted as the external validation set. Based on the median risk score of the TARGET model, 30 patients were regarded as the high-risk group, while the other 23 were at low risk (Fig. 5A). The risk scores, survival time, and survival status of patients are shown in Fig. 5B. High-risk patients had higher mortality and shorter overall survival time. Besides, the Kaplan-Meier curve results also showed a lower survival possibility in the high-risk group (Fig. 5C, P = 0.024). AUC values of external validation also showed an optimistic prediction, and the AUC was 0.801 for 1 year, 0.787 for 3 years, and 0.744 for 5 years (Fig. 5D).

Fig. 5
figure 5

Validation of the risk model in the GSE21257. (A) Distribution of patients in the GSE21257 based on the median risk score of the TARGET cohort. (B) The survival status for each patient. (C) Kaplan–Meier curves. (D) Time-dependent ROC curves for osteosarcoma

Full size image

Independent prognostic value of the risk model

Univariate and multivariable Cox regression models were conducted to assess independent prognostic factors of the gene-based risk score and clinical characteristics. In the TARGET cohort, univariate Cox analysis showed the metastatic status, primary tumor site, and risk score were significantly associated with prognosis (Fig. 6A).

Fig. 6
figure 6

Construction of the predictive model. (A) Independence detection of the constructed risk prediction model. (B) A prognostic model to predict overall survival in the TARGET cohort. (C) Calibration curves of the OS nomogram model

Full size image

Nomogram

According to the prognostic model and clinical factors (age, gender, and metastatic status), we developed a risk estimation nomogram in the TARGET cohort (Fig. 6B). 1-, 3-, and 5-year calibration curves showed that the nomogram consistently predicted the survival rate (Fig. 6C).

Immune microenvironment analysis

The immune infiltration of 22 immune cells was investigated by the CIBERSORT algorithm (Fig. 7A). As shown in Fig. 7B, there were apparent correlations between various immune cells in the prognostic model. Furthermore, we validated the correlation of MAPK1, PIK3CG, and EDIL3 expression and immune cell infiltration in the datasets of GSE162454 from the TISCH database (Fig. 7C). The results showed MAPK1, PIK3CG, and EDIL3 play a valuable role in the fibroblast cell(Fig. 7D). In osteosarcoma, immunotherapy often faces hurdles posed by cancer-associated fibroblasts (CAFs) that secrete dense extracellular matrix components and cytokines. Directly removing CAFs may prove ineffective and even promote tumor metastasis [19]. Therefore, we speculate that the metastasis and deterioration of osteosarcoma can be inhibited by influencing these three targets.

Fig. 7
figure 7

Immune microenvironment. (A) The immune infiltration of 22 immune cells between low- and high-risk group. (B) An immunocyte related heatmap. (C, D) single-cell cluster in TISCH database

Full size image

OncoPredict for drug susceptibility analysis

To explore suitable drugs for patients with high-risk scores, we transformed the gene expression of osteosarcoma tissues in the TARGET group into a drug sensitivity matrix using the OncoPredict algorithm (Fig. 8). Osteosarcoma tissues from high-risk group patients were more sensitive to 6 drugs than those osteosarcoma tissues from low-risk group patients, including those of Dihydrorotenone (mitochondrial inhibitor), MG-132 (autophagy activator), Sabutoclax(targeting drug, Bcl-2 inhibitor), Telomerase Inhibitor IX(Telomerase Inhibitor), Vorinostat(HDAC1 inhibitor), and VX-11e ERK inhibitor). These drugs may help in the treatment of osteosarcoma patients with high-risk scores.

Fig. 8
figure 8

OncoPredict for drug susceptibility analysis

Full size image

Verification the expression of two predictive genes

To verify the expression levels of MAPK1 and MYC in osteosarcoma, we performed Western blotting and RT-PCR analysis on the osteoblast cell line hFOB1.19 and two osteosarcoma cell lines (U-2OS and MG-63). The results showed that the expression levels of MYC and MAPK1 were significantly upregulated in two osteosarcoma cell lines (U-2OS and MG-63) compared with the osteoblast cell line hFOB1.19 (Fig. 9A, B).

Fig. 9
figure 9

The expression levels of two genes between osteosarcoma cell lines and osteoblasts. (A) The qRT-PCR results of MYC and MAPK1. (B) Western blotting results of MYC and MAPK1 expression

Full size image

Nilotinib can decrease osteosarcoma cell viability and down-regulate MAPK1 expression

We further investigated the sensitivity of MAPK1 to chemotherapeutic agents. The results showed that MAPK1 was sensitive to Nilotinib (Fig. 10A, p = 0.035). Nilotinib is an anti-vascular targeted agent that promotes apoptosis of several sarcoma cell lines, thereby inhibiting the metastasis of osteosarcoma [20]. However, there are few reports in the literature. To evaluate the therapeutic effect of Nilotinib on osteosarcoma, two types of osteosarcoma cell lines (U-2OS and MG-63) were treated with different doses of Nilotinib (10, 20, and 30 µM). Then the U-2OS and MG-63 cell survival rate was assessed via the CCK-8 kit. We found that Nilotinib has a dose-dependent cytotoxic effect on osteosarcoma cell lines (Fig. 10B, C). Subsequently, we used the western blotting analysis to detect the expression level of MAPK1 in the U-2OS and MG-63 osteosarcoma cell lines in response to Nilotinib (30µM). The results showed that the expression of MAPK1 was down-regulated in both osteosarcoma cell lines after Nilotinib treatment (Fig. 10D, E, F).

Fig. 10
figure 10

Nilotinib can reduce cell viability and MAPK1 expression in osteosarcoma cells. (A) Scatter plot of relationship between MAPK1 expression and drug sensitivity. (B) Evaluation of MG-63 osteosarcoma cell viability using CCK-8 assay after exposure to nilotinib for 24 h. (C) Evaluation of U-2OS osteosarcoma cell viability using CCK-8 assay after exposure to nilotinib for 24 h. (D, E, F) The expression level of MAPK1 protein in osteosarcoma cells. GAPDH serves as an internal standard. The gels have been run under the same experimental conditions. All experiments were repeated in triplicates (n = 3). The obtained data are represented as mean ± SE. Significance: *p-value < 0.05, **p-value < 0.01, ***p-value < 0.001

Full size image

Nilotinib promotes apoptosis in two osteosarcoma cells

The results of Western blotting analysis demonstrated that the expression levels of Bax and CASP3 were enhanced. In contrast, the expression of Bcl-2 decreased after applying Nilotinib to treat MG-63 osteosarcoma cells. What’s more, the effect of the high-concentration group (30µM) on apoptosis-related proteins in MG-63 osteosarcoma cells was more evident than the low-concentration group (10µM) (Fig. 11A, B). In addition, we applied TUNEL to evaluate the treatment of Nilotinib for U-2OS osteosarcoma cells. As expected, the number of TUNEL-positive cells was upregulated with increased Nilotinib concentration (Fig. 11C, D).

Fig. 11
figure 11

(A, B) The protein expression levels of Bcl-2, Bax, and CASP3 in osteosarcoma cells. (C, D) TUNEL staining was used to detect osteosarcoma cell apoptosis (bar: 50 μm; nuclei: blue; positive cells: green). All experiments were repeated in triplicates (n = 3). The obtained data are represented as mean ± SE. Significance: *p-value < 0.05, **p-value < 0.01, ***p-value < 0.001

Full size image

Discussion

Osteosarcoma, one of the common malignant bone tumors, occurs predominantly in the long bone epiphysis of children and adolescents and often develops metastasis [21]. Effective therapeutic strategies, including surgery, radiotherapy, and chemotherapy, are considered against osteosarcoma [22, 23]. However, the prognosis of osteosarcoma patients is still poor, and the 5-year survival rate is low [24]. Osteosarcoma patients with the same clinical risk factors may significantly differ in prognosis and treatment [25]. Therefore, it is of great significance for the early diagnosis, targeted therapy, and prognosis analysis of osteosarcoma to deeply understand the molecular pathological mechanism and screen key biomarkers related to the occurrence and development of osteosarcoma. At present, many risk scoring systems and prognosis predictions have been developed in clinical application and improvement of patient prognosis management.For example, The four pseudo-genetic markers developed by Xiaoqiang Zhang et al. apply to patients of different sex, age, and metastatic status. These four pseudogenes are involved in the regulation of malignant phenotype, immunity, and DNA/RNA editing, and have a good predictive effect on the treatment of osteosarcoma patients [26]. The result of this study proves that machine learning has a good prediction effect. In addition, the use of various algorithms will take your essay to the next level. Lai et al. used a variety of different algorithms to construct a new regulation network of gene-lncRNA-pathway-immunocyte [27]. miRNA is the main target of function. Therefore we decided to combine the fields of machine learning and biology to develop an entirely new gene predictive model.

Programmed cell death is regulated by various genes and plays an important role in the growth and development of organisms. It is also essential for maintaining tissue and organ homeostasis and is involved in a variety of pathological processes. In addition to apoptosis, iron death, necrotization and cell pyrodeath also contribute to the occurrence and development of cancer [28].Anoikis is an important defense of the organism. Once normal epithelial cells lose contact with the extracellular matrix (ECM), they rapidly undergo apoptosis [29]. However, a common feature of tumor development and growth is the ability of transformed cells to survive under independent growth conditions [30]. This resistance to Anoikis has been shown to be associated with loss of intracellular environmental homeostasis, cancer growth, and metastasis, and this acquired ability is known as lost-nest apoptosis resistance [31]. Cancer cells with Anoikis resistance can spread to distant tissues or organs through the peripheral circulatory system and cause cancer metastasis [32]. In this context, Anoikis resistance is a natural molecular prerequisite for the spread of invasive cancer metastases [33]. Studying the molecular mechanism controlling Anoikis resistance will help search for effective therapies for malignant tumors.

In our study, we comprehensively assessed the Anoikis-related genes in osteosarcoma. We found 28 Anoikis-related genes were differentially expressed between metastasis and non-metastasis groups. Then we used Lasso and SVM-RFE algorithms to screen out feature genes. The two clusters generated by consensus clustering analysis based on the feature genes showed significant differences in survival probability. Next, we constructed a 3-gene risk signature by univariable Cox and multivariate Cox regression analyses. Further, we evaluated the prognostic value of these Anoikis-related gene regulators in training and validation cohorts.

Among the three prognosis gene signatures, MAPK1 was a significant target in osteosarcoma treatment [34]. MAPK1/3 kinase can attenuate Mitophagy and promotes breast cancer bone metastasis [35]. In addition, Mitophagy which plays an important role in carcinogenesis and tumor progression, occurs through an alternative autophagy pathway, requiring the MAPK1 and MAPK14 signaling pathways [36, 37]. Previous studies have shown that ezrin is required for metastasis of osteosarcoma in mouse models, and high expression levels are often associated with adverse outcomes in dogs and patients with osteosarcoma [38]. Ezrin’s ability to attach cell membranes to the actin cytoskeleton allows the cell to interact directly with its microenvironment, thereby facilitating signal transduction through growth factor receptors and adhesion molecules. Furthermore, in mouse models, Ezrin-mediated metastasis survival was found to be partially dependent on MAPK signaling [38, 39]. Therefore, we speculate that MAPK1 may be an effective target to influence the resistance of osteosarcoma patients to anoikis by altering ezrin expression, thereby affecting the prognosis of osteosarcoma patients. MYC is a transcription factor that dimerizes with MAX to bind DNA and regulate gene expression [40]. It has been known that MYC promotes cell growth and proliferation in normal cells, but it also contributes to the genesis of many human cancers [41,42,43]. MYC could mediate cancer cell energy metabolism and may be a new anticancer therapy [44]. Moreover, previous research also proposes that therapies targeting the MYC pathway will be vital to reversing cancerous growth and restoring antitumor immune responses in patients with MYC-driven cancers [45]. Inhibition of PML could lead to a remarkable growth arrest associated with a decrease in MYC kinase levels [46]. EDIL3 was identified as a novel regulator of epithelial-mesenchymal transition (EMT), contributing to angiogenesis, metastasis, and recurrence of hepatocellular carcinoma [47]. EDIL3 is also a strong and highly accurate diagnostic marker for breast cancer [48]. Moreover, in accordance with previous studies, we found that tumor-derived EDIL3 was involved in tumor-associated bone loss [49]. We believe these three genes may be essential to osteosarcoma’s occurrence, development, and prognosis.

Targeting the tumor immune and bone microenvironment could open up new therapeutic opportunities for patients [50]. According to the immune cell Infiltration results, this new scoring system is closely related to the tumor immune microenvironment. We then used single-cell sequencing analysis to verify our results, showing that the three genes characteristic of the novel scoring model are closely related to immune cells.

Currently, there are many factors that contribute to cancer chemotherapy resistance. For example, the regulation of circular RNA on downstream targets [51, 52]. We want to find new therapeutic targets and drugs. Subsequently, we evaluated the association between relevant prognostic genes and chemotherapeutic drug sensitivity to apply our findings to clinical treatment. From the results, we found that MAPK1 is an important therapeutic target sensitive to various chemotherapy drugs.

We found that the expression level of MAPK1 was significantly decreased in osteosarcoma cells treated with Nilotinib. What’s more, the expression level of MAPK1 decreased in a dose-dependent manner. In addition, through a series of studies, we found that Nilotinib could significantly increase the apoptosis of osteosarcoma cells by down-regulating MAPK1 expression. We suggest that Nilotinib may provide a better therapeutic effect for osteosarcoma patients with elevated MAPK1 expression. MAPK1 may be a key anoikis-related target for the treatment of osteosarcoma.

However, there are some limitations to our study. Firstly, we constructed and validated our risk score model based on public databases, and the sample size was not rich enough. Therefore, our model needs to be further validated based on our clinical data in the future. Moreover, the mechanism of these three predictive genes could be more precise. We will conduct comprehensive functional experiments and multi-omics analysis in our future research.

Conclusion

In our study, we looked forward to exploring rational prognostic predictors for osteosarcoma patients with metastasis. To predict prognosis, we developed and validated an Anoikis-based risk score system for osteosarcoma patients with metastasis. The AUC value also showed this Anoikis-related risk score system had a good prediction performance. Our study provides a novel gene signature for predicting the prognosis of osteosarcoma patients with metastasis. It opens an avenue for future studies of the relationships between Anoikis-related genes and immunity in these patients. Finally, MAPK1 may be a vital biotherapeutic target.

Data availability

The data that support the findings of this study are openly available in the Therapeutically Applicable Research to Generate Effective Treatments (TARGET; https://ocg.cancer.gov/programs/target) and GEO (https://www.ncbi.nlm.nih.gov/geo/) databases at GSE21257 and GSE162454.

References

  1. Kansara M, Teng MW, Smyth MJ, Thomas DM. Translational biology of osteosarcoma. Nat Rev Cancer. 2014;14(11):722–35.

    Article  CAS  PubMed  Google Scholar 

  2. Isakoff MS, Bielack SS, Meltzer P, Gorlick R. Osteosarcoma: current treatment and a collaborative pathway to Success. J Clin Oncol. 2015;33(27):3029–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Kempf-Bielack B, Bielack SS, Jurgens H, Branscheid D, Berdel WE, Exner GU, Gobel U, Helmke K, Jundt G, Kabisch H, et al. Osteosarcoma relapse after combined modality therapy: an analysis of unselected patients in the Cooperative Osteosarcoma Study Group (COSS). J Clin Oncol. 2005;23(3):559–68.

    Article  PubMed  Google Scholar 

  4. Lilienthal I, Herold N. Targeting Molecular Mechanisms Underlying Treatment Efficacy and Resistance in Osteosarcoma: A Review of Current and Future Strategies. Int J Mol Sci 2020, 21(18).

  5. Taddei ML, Giannoni E, Fiaschi T, Chiarugi P. Anoikis: an emerging hallmark in health and diseases. J Pathol. 2012;226(2):380–93.

    Article  CAS  PubMed  Google Scholar 

  6. Kakavandi E, Shahbahrami R, Goudarzi H, Eslami G, Faghihloo E. Anoikis resistance and oncoviruses. J Cell Biochem. 2018;119(3):2484–91.

    Article  CAS  PubMed  Google Scholar 

  7. Wang YN, Zeng ZL, Lu J, Wang Y, Liu ZX, He MM, Zhao Q, Wang ZX, Li T, Lu YX, et al. CPT1A-mediated fatty acid oxidation promotes colorectal cancer cell metastasis by inhibiting anoikis. Oncogene. 2018;37(46):6025–40.

    Article  CAS  PubMed  Google Scholar 

  8. Jin L, Chun J, Pan C, Kumar A, Zhang G, Ha Y, Li D, Alesi GN, Kang Y, Zhou L, et al. The PLAG1-GDH1 Axis promotes Anoikis Resistance and Tumor Metastasis through CamKK2-AMPK Signaling in LKB1-Deficient Lung Cancer. Mol Cell. 2018;69(1):87–e9987.

    Article  CAS  PubMed  Google Scholar 

  9. Rouillard AD, Gundersen GW, Fernandez NF, Wang Z, Monteiro CD, McDermott MG, Ma’ayan A. The harmonizome: a collection of processed datasets gathered to serve and mine knowledge about genes and proteins. Database (Oxford) 2016, 2016.

  10. Zhao YD, Li MC, Konaté MM, Chen L, Das B, Karlovich C, Williams PM, Evrard YA, Doroshow JH, McShane LM. TPM, FPKM, or normalized counts? A comparative study of quantification measures for the analysis of RNA-seq data from the NCI patient-derived models repository. J Transl Med 2021, 19(1).

  11. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized Linear models via Coordinate Descent. J Stat Softw. 2010;33(1):1–22.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Huang ML, Hung YH, Lee WM, Li RK, Jiang BR. SVM-RFE based feature selection and Taguchi parameters optimization for multiclass SVM classifier. ScientificWorldJournal. 2014;2014:795624.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Peng C, Li L, Luo G, Tan S, Xia R, Zeng L. Integrated analysis of the M2 macrophage-related signature associated with prognosis in ovarian cancer. Front Oncol. 2022;12:986885.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Powers RK, Goodspeed A, Pielke-Lombardo H, Tan AC, Costello JC. GSEA-InContext: identifying novel and common patterns in expression experiments. Bioinformatics. 2018;34(13):i555–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform 2021, 22(6).

  17. Reinhold WC, Sunshine M, Liu H, Varma S, Kohn KW, Morris J, Doroshow J, Pommier Y. CellMiner: a web-based suite of genomic and pharmacologic tools to explore transcript and drug patterns in the NCI-60 cell line set. Cancer Res. 2012;72(14):3499–511.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Sun D, Wang J, Han Y, Dong X, Ge J, Zheng R, Shi X, Wang B, Li Z, Ren P, et al. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment. Nucleic Acids Res. 2021;49(D1):D1420–30.

    Article  CAS  PubMed  Google Scholar 

  19. Wang H, Chen Y, Wei R, Zhang JL, Zhu JH, Wang WB, Wang ZF, Wupur Z, Li YJ, Meng H. Synergistic chemoimmunotherapy augmentation via sequential nanocomposite hydrogel-mediated reprogramming of Cancer-Associated fibroblasts in Osteosarcoma. Adv Mater 2023.

  20. Alemany R, Moura DS, Redondo A, Martinez-Trufero J, Calabuig S, Saus C, Obrador-Hevia A, Ramos R, Villar VH, Valverde C, et al. Nilotinib as Coadjuvant Treatment with Doxorubicin in patients with sarcomas: a phase I trial of the Spanish Group for Research on Sarcoma. Clin Cancer Res. 2018;24(21):5239–49.

    Article  CAS  PubMed  Google Scholar 

  21. Anderson ME. Update on Survival in Osteosarcoma. Orthop Clin North Am. 2016;47(1):283–92.

    Article  PubMed  Google Scholar 

  22. Miwa S, Shirai T, Yamamoto N, Hayashi K, Takeuchi A, Igarashi K, Tsuchiya H. Current and emerging targets in Immunotherapy for Osteosarcoma. J Oncol. 2019;2019:7035045.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Xie L, Ji T, Guo W. Anti-angiogenesis target therapy for advanced osteosarcoma (review). Oncol Rep. 2017;38(2):625–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Harrison DJ, Geller DS, Gill JD, Lewis VO, Gorlick R. Current and future therapeutic approaches for osteosarcoma. Expert Rev Anticancer Ther. 2018;18(1):39–50.

    Article  CAS  PubMed  Google Scholar 

  25. Ingley KM, Maleddu A, Grange FL, Gerrand C, Bleyer A, Yasmin E, Whelan J, Strauss SJ. Current approaches to management of bone sarcoma in adolescent and young adult patients. Pediatr Blood Cancer. 2022;69(2):e29442.

    Article  PubMed  Google Scholar 

  26. Liu F, Xing L, Zhang XQ, Zhang XQ. A four-pseudogene classifier identified by machine learning serves as a novel prognostic marker for survival of Osteosarcoma. Genes-Basel 2019, 10(6).

  27. Zhang XQ, Ren LH, Yan XY, Shan Y, Liu L, Zhou J, Kuang QY, Li MQ, Long H, Lai WL. Identification of immune-related lncRNAs in periodontitis reveals regulation network of gene-lncRNA-pathway-immunocyte. Int Immunopharmacol 2020, 84.

  28. Zhu XD, Li SL. Ferroptosis, Necroptosis, and Pyroptosis in Gastrointestinal Cancers: The Chief Culprits of Tumor Progression and Drug Resistance. Adv Sci 2023, 10(26).

  29. Gilmore AP. Anoikis. Cell Death Differ. 2005;12(Suppl 2):1473–7.

    Article  CAS  PubMed  Google Scholar 

  30. Mullard A. Addressing cancer’s grand challenges. Nat Rev Drug Discov. 2020;19(12):825–6.

    Article  CAS  PubMed  Google Scholar 

  31. Fanfone D, Wu Z, Mammi J, Berthenet K, Neves D, Weber K, Halaburkova A, Virard F, Bunel F, Jamard C et al. Confined migration promotes cancer metastasis through resistance to anoikis and increased invasiveness. Elife 2022, 11.

  32. Simpson CD, Anyiwe K, Schimmer AD. Anoikis resistance and tumor metastasis. Cancer Lett. 2008;272(2):177–85.

    Article  CAS  PubMed  Google Scholar 

  33. Luo M, Huang M, Yang N, Zhu Y, Huang P, Xu Z, Wang W, Cai L. Impairment of rigidity sensing caused by mutant TP53 gain of function in osteosarcoma. Bone Res. 2023;11(1):28.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Liu JF, Chen PC, Chang TM, Hou CH. Monocyte chemoattractant Protein-1 promotes cancer cell migration via c-Raf/MAPK/AP-1 pathway and MMP-9 production in osteosarcoma. J Exp Clin Cancer Res. 2020;39(1):254.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Deng R, Zhang HL, Huang JH, Cai RZ, Wang Y, Chen YH, Hu BX, Ye ZP, Li ZL, Mai J, et al. MAPK1/3 kinase-dependent ULK1 degradation attenuates mitophagy and promotes breast cancer bone metastasis. Autophagy. 2021;17(10):3011–29.

    Article  CAS  PubMed  Google Scholar 

  36. Panigrahi DP, Praharaj PP, Bhol CS, Mahapatra KK, Patra S, Behera BP, Mishra SR, Bhutia SK. The emerging, multifaceted role of mitophagy in cancer and cancer therapeutics. Semin Cancer Biol. 2020;66:45–58.

    Article  CAS  PubMed  Google Scholar 

  37. Hirota Y, Yamashita S, Kurihara Y, Jin X, Aihara M, Saigusa T, Kang D, Kanki T. Mitophagy is primarily due to alternative autophagy and requires the MAPK1 and MAPK14 signaling pathways. Autophagy. 2015;11(2):332–43.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Khanna C, Wan XL, Bose S, Cassaday R, Olomu O, Mendoza A, Yeung C, Gorlick R, Hewitt SM, Helman LJ. The membrane-cytoskeleton linker ezrin is necessary for osteosarcoma metastasis. Nat Med. 2004;10(2):182–6.

    Article  CAS  PubMed  Google Scholar 

  39. Wan XL, Mendoza A, Khanna C, Helman LJ. Rapamycin inhibits ezrin-mediated metastatic behavior in a murine model of osteosarcoma. Cancer Res. 2005;65(6):2406–11.

    Article  CAS  PubMed  Google Scholar 

  40. Conacci-Sorrell M, McFerrin L, Eisenman RN. An overview of MYC and its interactome. Cold Spring Harb Perspect Med. 2014;4(1):a014357.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Stine ZE, Walton ZE, Altman BJ, Hsieh AL, Dang CV. MYC, metabolism, and Cancer. Cancer Discov. 2015;5(10):1024–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Dang CV. MYC on the path to cancer. Cell. 2012;149(1):22–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Zhong Y, Yang L, Xiong F, He Y, Tang Y, Shi L, Fan S, Li Z, Zhang S, Gong Z, et al. Long non-coding RNA AFAP1-AS1 accelerates lung cancer cells migration and invasion by interacting with SNIP1 to upregulate c-Myc. Signal Transduct Target Ther. 2021;6(1):240.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Dang CV, Le A, Gao P. MYC-induced cancer cell energy metabolism and therapeutic opportunities. Clin Cancer Res. 2009;15(21):6479–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Dhanasekaran R, Deutzmann A, Mahauad-Fernandez WD, Hansen AS, Gouw AM, Felsher DW. The MYC oncogene - the grand orchestrator of cancer growth and immune evasion. Nat Rev Clin Oncol. 2022;19(1):23–36.

    Article  CAS  PubMed  Google Scholar 

  46. Arreal L, Piva M, Fernandez S, Revandkar A, Schaub-Clerigue A, Villanueva J, Zabala-Letona A, Pujana M, Astobiza I, Cortazar AR, et al. Targeting PML in triple negative breast cancer elicits growth suppression and senescence. Cell Death Differ. 2020;27(4):1186–99.

    Article  CAS  PubMed  Google Scholar 

  47. Xia H, Chen J, Shi M, Gao H, Sekar K, Seshachalam VP, Ooi LL, Hui KM. EDIL3 is a novel regulator of epithelial-mesenchymal transition controlling early recurrence of hepatocellular carcinoma. J Hepatol. 2015;63(4):863–73.

    Article  CAS  PubMed  Google Scholar 

  48. Wei YX, Han JH, Shen HM, Wang YY, Qi M, Wang L, Li J. Highly sensitive fluorescent detection of EDIL3 overexpressed exosomes for the diagnosis of triple-negative breast cancer. Nanotechnology 2022, 33(42).

  49. Kun Z, Xin G, Tao W, Chenglong Z, Dongsheng W, Liang T, Tielong L, Jianru X. Tumor derived EDIL3 modulates the expansion and osteoclastogenesis of myeloid derived suppressor cells in murine breast cancer model. J Bone Oncol. 2019;16:100238.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Li SL, Kang Y, Zeng Y. Targeting tumor and bone microenvironment: novel therapeutic opportunities for castration-resistant prostate cancer patients with bone metastasis. Bba-Rev Cancer 2024, 1879(1).

  51. Guo X, Gao CY, Yang DH, Li SL. Exosomal circular RNAs: a chief culprit in cancer chemotherapy resistance. Drug Resist Update 2023, 67.

  52. Li SL, Liu F, Zheng K, Wang W, Qiu ED, Pei Y, Wang S, Zhang JM, Zhang XJ. CircDOCK1 promotes the tumorigenesis and cisplatin resistance of osteogenic sarcoma via the miR-339-3p/IGF1R axis. Mol Cancer 2021, 20(1).

Download references

Acknowledgements

We like to acknowledge the TARGET and the GEO (GSE21257) network for providing data.

Funding

Not applicable.

Author information

Author notes

  1. Zhouwei Wu, Jiapei Yu and Tao Han contributed equally to this work.

Authors and Affiliations

  1. Department of Orthopedics, the Second Affiliated Hospital, Yuying Children’s Hospital of Wenzhou Medical University, Wenzhou, 325000, China

    Zhouwei Wu, Jiapei Yu, Yiting Tu, Fang Su, Shi Li & Yixing Huang

  2. Department of Orthopedics, the Shaoxing People’s Hospital (Shaoxing Hospital of Zhejiang University), Shaoxing, 312000, China

    Tao Han

  3. Key Laboratory of Orthopedics of Zhejiang Province, Wenzhou, 325000, China

    Zhouwei Wu, Jiapei Yu, Tao Han, Yiting Tu, Shi Li & Yixing Huang

  4. Department of Orthopaedics, The Second Affiliated Hospital, Yuying Children’s Hospital of Wenzhou Medical University, 109 West Xueyuan Road, Wenzhou, 325027, Zhejiang Province, China

    Shi Li

Authors
  1. Zhouwei Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Jiapei Yu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Tao Han
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Yiting Tu
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Fang Su
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Shi Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Yixing Huang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

HYX and LS conceived this idea. All authors designed, performed, and analyzed the experiments. WZW, YJP and HT conducted the data analysis and wrote the manuscript. All authors reviewed and approved the final manuscript.

Corresponding authors

Correspondence to Shi Li or Yixing Huang.

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.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Supplementary Material 3

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

Wu, Z., Yu, J., Han, T. et al. System analysis based on Anoikis-related genes identifies MAPK1 as a novel therapy target for osteosarcoma with neoadjuvant chemotherapy. BMC Musculoskelet Disord 25, 437 (2024). https://doi.org/10.1186/s12891-024-07547-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12891-024-07547-2

Keywords

BMC Musculoskeletal Disorders

ISSN: 1471-2474

Contact us