Development and Verification of a Hypoxic Gene Signature for Predicting Prognosis, Immune Microenvironment, and Chemosensitivity for Osteosarcoma

Objective: Hypoxic tumors contribute to local failure and distant metastases. Nevertheless, the molecular hallmarks of hypoxia remain ill-defined in osteosarcoma. Here, we developed a hypoxic gene signature in osteosarcoma prognoses. Methods: With the random survival forest algorithm, a prognostic hypoxia-related gene signature was constructed for osteosarcoma in the TARGET cohort. Overall survival (OS) analysis, receiver operating characteristic (ROC) curve, multivariate cox regression analysis, and subgroup analysis were utilized for assessing the predictive efficacy of this signature. Also, external validation was presented in the GSE21257 cohort. GSEA was applied for signaling pathways involved in the high- and low-risk samples. Correlation analyses between risk score and immune cells, stromal/immune score, immune checkpoints, and sensitivity of chemotherapy drugs were performed in osteosarcoma. Then, a nomogram was built by integrating risk score, age, and gender. Results: A five-hypoxic gene signature was developed for predicting survival outcomes of osteosarcoma patients. ROC curves confirmed that this signature possessed the well predictive performance on osteosarcoma prognosis. Furthermore, it could be independently predictive of prognosis. Metabolism of xenobiotics by cytochrome P450 and nitrogen metabolism were activated in the high-risk samples while cell adhesion molecules cams and intestinal immune network for IgA production were enriched in the low-risk samples. The low-risk samples were characterized by elevated immune cell infiltrations, stromal/immune scores, TNFRSF4 expression, and sensitivity to cisplatin. The nomogram accurately predicted 1-, 3-, and 5-years survival duration. Conclusion: These findings might offer an insight into the optimization of prognosis risk stratification and individualized therapy for osteosarcoma patients.


INTRODUCTION
Osteosarcoma represents the main primary bone malignancy, occurring in 2.4% of children and adolescents . This disease typically presents in distal femoral metaphysis, followed by the upper extremities, pelvis, and the like (Kelley et al., 2020). Surgical resection and chemotherapy distinctly prolong survival duration of patients with local osteosarcoma. Non-metastatic patients reach the 60-70% of the 5-years survival rate following diagnosis . However, osteosarcoma is characterized by high metastatic and recurrent risks. For metastatic or relapsed subjects, the 5-years survival rate is only 20% (Shi et al., 2020). It is urgent to probe biomarkers for assisting clinicians to predict survival outcomes as well as offer a basis upon personalized medicine.
The intratumoral hypoxia is a noticeable feature in solid tumors, and hypoxia is in relation to aggressive phenotypes and metastases (Song et al., 2020). For osteosarcoma, hypoxia is involved in facilitating excessive tumor proliferative and invasive capacities through inducing serial molecular reactions to promote the formation of the neoplastic microenvironment . Under hypoxic conditions, various genes can be activated in osteosarcoma cells. For instance, hypoxia induces chemoresistance through down-regulation of SKA1 expression in osteosarcoma (Ma et al., 2017). Hypoxia-induced WSB1 accelerates invasion and metastasis of osteosarcoma cells (Cao et al., 2015). Nevertheless, no hypoxic gene signature has been established yet for osteosarcoma patients' prognosis. Transcriptomic analysis may uncover the entire expression patterns as well as identify cancer-related markers. Here, this study effectively developed a robust prognosis-related hypoxic gene signature for osteosarcoma, which might provide a basis for an in-depth understanding of risk stratification and therapeutic effects of osteosarcoma patients.

Data Collection
The gene expression profiles of 84 cases of osteosarcoma patients with complete follow-up information were downloaded from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET; https://ocg.cancer.gov/programs/target) database, which were utilized as the training set. The GSE21257 dataset containing expression profiling of 53 osteosarcoma subjects was retrieved from the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/) (Buddingh et al., 2011), which was used as the external validation set. Genes up-regulated in response to hypoxia were retrieved from the "HALLMARK_HYPOXIA" gene set in the Molecular Signatures database v7.2 with the gene set enrichment analysis (GSEA) software (Subramanian et al., 2005). Figure 1 depicts the workflow of this study.

Establishment of a Prognostic Gene Signature
Survival-related hypoxic genes with p < 0.05 were screened for osteosarcoma samples from the TARGET cohort with the univariate cox regression analyses. Then, the importance of these prognostic genes was ranked using the random survival forest algorithm. The number of Monte Carlo iterations was set as 100 and the number of steps forward was set as 5. The genes with relative importance >0.25 were feature factors. Multivariate cox regression analyses were then presented. The risk scoring model was established as follows: risk score coefficient of gene 1 * expression of gene 1 + coefficient of gene 2 * expression of gene 2 +. . .+ coefficient of gene n * expression of gene n. The specimens were classified into high and low risk groups on the basis of the mean value of risk scores. Expression patterns of genes in the signature were visualized into a heat map. Survival status was evaluated between the two groups. Kaplan-Meier curves of overall survival (OS) were then depicted between groups and survival differences were compared via log-rank test. The time-dependent receiver operating characteristic (ROC) was analyzed for examining the sensitivity and specificity of the gene signature for predicting OS. The area under the curve (AUC) was determined for assessing the prognostic accuracy. Furthermore, the AUC of this signature was compared with other gene signatures developed by Wen C. et al. (2020), Wu Z.-l. et al. (2020, Xiao et al. (2020) and Zhang et al. (Yiqi et al., 2020). Univariate and multivariate cox regression analyses were carried out based on risk score, age, and gender as variables. p values, hazard ratio (HR) as well as 95% confidence interval (CI) were calculated. Moreover, survival analysis of each gene in this signature was presented for osteosarcoma subjects in the TARGET cohort.

Subgroup Analysis
Osteosarcoma samples were classified into different subgroups according to age (<16 and ≥16) and gender (female and male). In each subgroup, OS differences were compared between high and low risk samples with Kaplan-Meier curves and logrank tests.

GSEA
Osteosarcoma specimens in the TARGET cohort were classified into high and low risk groups. GSEA was carried out to explore the underlying pathways involved in the two groups. The | normalized enrichment score (NES)| > 1.5 and false discovery rate (FDR) < 0.05 were used to determine statistical significance.

Single-Sample Gene Set Enrichment Analysis
The enrichment scores of activated B cells, activated CD4 T cells, activated CD8 T cells, central memory CD4 T cells, central memory CD8 T cells, effector memory CD4 T cells, effector memory CD8 T cells, gamma delta T cells, immature B cells, memory B cells, regulatory T cells, T follicular helper cells, type 1 T helper cells, type 17 T helper cells, type 2 T helper cells, activated dendritic cells, CD56bright natural killer cells, CD56dim natural killer cells, eosinophils, immature dendritic cells, macrophages, mast cells, MDSCs, monocytes, natural killer cells, natural killer T cells, neutrophils, and plasmacytoid dendritic cells were determined in high and low risk samples from the TARGET cohort with the ssGSEA algorithm.

Estimation of Stromal and Immune Cells in Malignant Tumours Using Expression Data
With the ESTIMATE algorithm, stromal and immune scores were inferred in high and low risk samples from the TARGET cohort based on the gene expression profiling (Yoshihara et al., 2013).
Correlation Between Risk Score and Immune Checkpoint Expression were assessed in the osteosarcoma specimens from the TARGET cohort.

Assessment of Chemotherapy Drug Sensitivity
Through the Genomics of Drug Sensitivity in Cancer (GDSC) database (www.cancerRxgene.org) (Yang et al., 2013), the IC50 values of chemotherapy drugs (cisplatin, doxorubicin, methotrexate, and paclitaxel) were estimated in each osteosarcoma specimen with the pRRophetic package (Geeleher et al., 2014).

Nomogram
A nomogram was established by integrating gender, age, and risk score with the "rms" package for osteosarcoma prognoses. The calibration plots were depicted for verifying the accuracy of prediction of 1-, 3-, and 5-years survival. Also, the predictive efficacy of this nomogram was evaluated by time-dependent ROC analysis.

Colony Formation Assay
In total, 500 Saos-2 and U2OS cells were planted onto a 6-well dish. After culturing for 2 weeks, the colonies were washed by PBS and dyed with a dyeing solution containing 0.1% crystal violet and 20% methanol. Thereafter, the colonies were counted.

Transwell Assays
Transwell assays were performed for measuring migration and invasion of Saos-2 and U2OS cells utilizing 24-well inserts (8 μm pore, Corning, United States). For migration, Saos-2 and U2OS cells re-suspended in 100 μl serum-free DMEM medium were added to the upper chamber without coating membrane. For invasion, re-suspended Saos-2 and U2OS cells were added to the upper chamber with Matrigel-coated membrane. DMEM plus 600 μl 10% FBS was added to the lower chamber. Migrated or invaded cells from the bottom surfaces were fixed by 4% paraformaldehyde, followed by being stained by 0.1% crystal violet for 20 min. Finally, the cells were investigated and counted under a microscope (Leica, Germany).

Statistical Analysis
All statistics were carried out by applying the R software (version 3.6.5; https://www.R-project.org) and its packages. Comparisons between two groups were presented via Wilcoxon or Student's t tests. Multiple comparisons were presented through one-or two-way ANOVA test. p values <0.05 were set as statistical significances.

Establishment of a Prognostic Hypoxic Gene Signature in Osteosarcoma
In total, 197 hypoxic genes were obtained for this study (Supplementary Table S1). The chromosomal locations of these genes are shown in Figure 2A. Through univariate cox regression analysis, 26 hypoxic genes exhibited significant correlations to prognosis of osteosarcoma patients in the Frontiers in Molecular Biosciences | www.frontiersin.org January 2022 | Volume 8 | Article 705148 TARGET cohort ( Table 1). The most important hypoxic genes related to osteosarcoma prognosis were selected using the random survival forest algorithm ( Table 2). When the number of classification trees was 5, the error rate was the lowest ( Figure 2B). The feature importance ranking of the top five genes is shown in Figure 2C, including TES (chr7: Following multivariate cox regression analysis, a five-hypoxic gene signature was constructed in osteosarcoma ( Figure 2D). The scoring formula was as follows: risk score 0.081751 * EGFR expression + (−0.3013) * CAVIN1 expression +0.812177 * MXI1 expression + (−0.46302) * SDC3 expression + (−0.37064) * TES expression. According to the median value of risk scores, osteosarcoma patients were separated into high-and low-risk subgroups (n 42). The heat map visualized the expression patterns of TES, SDC3, MXI1, CAVIN1, and EGFR in high-and low-risk specimens ( Figure 2E). The proportions of dead patients were separately 50 and 14% in high-and low-risk groups ( Figure 2F). Moreover, high-risk scores were indicative of unfavorable survival outcomes of osteosarcoma patients ( Figure 2G).

The Hypoxic Gene Signature Robustly and Independently Predicts Osteosarcoma Prognosis
To evaluate the predictive efficacy of this hypoxic gene signature on survival outcomes, time-independent ROC curves were conducted and the AUCs of 1-, 3-, and 5-years OS were 0.789, 0.826, and 0.789, respectively ( Figure 3A). In comparison to other prognostic models, this signature had the higher accuracy and sensitivity for predicting prognosis of osteosarcoma patients ( Figure 3B). Through univariate cox regression analyses, risk score displayed a significant association with osteosarcoma prognosis (p < 0.001, HR: 1.083, 95%CI: 1.035-1.132; Figure 3C). Multivariate cox regression analyses confirmed that risk score was independently predictive of clinical outcomes (p < 0.001, HR: 1.083, 95%CI: 1.033-1.136) with age, gender, and risk score as variables ( Figure 3D).

External Validation of the Predictive Performance of Hypoxic Gene Signature on Osteosarcoma Prognosis
With the same formula, the risk score of each subject in the GSE21257 cohort was determined. Figure 4A depicted the expression of MXI1, SDC3, TES, EGFR, and CAVIN1 in high and low risk osteosarcoma specimens. The percentages of dead patients were 57 and 26% in the high-and low-risk groups, respectively ( Figure 4B). Furthermore, unfavorable survival outcomes were found in high-risk subjects (p 1.757e-02; Figure 4C), which exhibited the consistency with the results in the TARGET cohort. The AUCs of 1-, 3-, and 5-years survival were 0.668, 0.679, and 0.746, respectively ( Figure 4D). Under univariate cox regression analyses, the risk score was in relation to desirable clinical outcomes of osteosarcoma (p 0.005, HR 2.467, and 95% CI: 1.314-4.631; Figure 4E). Nevertheless, age and gender were not significantly correlated to osteosarcoma prognosis. Multivariate cox regression analyses demonstrated that this signature was capable of independently predicting prognosis of osteosarcoma patients (p 0.005, HR 2.485, and 95% CI: 1.317-4.689; Figure 4F). Subgroup analysis was applied to evaluate the predictive performance of the gene signature on osteosarcoma prognosis. First, patients from the TARGERT cohort were classified into age <16 and ≥16 subgroups. We found that high risk scores were in relation to worse survival duration both in two subgroups (p 0.001 and 0.049; Figures 5A,B). Also, both in female (p 0.032; Figure 5C) and male (p 0.002; Figure 5D) subgroups, high-risk scores were indicative of shorter survival outcomes of osteosarcoma patients.

GSEA Identifies Involved Signaling Pathways
GSEA results revealed that altered genes were distinctly enriched in several common pathways. Metabolism of xenobiotics by cytochrome P450 (NES 1.56 and NOM p-value 0.045) and nitrogen metabolism (NES 1.52 and NOM p-value 0.039) were mainly activated in the highrisk osteosarcoma samples ( Figures 6A,B). In the meanwhile, low-risk samples were distinctly related to cell adhesion molecules cams (NES -1.74 and NOM p-value 0.035) and intestinal immune network for IgA production (NES -1.78 and NOM p-value 0.041; Figures 6C,D).

Differences in Immune Status Between High-and Low-Risk Osteosarcoma
The enrichment scores of immune cells were quantified in osteosarcoma specimens by applying the ssGSEA algorithm. Our data demonstrated that higher infiltrations of activated (C) Kaplan-Meier curves of OS between high-and low-risk subjects. The differences were compared by log-rank tests. (D) Time-independent ROC curves of 1-, 3-, and 5-years survival based on risk score. (E, F) Univariate and multivariate cox regression analyses for the associations between prognosis and age, gender, and risk scores.
Frontiers in Molecular Biosciences | www.frontiersin.org January 2022 | Volume 8 | Article 705148 8 B cell s, activated CD8 T cells, central memory CD4 T cells, central memory CD8 T cells, regulatory T cells, type 1 T helper cells, CD56bright natural killer cells, macrophages, MDSC, natural killer cells, and natural killer T cells were found in the low-risk osteosarcoma in comparison to the high-risk osteosarcoma in the TARGET cohort ( Figure 8A). Also, there were higher stromal scores (p 0.00024) as well as immune scores (p 0.045) in the low-risk specimens than the high-risk specimens ( Figure 8B). The correlations between risk scores and immune checkpoints were assessed in osteosarcoma subjects. As a result, risk scores were positively correlated to TNFSF4 ( Figure 8C).

Correlations Between Risk Score and Sensitivity to Chemotherapy Drugs
The IC50 values of chemotherapy drugs were determined in osteosarcoma patients from the TARGET cohort. Lowered IC50 values of cisplatin were detected in the high-risk compared to low-risk specimens (p 0.015; Figure 9A), indicating that high-risk patients might be more sensitive to cisplatin. However, there were no significant differences in IC50 values of doxorubicin, methotrexate, and paclitaxel between groups ( Figures 9B-D).

Construction and Verification of a Personalized Prognostic Nomogram
To offer an applicable quantitative tool for clinicians, a nomogram combining risk score, age, and gender was built for prediction of the probability of 1-, 3-, and 5-years survival of osteosarcoma patients in the TARGET cohort ( Figure 10A). Each patient was determined based on each prognostic index. The higher number of total points was indicative of an undesirable prognosis for osteosarcoma patients. The AUC of OS time was 0.730, demonstrating that this nomogram exhibited high potential upon clinical utility ( Figure 10B). Compared to the ideal model, the nomogram displayed a similar performance at 1-, 3-, and 5-years survival ( Figures  10C-E). This prognostic nomogram model was externally verified in the GSE21257 osteosarcoma dataset ( Figure 10F). The AUC value (0.737) showed that this nomogram possessed the favorable predictive capacity for the OS in osteosarcoma patients ( Figure 10G). Our calibration curves showed excellent agreement between the nomogram prediction and actual observation at 1-, 3-, and 5-years survival in the GSE21257 dataset ( Figures 10H-J). The above data suggested the appreciable reliability of this nomogram.

Verification of the Expression of Hypoxic Genes From the Gene Signature
The expression of hypoxic genes from the gene signature was verified in osteosarcoma and normal cell lines. Our results confirmed that SCD3, EGFR, and MXI1 were remarkably upregulated while CAVIN1 and TES were prominently downregulated in U2OS osteosarcoma cells compared with  Figures 11A-G). We also compared their expression in hFOB1.19 normal cells and Saos-2 osteosarcoma cells. Compared with hFOB1.19 cells, SCD3, EGFR, and MXI1 displayed higher expression while CAVIN1 and TES presented lower expression in Saos-2 cells, indicating that the above genes might participate in osteosarcoma progression ( Figures 11H-N).

Silencing MXI1 Attenuates Proliferation, Migration, and Invasion of Osteosarcoma Cells
The biological functions of MXI1 were further investigated through in vitro experiments. First, the expression of MXI1 was remarkably silenced by its specific siRNAs both in U2OS and Saos-2 cells (Figures 12A,B). Compared with the si-NC group, MXI1 knockdown prominently attenuated proliferation ( Figures 12C,D), migration ( Figures 12E,F), and invasion ( Figures 12G,H) of U2OS and Saos-2 cells.

DISCUSSION
Hypoxia, a hallmark of cancer, may be induced by an insufficient oxygen supply, which is correlated to unfavorable survival outcomes, increased tumor heterogeneity, and drug resistance as well as immune escape (Vito et al., 2020). Here, we developed and verified a hypoxic gene signature for robustly predicting osteosarcoma prognoses, which could become a clinical stratified tool.
The hypoxic gene signature was built with the random survival forest algorithm, which consisted of CAVIN1, EGFR, SCD3, TES, and MXI1. High-risk scores were indicative of undesirable survival outcomes. Our ROCs confirmed the well predictive performance of this signature in osteosarcoma prognoses. The multivariate cox regression analyses were indicative of the predictive independency of the signature. Previously, a 28-gene hypoxic signature was proposed, which possessed the robust and independent prognostic efficacy in prostate cancer (Yang L. et al., 2018). Furthermore, a four-hypoxic gene model served as a prognostic factor for lung adenocarcinoma, which might predict the immunotherapeutic responses (Mo et al., 2020). The signaling pathways involved in risk scores were analyzed in osteosarcoma subjects. Metabolism of xenobiotics by cytochrome P450 and nitrogen metabolism were activated in the high-risk samples while cell adhesion molecules cams and intestinal immune network for IgA production were enriched in the low-risk samples. Evidence suggests that intestinal immune network for IgA production pathway is associated with cell proliferation and migration of hepatocellular carcinoma cells (Yang Z. et al., 2018).
However, no study has reported the role of this pathway in osteosarcoma progression. Our study indicated that low-risk samples might have enhanced immune activation compared with high-risk samples. In the five-gene signature, we separately evaluated the prognostic potential of CAVIN1, EGFR, SCD3, TES, and MXI1 in osteosarcoma. Among them, up-regulation of CAVIN1, EGFR, SCD3, and TES was related to favorable survival outcomes while MXI1 up-regulation was associated with undesirable prognoses. Previously, CAVIN1 emerges as a novel tumor suppressor in Ewing sarcoma The darker the color, the stronger the correlation. Ns: not significant; *p < 005; **p < 0.01; ***p < 0.001.
Frontiers in Molecular Biosciences | www.frontiersin.org January 2022 | Volume 8 | Article 705148 12 (Huertas-Martínez et al., 2017). Overactivity of EGFR is frequently observed in osteosarcoma and becomes an underlying therapeutic target (Ji and He, 2019). MXI1 is a key transcription factor involved in osteosarcoma progression (Xu et al., 2017). Additionally, our experimental results confirmed that SCD3, EGFR, and MXI1 were remarkably up-regulated while CAVIN1 and TES were prominently down-regulated in osteosarcoma cells compared with normal cells. The above findings highlighted the implications of these genes in osteosarcoma prognoses and progression. Our in vitro experiments demonstrated that silencing MXI1 attenuated proliferation, migration, and invasion of osteosarcoma cells, confirming the important role of MXI1 in osteosarcoma progression. Thus, our study is novel and better suited than previously published prognostic protocols for osteosarcoma (Wu Z.-l. et al., 2020;Qiu et al., 2021).
Cancer progression is highly correlated with the physiological state of osteosarcoma immune microenvironment that consists of innate as well as adaptive immune cells (Roma-Rodrigues et al., 2019).
Osteosarcoma possesses a unique tumor microenvironment triggered by distinct molecular characteristics (Dyson et al., 2019). The immune microenvironment profiles are preponderant on prognoses and anti-tumor therapeutic efficacy. Immunotherapy emerges as a promising therapeutic strategy against osteosarcoma . Nevertheless, only minor subjects experience the desired anti-osteosarcoma immunity as well as impressive clinical response. For maximizing treatment effects of osteosarcoma subjects, it is a priority to alter the osteosarcoma immune microenvironment (Riera-Domingo et al., 2020). Hypoxia participates in anti-tumor immune response through lowering cytolytic and migratory activation of effector cells and the production and release of effector cytokines as well as triggering immunosuppressive cells, production and release of immunosuppressive cytokines, and expressions of immune checkpoint inhibitors (Multhoff and Vaupel, 2020). Herein, we  found that there were higher infiltrations of activated B cells, activated CD8 T cells, central memory CD4 T cells, central memory CD8 T cells, regulatory T cells, type 1 T helper cells, CD56bright natural killer cells, macrophages, MDSC, natural killer cells, and natural killer T cells in the low-risk osteosarcoma than the high-risk osteosarcoma, indicating that this signature might reflect the immune microenvironment of osteosarcoma for bench observations. We noted that immunosuppressive cells and immuno-promoting cells were all significantly activated in low-risk groups, indicating the complex interactions between immunosuppressive cells and immuno-promoting cells during osteosarcoma progression. Infiltrating stromal and immune cells constitute the main fractions of the tumor microenvironment. Previous research has demonstrated that high stromal or immune scores indicated favorable survival outcomes of osteosarcoma (Hong et al., 2020). Their scores were determined in osteosarcoma tissues. The increase in stromal and immune scores was detected in the low-risk osteosarcoma specimens. Limited clinical activity of immune checkpoint inhibitors is investigated in osteosarcoma subjects (Wu C.-C. et al., 2020). Hence, it is of significance to gain the immunogenic potential of osteosarcoma. TNFRSF4 possesses potential as a target upon cancer immunotherapy (Buchan et al., 2018). Our data showed that risk score exhibited a positive correlation to immune checkpoint TNFRSF4 in osteosarcoma.
Chemotherapy resistance is the main issue in osteosarcoma therapy, which leads to undesirable prognoses (Prudowsky and Yustein, 2020). Alleviating hypoxia through tumor reoxygenation sensitizes the chemotherapy in osteosarcoma (Fu et al., 2021). Currently, there is a lack of clinically relevant molecular biomarkers that are predictive of the responses to chemotherapies.
Cisplatin represents the first-line chemotherapy drug regarding osteosarcoma treatment, but chemoresistance limits the effectiveness of cisplatin (Wen J.-F. et al., 2020). Our data demonstrated that low-risk patients were more sensitive to cisplatin, indicating that this risk score might aid in overcoming cisplatin chemoresistance. The nomogram, as an appropriate clinical tool, has been widely utilized for quantitatively determining a person's prognosis in the clinical setting through integration of several prognostic factors (Wu and Zhang, 2020). Here, the nomogram model combined risk score, age, and gender was established for predicting the probability of 1-, 3-, and 5-years survival time in osteosarcoma subjects. By calculating all the points, the nomogram yielded the numerical possibility of osteosarcoma patients concerning clinical outcomes. Following verification by ROCs and calibration plots, the nomogram model could be well predictive of 1-, 3-, and 5-years survival possibilities of osteosarcoma subjects. Nevertheless, the predictive efficacy of this nomogram requires to be verified in a larger osteosarcoma cohort.

CONCLUSION
Collectively, we developed a hypoxic gene signature for risk stratification in osteosarcoma. After external verification, this risk score may be independently predictive of patients' survival outcomes. Also, the risk score was in relation to osteosarcoma immune microenvironment and sensitivity to the chemotherapy drug cisplatin. A nomogram that combined the risk score, age, and gender became a promising prognostic classifier for osteosarcoma. Thus, our findings offered insights into personalized prognostication and therapy as well as follow-up scheduling.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
GW conceived and designed the study. FW, JX, MJ, and XJ conducted most of the experiments and data analysis and wrote the manuscript. JL, XL, ZC, JN and ZM participated in collecting