3.1. Explore differentially expressed MRGs in OS
The "limma" R package was utilized to pinpoint DEGs between osteosarcoma (OS) and normal skeletal muscle tissue.The selection criteria were p < 0.05 and |FoldChange| > 1.5 (|log2 FoldChange| > 0.585). The analysis revealed 12,578 key genes exhibited differential expression (Fig. 2A,B). Next, these genes were allowed to intersect with 103 mitophagy-related genes obtained from the KEGG website to get 95 overlapping genes (Fig. 2C). GO enrichment analyses emphasized the involvement of these genes in processes like mitochondrion disassembly, autophagy of mitochondrion, organelle disassembly, mitochondrial outer membrane, organelle outer membrane, ubiquitin protein ligase binding, ubiquitin-like protein ligase binding, and GDP binding. Additionally, the KEGG enrichment analysis demonstrated significant enrichment of these genes in processes related to mitophagy, autophagy, and neurodegeneration pathways, providing further evidence for the functional relevance of the selected genes (Fig. 2D-H).
3.2. Construction of prognosis prediction model
Clinical data on osteosarcoma patients was collected to further identify prognostic genes. Univariable Cox regression analysis was conducted, and 10 key genes were selected (Fig. 3A). Subsequently, LASSO regression analysis was performed to further shortlist the selected genes to 7 for model construction (Fig. 3B,C). The bar graph presented the coef values of model genes (Fig. 3D), while the chromosomal ring map showed the location of genes on the chromosomes (Fig. 3E). Patients were divided into high- and low-risk groups based on the median risk score. The survival status plot revealed that patients' survival time progressively decreased, and their mortality rate increased as the risk score increased (Fig. 3F,G). A heatmap visualized gene expression levels in OS patient samples (Fig. 3H). The KM survival curve demonstrated that patients in the high-risk group had a poor prognosis compared to those in the low-risk group (Fig. 3I). The AUC values for the ROC curves at 1, 3, and 5 years were 0.925, 0.936, and 0.948, respectively (Fig. 3J).
3.3. Validation of the prognostic model
Validation across multiple datasets (TARGET testing, total, and GSE21257) confirmed that as risk scores increased, patient survival time decreased, and mortality rates rose (Fig. 4A, B). Patients in the high-risk group had a worse prognosis than those in the low-risk group, as shown by the KM survival curve (Fig. 4C). Furthermore, the ROC curve of all groups indicated that the model exhibited robust predictive performance and stability for prognosis in OS patients (Fig. 4D).
3.4. Construction of nomogram and subgroup analysis
A nomogram integrating risk scores, age, gender, and metastasis was developed to enhance survival prediction accuracy (Fig. 5A). Furthermore, the calibration curves of the nomogram demonstrated strong predictive accuracy and performance validation, confirming the nomogram's robustness for prediction and validation (Fig. 5B). Based on the clinicopathological characteristics, the patients were placed into six groups: patients with ≥ 12 years of age, patients with < 12 years of age, male patients, female patients, non-metastatic group, and metastatic group. Among these six subgroups, the overall survival (OS) of patients in the high-risk group was lower than that of the low-risk group (p < 0.05; Fig. 5C-H). This suggests that the model possesses a remarkable capacity to differentiate between various patient outcomes.
3.5. GO, KEGG, and GSEA enrichment analysis
DEGs between high- and low-risk groups underwent GO, KEGG, and GSEA analyses. The GO enrichment analysis showed that the DEGs primarily regulated Ca2+ concentration, Ca2+ concentration homeostasis, signal transduction and ion channels (Fig. 6A,B). The GSEA analysis revealed that genes in the high-risk group were mainly enriched in chemical carcinogenesis, citrate cycle (TCA cycle), and metabolism of xenobiotics by cytochrome P450. In contrast, genes in the low-risk group were primarily enriched in Allograft rejection, chemokine signaling pathway, and transcriptional misregulation in cancer (Fig. 6C,D).
3.6. Exploration of immune cell infiltration in the TME
TME, the environment surrounding tumor cells, encompasses various components. These included blood vessels, immune cells, fibroblasts, bone marrow-derived inflammatory cells, signaling molecules, and the extracellular matrix surrounding tumor cells[16–18]. Using the Cibersort method, we found macrophages, particularly M0 and M2, to be predominant in osteosarcoma samples. Notably, the high-risk group exhibited a higher proportion of M2 macrophages, reflecting their immunosuppressive and tumor-promoting roles (Fig. 7A). The Xcell algorithm further revealed a reduction in activated dendritic cells (aDC), CD8+ T cells, cytotoxic cells, natural killer (NK) cells, T cells, T follicular helper (Tfh), T helper (Th) 1 cells, and T regulatory (Treg) cells in the high-risk group compared to the low-risk group (Fig. 7B). Correlation analysis among immune cells showed the association between 21 types of immune cells. Moreover, a high expression correlation was observed between T cells, aDC, and cytotoxic cells, indicating potential relevant associations among these cells (Fig. 7C). Afterwards, the correlation between immune cells, risk score, and the expression level of model genes was explored. An increase in the MRAS expression led to a higher infiltration of multiple immune cells, confirming that MRAS acts as a protective factor in the progression of osteosarcoma. The increasing presence of immune cells that target tumors might lead to more efficient surveillance and elimination of tumor cells, thereby impeding the growth and spread of tumors (Fig. 7D). Furthermore, in the high-risk group, it was observed that lower scores in antigen-presenting cell (APC) co-inhibition, CCR, checkpoint, inflammation-promoting, major histocompatibility complex (MHC) class I, T cell co-stimulation, T cell co-inhibition, type I interferon (IFN) response, and parainflammation pathways compared to the low-risk group. This reflects potential functional defects or suppressive states of the immune system in high-risk patients across multiple levels (Fig. 7E). The TME score analysis also revealed decreased ESTIMATE, stromal, and immune scores, along with increased tumor purity in the high-risk group (Fig. 7F-I). These findings suggest insufficient infiltration of stromal cells and immune cells in the TME, potentially leading to faster tumor progression and a higher risk of metastasis.
3.7. Validation of Model Genes in OS cell line and OS Clinical Samples
RT-qPCR analysis further validated the expression of model genes in three osteosarcoma cell lines (HOS, MG63, and U2OS). Normal osteoblast cells (hFOB1.19)served as the control group. Significantly elevated expression was observed for SP1, E2F1, and RPS27A in OS cell lines, while RAB5C, AMBRA1, MRAS, and SMURF1 expression decreased in osteosarcoma cell lines(Fig. 8A-G). A total of 25 osteosarcoma tumors and corresponding tissue samples were obtained from the Second Affiliated Hospital of Nanchang University between January 2021 and December 2023. The expression levels in human tissue samples were further validated, and the results were consistent with the validation of OS cell lines. This suggests that the developed risk prediction model is effective and reliable (Fig. 8H-N).
3.8. Overexpression of MRAS inhibits proliferation, migration, and invasion and promotes apoptosis in osteosarcoma cells.
The lentiviral transfection was utilized to upregulate MRAS expression (Figs. 9A and 10A) and assess its effects on cell proliferation, migration, invasion, and apoptosis. CCK-8 and Edu assays demonstrated that MRAS overexpression reduced HOS and MG63 cell proliferation (Figs. 9B and 10B). Wound healing assay also indicated that MRAS overexpression significantly suppressed the migration ability of HOS cells (Figs. 9C and 10C). The Transwell experiment proved that overexpressing MRAS strongly inhibited the migration and invasion of HOS cells (Figs. 9D,E and 10D,E). Flow cytometry analysis revealed that overexpression of MRAS inhibited apoptosis in HOS cells (Figs. 9F and 10F). In conclusion, overexpression of the MRAS gene can suppress cell proliferation, migration, and invasion and simultaneously promote apoptosis.
3.9. Pan-cancer analysis
Through comprehensive pan-cancer analysis, it was found that MRAS expression is elevated in cholangiocarcinoma (CHOL), kidney renal papillary cell carcinoma (KIRP), and liver hepatocellular carcinoma (LIHC). Conversely, its expression is lower in 13 different tumor types, including bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), colon adenocarcinoma (COAD), glioblastoma (GBM), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), pheochromocytoma and paraganglioma (PCPG), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), and uterine corpus endometrial carcinoma (UCEC) (Fig. 11A). Kaplan-Meier survival analysis was subsequently employed to assess the relationship between MRAS and patient outcomes, specifically disease-specific survival (DSS), disease-free interval (DFI), and overall survival (OS). In terms of DSS, MRAS served as a prognostic protective factor in LUAD, KIRP, and skin cutaneous melanoma (SKCM), but as a prognostic risk factor in UCEC (Fig. 11B). For DFI, MRAS was beneficial in thymoma (THYM) and SKCM, but detrimental in adrenocortical carcinoma (ACC) and UCEC (Fig. 11C). Regarding OS, MRAS was protective in LUAD and SKCM, but a risk factor in UCEC (Fig. 11D). Tumor mutational burden (TMB), a biomarker indicating responsiveness to immune checkpoint inhibitors, especially those targeting the PD-1/PD-L1 pathway, was also investigated in relation to MRAS expression. Notably, MRAS expression correlated with TMB in 10 tumor types. Positive correlations were observed in THYM and KIRP, while negative correlations were found in UCEC, PRAD, LIHC, stomach adenocarcinoma (STAD), pancreatic adenocarcinoma (PAAD), brain lower grade glioma (LGG), head and neck squamous cell carcinoma (HNSC), and esophageal carcinoma (ESCA) (Fig. 11E). Additionally, the relationship between MRAS expression and microsatellite instability (MSI), a hypermutation phenotype linked to PD-1 blockade efficacy, was examined. Among 11 tumor types, only BRCA and COAD showed a positive correlation between MRAS expression and MSI. All other tumor types, including testicular germ cell tumors (TGCT), acute myeloid leukemia (LAML), UCEC, STAD, LUSC, KICH, HNSC, ESCA, and CHOL, exhibited negative correlations (Fig. 11F).
3.10. Molecular docking analysis
To identify potential small molecule compounds interacting with MRAS, molecular docking analysis was conducted. This method aids in discovering therapeutic targets and developing innovative treatments. Among 2614 anticancer small molecules screened using MOE software, the top five scoring compounds were selected for molecular docking: LMP-400 (-17.4463), RGB-286638 (-16.9881), INH 34 (-15.9679), canertinib dihydrochloride (-15.0738), and cytidine 5'-diphosphocholine (-15.0732) (Fig. 12).