Development and Validation of a Robust Apoptosis-Related Prognostic Classifier in Patients with Osteosarcoma

Introduction

Osteosarcoma is the prevailing primary bone tumor with high-ly malignant, which is tend to occur in teenagers. The 5-year survival rate is between 50 percent and 60 percent due to local progression and early metastasis [1]. The therapeutic effect of osteosarcoma patients depends mainly upon the time of diagno¬sis. However, the progression of osteosarcoma is often severe at the time of diagnosed due to the lack of specific symptoms in early stages[2,3]. New biomarkers for osteosarcoma can not only conduce to early screening, diagnosis and predicting prog¬nosis, but also provide a new perspective for the treatment of osteosarcoma.

Apoptosis is one of the main regulative ways of cell death and it can be influenced by genes or directly affected by cytokines. Inhibition of apoptosis is associated with tumorgenesis while ac-tivation of apoptosis can be used as a way of oncotherapy [4, 5]. For example, Chen et al has proved that apoptosis is related to the occurrence of gastric cancer, and Wahba et al has proved that apoptosis is related to the treatment of epithelial ovarian cancer [6,7]. There are many studies focus on the relationship between osteosarcoma and apoptosis by single gene or gene family[8-10]. However, apoptosis is a complex biological process, which in¬evitably involves many regulatory mechanisms [11]. Compre¬hensive studies of multiple pathways and biological processes associated with apoptosis are still lacking.

In this study, all genes were extracted form 153 apoptosis-re¬lated pathways and biological processes, and the datasets of TARGET and GEO (https://www.ncbi.nlm.nih.gov/gds/) were used to construct and validate ARS. ESTIMATE algorithm and CIBERSROT were performed to further investigated the relationship among ARS, immunity and tumor composition. In addi¬tion, the expression levels of the immune checkpoint in different ARS groups were also be analyzed. Finally, a nomogram was constructed and tested to facilitate clinical application. Our find- ings put forward a new perspective for treatment and exploring the potential genes of apoptosis in osteosarcoma.

Materials and Methods

Data Acquisition and Processing FPKM and clinicopathological information of 89 samples with osteosarcoma were downloaded from TARGET datasets of UCSC Xena (http://xena.ucsc.edu/; Accessed 6 November 2021) [12]. The samples without survival information were removed. Ultimately, the training cohort with 84 osteosarcoma samples and corresponding clinicopathological and prognostic informa- tion were brought into the analysis.

The validation cohort with 53 osteosarcoma samples was de¬rived from microarray dataset GSE21257 of GEO DataSets. This dataset was produced by Illumina human-6 v2.0 expression beadchip (using nuIDs as identifier) and contained survival information of each sample.

The data processing is shown in Figure 1. All microarray  RNA-seq data included in our study were normalized and log2 transformed. 

Figure 1: The flow chart of our study.

Apoptosis-Related Pathways and Biological Processes

153 apoptosis-related pathways or biological processes were extracted from hallmark gene sets (H collection), curated gene sets (C2 collection), and ontology gene sets (C5 collection) in MSigDB (http://software.broadinstitute.org/gsea/index.jsp) [13, 14]. All these pathways, biological processes, and their corresponding genes are fused into Table S1.

Signature Establishment in Training Cohort

All genes in apoptosis-related pathways or biological processes were extracted and prognosis-related genes were filtered by the univariate Cox regression with the same criterion of P-value < 0.05. Then, LASSO regression model was performed to find the most robust prognosis-related markers. ARS of each sample was established by the formula:

ARS = ∑ 1 Coeficient(mRNAi) × Expression(mRNAi)

The Power of ARS in Different Cohorts

The training cohort was split into a high-score group and a low-score group with the median value of ARS. The survival state and ARS in different group are depicted by risk factor curve, survival status scatter plot, and Kaplan-Meier curve. To verify the predictive capacity of ARS in the training cohort, AUC of 3- and 5-years were calculated and plotted separately using the survivalROC package. And then, the reliability and applicability of ARS were further validated in the validation cohort.

Relationship between ARS and Clinicopathological Features

To explore the relationship between ARS and different clinico-pathological features containing age, sex, location and metas-tasis, the Kaplan-Meier curve was plotted to evaluate the pre-dictive capacity of ARS in each subgroup based on the median value. ESTIMATE, Immune Profile and Immune Checkpoint Molecules To estimate the proportion of stromal and immune cells in each sample, ESTIMATE algorithm (https://bioinformatics.mdander-son.org/estimate/) was used to calculate the stromal score and immune score respectively [15]. To determine the depth of immune cell infiltration, the CIBERSORT package was performed to calculate the expression of 22 immune cells in different ARS groups [16]. Besides, the expression levels of the immune checkpoint (PD-1, PD-L1, and CTLA4) in different ARS groups were compared.

Construction of a Prognostic Nomogram

All TARGET osteosarcoma samples were used to establish the nomogram. The univariate and multivariate cox regression anal¬yses were performed to selected independent risk factors from ARS and clinicopathological features by the criterion of P-value < 0.05. The nomogram was constructed by the “RMS” package and the stability of the nomogram performance was evaluated by the calibration curve and tROC.

Statistical Analysis

All statistical analyses were performed using the R software (www.r-project.org). Clinicopathological characteristics were compared within the training group and validation group using the Chi-square test, Fisher exact probability test, and Studentâ??s t test.

Results

Data Processing

The expression data and clinicopathological information of samples in training cohort and validation cohort were download respectively. The clinicopathological characteristics of the two cohorts are detailed in Table 1. All 153 apoptosis-related pathways or biological processes were downloaded and corresponding 1807 nonredundant genes were used to take intersection with training and validation cohorts. Ultimately, 1476 mutual genes were identified.

                                   Table 1: characteristics of training and validation cohorts

characteristics	No.(%)		P-value
	TARGET (n=84)	GSE21257 (n=53)
age(mean±sd)	14.99±4.82	18.71±12.20	0.013
sex(%)			0.440
Female	37(44.0)	19(35.8)
Male	47(56.0)	34(64.2)
metastasis state(%)			<0.001
metastasis	63(75.0)	0(0.0)
non-metastasis	21(25.0)	0(0.0)
location(%)			0.145
foot	76(90.5)	44(83.0)
hand	6(7.1)	8(15.1)
pelvis	2(2.4)	0(0.0)
unkown	0(0.0)	1(1.9)

The Construction of ARS

Univariate Cox regression analysis was performed for 1476 mutual genes, and 180 genes with P-value < 0.05 in univariate analysis were included in the LASSO regression analysis (Table S2) to construct the ARS. The ten-fold cross-validation was per-formed to determine the penalty parameter (λ) of the model. A total of 22 genes (RPS6, IFITM3, GRN, ATF4, MYC, DYNLL2, G6PD, BNIP3, PTGIS, BCL10, TRIM32, MAGEA3, PDK2, NMNAT1, EN1, TNFRSF11B, PPARG, SYNGAP1, GAL, GRIK2, MCF2, TERT) were included in the LASSO model (Figure 2). The coefficient of each gene was shown in Table S3.

 Figure 2: The LASSO regression model was used to identify the most robust markers, with an optimal λ value.

The Power of ARS in Different Cohorts In training cohort, the high-score group exhibited worse overall survival data as compared to the low-score group, with P-value < 0.01 (Figure 3E). Meanwhile, the distribution patterns of risk scores and survival status were plotted (Figure 3A and 3C). With the increasing ARS score, the overall survival time decreased and mortality increased. The 3- and 5-years AUC (AUC = 0.937, 0.947) of ARS was shown in Figure 3G, calculated by the “sur-vivalROC” package.

These results indicate that the higher ARS represent the worse prognosis, and ARS has a good ability to predict the prognosis of osteosarcoma patients. The same result was also shown in the validation cohort. we calculate the ARS in validation cohort ac¬cording to the LASSO formula. The high-score group exhibited worse overall survival as compared with the low-score group (P-value < 0.05; Figure 3F). The trend of survival time and mor¬tality was the same as in the training cohort (Figure 3B and 3D). We found that the 3- and 5-years AUC (AUC = 0.771, 0.737) of ARS could still accurately predict survival state of osteosarcoma patients in validation cohort

(Figure 3H).

Figure 3: The performance of ARS in different cohorts. (A, B) The distribution of ARS in training and validation cohorts. (C, D) The distributions of overall survival status, overall survival and risk score in training and validation cohorts. (E, F) Kaplan-Meier curves for the overall survival of the high- and low-score groups in training and validation cohorts. (G, H) AUC plots showed that ARS was an accurate variable for survival prediction in training and validation cohorts.

Relationship Between ARS and Clinicopathological Features

To explore the relationship between ARS and clinicopatholog-ical features including age, sex, location and metastasis, we plotted the Kaplan-Meier curve in a different subgroup of clin¬icopathological features. As shown in Figure 4, the high-score group exhibited worse overall survival compared with the low-score group in each subgroup except in location of hand.

Figure 4: The performance of ARS in different clinicopathological including age, gender, location and metastasis state subgroups.

ESTIMATE, Immune Profile and Immune Checkpoint Molecules

The ESTIMATE algorithm revealed the immune (-1439.104 to 2560.632) and stromal scores (-689.9769 to 1927.611) of the training cohort (Figure 5A). It can be seen from the heatmap that there are significant differences in tumor purity, immune score and stromal score between the two groups. To further elucidate the relationship between ARS and immune/stromal score, both of the scores in ARS subgroups were compared. The Figure 5B-C shows that both immune scores and stromal scores in the low-score group were higher than that in the high-score group. This verified that immune score and ratio of tumor to stroma were significantly associated with ARS.

Figure 5: (A) A heat map of the estimate results. (B, C) The expression levels of immune score and stromal score in high- and low-score groups of ARS. CIBERSORT package was performed to calculate the expres-sion of 22 immune cells in different ARS groups. The result showed that NK cells resting is significantly more prevalent in the high-score group, while Macrophages M1 is significantly more prevalent in the low-score group (Figure 6). To analyze the expression levels of immune checkpoint proteins in ARS subgroup, the expression of PD1, PDL1, and CTLA4 was calculated. The samples with low-score had higher expres¬sion of CTLA4 and PDL1 (Figure 6B-C), suggesting that patients with low-score ARS may respond better to immune check¬point inhibitors targeting CTLA4 and PDL1.

Figure 6: (A) Violin plot showed the different proportions of tumor-infiltrating cells between high- and low-score groups. (B-D) The expression levels of immune checkpoint molecules in high- and low-score groups.

Construction of a Prognostic Nomogram

The clinicopathological features and ARS of training cohort were used to construct a prognostic nomogram. Ultimately, the metastasis state and ARS were identified as the independent fac-tors for prognosis of osteosarcoma according to univariate and multivariate cox regression analysis by the criterion of P-value < 0.05. A prognosis-related nomogram was constructed by the “rms” package with independent factors (Figure 7C). The tAUC of clinicopathological features, ARS and nomogram were shown in Figure 7B. The figure shown that the prediction ability of ARS (AUC = 0.887, 0.976, 0.937, 0.937, 0.947) and nomogram (AUC = 0.932, 0.984, 0.939, 0.939, 0.948) was higher than all other clinicopathological features. Then, the calibration curves of 1-, 3- and 5-year also shown that the prediction ability of the no¬mogram was very stable in different time points (Figure 7D-F).

Figure 7: (A) tROC analysis showed that ARS was an accurate variable for survival prediction in training cohort. (B) Univariate and multivariate cox regression analysis indicated that ARS was the independent risk factor among various features. (C) A nomogram was constructed to quantify risk assessment for individual patients. (D-F) Calibration analysis indicated a high accuracy of survival prediction in 1-, 3- and 5-year.

Discussions

Although the therapeutic effect of surgery combined with che-motherapy is satisfactory in most patients with osteosarcoma, the treatment regimen did not significantly improve overall survival in those patients with metastasis or recurrence and the overall 5-year survival rate of those patients is only 20 percent [17-19]. It has been proved that preoperative response to chemothera¬py can predict the overall survival of patients, but its accuracy can be further improved [20]. It is still of great significance to find new genes as therapeutic targets and prognostic indicators. Currently, biotherapy is a main direction of tumor therapy, and inducing tumor cell apoptosis is the foundation for biotherapy to achieve therapeutic effect. There are also many potential mechanisms of apoptosis, such as death receptor-dependent pathway, mitochondrial-dependent pathway and caspase-independent apoptosis [21].

In this study, 153 apoptosis-related pathways and biological pro¬cesses and corresponding genes were used to explore the role of apoptosis in osteosarcoma.

In order to meet the requirements of subsequent analysis, the TARGET cohort was used to construct ARS consisting of 22 apoptosis-related genes by the univariate cox regression and LASSO regression. Next, we found that ARS had a very good predictive accuracy for overall survival of osteosarcoma in both training and validation cohorts by using the Kaplan-Meier curve and AUC.

For the sake of illustrating the applicability of ARS in different situations, subgroup analysis was conducted. The performance of ARS was very stable in each subgroup of age, sex, metasta¬sis and location, except for the hand subgroup. The small sam¬ple size of the subgroup may seriously affect the results of Ka-plan-Meier curve. Nevertheless, the above results can be used to demonstrate that ARS has the ability to distinguish the prognosis of patients with osteosarcoma in different subgroups.

Recently, great progress has been made in immunotherapy. Sip-uleucel-t has been proved to treat castration-resistant prostate cancer, and the HPV vaccine can prevent and treat infections caused by the HPV virus as well as cancers induced by it [22,23]. ESTIMATE algorithm was performed to explore the relation¬ship among ARS, immune infiltration and proportion of stroma in osteosarcoma. We found a correlation between ARS and im¬mune infiltration as well as proportion of stroma in osteosarco¬ma. CIBERSORT was used to further analyze which immune cell infiltrates were different in different ARS group. The result showed that NK cells resting is significantly more prevalent in the high-score group, while Macrophages M1 is significantly more prevalent in the low-score group. Natural killer (NK) cells have significant capability in tumor immune-surveillance but the penetration of NK cells in tumor is a huge obstacle for cancer immunotherapy. Verhoeven et al had shown that NK cells can recognize and lyse Ewing sarcoma cells through NKG2D and DNAM-1 receptor dependent pathways. In addition, cetuximab can enhance the the cytolytic activity of resting NK cells in os-teosarcoma[24-26]. These results suggest that NK cells may contribute to anticancer activity in osteosarcoma with high-score ARS.

Macrophages play an important role in cancer develop¬ment and metastasis and Macrophages M1 also can phagocytose tumor cells. It can be used as drug carriers for tumor therapy because it can anchor tumor cells. Cersosimo et al had shown that increased infiltration of macrophages M2 is associated with the metastasis and prognosis of osteosarcoma[27,28]. According to the results of our study, Macrophages M1 may contribute to anticancer activity in osteosarcoma with low-score ARS. We also analyzed the relationship between immune checkpoint inhibitors and ARS, and the results showed that the expression of CTLA4 and PDL1 was higher in the low-score group, indicating that patients with low-score ARS may benefit from the treatment of immune checkpoint inhibitors.

Among the 22 genes of ARS, some genes, such as IFITM3, MYC, DYNLL2 and G6PD, had been proved associated with the occurrence and development of osteosarcoma while other genes, such as RPS6, GRN, GRN and PTGIS, has not been stud¬ied in osteosarcoma and further research is needed [29-33].

Several limitations exist in this study. On the one hand, the number of osteosarcoma samples is not enough in both the training and validation groups, and larger-scare samples are needed to verify the stability of the model. On the other hand, the rela¬tionship between ARS and immune cells or immune checkpoints requires further experimental verification.

Conclusion

In conclusion, a prognostic classifier based on 22 apoptosis-related genes and a nomogram was constructed to predict the overall survival of patients with osteosarcoma. The 22 genes of ARS provide a new direction for exploring the mechanisms of apoptosis in osteosarcoma. ARS also provides a reference for selecting suitable patients for immunotherapy and targeted therapy.

Declaration

Ethics approval and consent to participate

The study is the bioinformatics analysis article. The data in our study were obtained from the free online databases and the ethical approval was not necessary. There are no animal and human experiments involved in the study. There are no human subjects in the article and informed consent was not necessary.

Consent for Publication

Not applicable

Availability of Data and Materials

All data used in the study can be downloaded from UCSC Xena (http://xena.ucsc.edu/; Accessed 6 November 2021), GEO Data-Sets (https://www.ncbi.nlm.nih.gov/gds/; Accessed 6 September 2021) and the Molecular Signatures Database (http://software. broadinstitute.org/gsea/index.jsp; Accessed 9 October 2021).

Competing interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this ar¬ticle.

Funding

Not applicable

Authors Contributions

Yinquan Zhang and Zhifeng Zhang performed the study conceptualization.

Zhifeng Zhang and Yi Wang performed the data collection and analysis.

Fengmei Chen performed writing the original manuscript draft.

Zhifeng Zhang contributed to the manuscript review and editing. Yi Wang contributed to the drawing pictures.

Zhifeng Zhang, Yi Wang and Fengmei Chen contributed equally to this work.

Acknowledgement

Not applicable.

Abbreviations

FPKM: The fragments per kilobase per million mapped reads; TARGET: Therapeutically Applicable Research to Generate Ef-fective Treatments; GEO: Gene Expression Omnibus; MSigDB: Molecular Signatures Database; ARS: apoptosis-related risk score; LASSO: Least Absolute Shrinkage and Selection Oper-ator; AUC: area under the curve; tROC: time-dependent receiv¬er operating characteristic; ICI: immune checkpoint inhibitors; PD1: programmed cell death protein 1; PDL1: programmed death-ligand 1; CTLA4: cytotoxic T-lymphocyte-associated protein 4.

References

1. Li, J., Yang, Z., Li, Y., Xia, J., Li, D., Li, H., ... & Zhang, Y.(2016). Cell apoptosis, autophagy and necroptosis in osteo­sarcoma treatment. Oncotarget, 7(28), 44763.
2. Ritter, J., & Bielack, S. S. (2010). Osteosarcoma. Annals of oncology, 21, vii320-vii325.
3. Kansara, M., Teng, M. W., Smyth, M. J., & Thomas, D. M. (2014). Translational biology of osteosarcoma. Nature Re­views Cancer, 14(11), 722-735.
4. Johnstone, R. W., Ruefli, A. A., & Lowe, S. W. (2002). Apoptosis: a link between cancer genetics and chemothera­py. Cell, 108(2),153-164.
5. Maiuri, M. C., Zalckvar, E., Kimchi, A., & Kroemer, G. (2007). Self-eating and self-killing: crosstalk between au­tophagy and apoptosis. Nature reviews Molecular cell biol­ogy, 8(9), 741-752.
6. Chen, X., Zhao, Y., Luo, W., Chen, S., Lin, F., Zhang, X., ... & Liang, G. (2020). Celastrol induces ROS-mediated apop­tosis via directly targeting peroxiredoxin-2 in gastric cancer cells. Theranostics, 10(22), 10290.
7. Wahba, J., Natoli, M., Whilding, L. M., Parente-Pereira, A. C., Jung, Y., Zona, S., ... & Ghaem-Maghami, S. (2018). Chemotherapy-induced apoptosis, autophagy and cell cycle arrest are key drivers of synergy in chemo-immunotherapy of epithelial ovarian cancer. Cancer Immunology, Immuno­therapy, 67(11), 1753-1765.
8. Liu, K., Ren, T., Huang, Y. I., Sun, K., Bao, X., Wang, S., ...& Guo, W. (2017). Apatinib promotes autophagy and apop­tosis through VEGFR2/STAT3/BCL-2 signaling in osteo-sarcoma. Cell Death & Disease, 8(8), e3015-e3015.
9. Wang, Z., Yin, F., Xu, J., Zhang, T., Wang, G., Mao, M., ... & Cai, Z. (2019). CYT997 (Lexibulin) induces apoptosis and autophagy through the activation of mutually reinforced ER stress and ROS in osteosarcoma. Journal of Experimental & Clinical Cancer Research, 38(1), 1-15.
10. He, C., Liu, C., Wang, L., Sun, Y., Jiang, Y., & Hao, Y. (2019). Histone methyltransferase NSD2 regulates apop­tosis and chemosensitivity in osteosarcoma. Cell death & disease, 10(2), 1-13.
11. Hengartner, M. O. (2000). The biochemistry of apoptosis. Nature, 407(6805), 770-776.
12. Goldman, M. J., Craft, B., Hastie, M., Repečka, K., Mc-Dade, F., Kamath, A., ... & Haussler, D. (2020). Visualizing and interpreting cancer genomics data via the Xena plat­form. Nature biotechnology, 38(6), 675-678.
13. Liberzon, A., Birger, C., Thorvaldsdóttir, H., Ghandi, M., Mesirov, J. P., & Tamayo, P. (2015). The molecular signa­tures database hallmark gene set collection. Cell systems, 1(6), 417-425.
14. Liberzon, A., Subramanian, A., Pinchback, R., Thorvalds­dóttir, H., Tamayo, P., & Mesirov, J. P. (2011). Molecular signatures database (MSigDB) 3.0. Bioinformatics, 27(12),1739-1740.
15. Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., ... & Mesirov, J. P. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proceed­ings of the National Academy of Sciences, 102(43), 15545-15550.
16. Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J.,Feng, W., Xu, Y., ... & Alizadeh, A. A. (2015). Robust enu­meration of cell subsets from tissue expression profiles. Na­ture methods, 12(5), 453-457.
17. Link, M. P., Goorin, A. M., Miser, A. W., Green, A. A., Pratt,C. B., Belasco, J. B., ... & Vietti, T. J. (1986). The effect of adjuvant chemotherapy on relapse-free survival in patients with osteosarcoma of the extremity. New England Journal of Medicine, 314(25), 1600-1606.
18. Meyers, P. A., Healey, J. H., Chou, A. J., Wexler, L. H.,Merola, P. R., Morris, C. D., ... & Gorlick, R. (2011). Ad­dition of pamidronate to chemotherapy for the treatment of osteosarcoma. Cancer, 117(8), 1736-1744.
19. Bernthal, N. M., Federman, N., Eilber, F. R., Nelson, S. D., Eckardt, J. J., Eilber, F. C., & Tap, W. D. (2012). Long-term results(> 25 years) of a randomized, prospective clinical trial evaluating chemotherapy in patients with high-grade, operable osteosarcoma. Cancer, 118(23), 5888-5893.
20. Collins, M., Wilhelm, M., Conyers, R., Herschtal, A., Whel­an, J., Bielack, S., ... & Thomas, D. M. (2013). Benefits and adverse events in younger versus older patients receiving neoadjuvant chemotherapy for osteosarcoma: findings from a meta-analysis. Journal of Clinical Oncology, 31(18), 2303-2312.
21. Vermeulen, K., Van Bockstaele, D. R., & Berneman, Z. N. (2005). Apoptosis: mechanisms and relevance in cancer. Annals of hematology, 84(10), 627-639.
22. Markowitz, L. E., Liu, G., Hariri, S., Steinau, M., Dunne,E. F., & Unger, E. R. (2016). Prevalence of HPV after in­troduction of the vaccination program in the United States. Pediatrics, 137(3).
23. Kantoff, P. W., Higano, C. S., Shore, N. D., Berger, E. R.,Small, E. J., Penson, D. F., ... & Schellhammer, P. F. (2010). Sipuleucel-T immunotherapy for castration-resistant pros­tate cancer. New England Journal of Medicine, 363(5), 411-422.
24. Pahl, J. H., Ruslan, S. E. N., Buddingh, E. P., Santos, S. J.,Szuhai, K., Serra, M., ... & Lankester, A. C. (2012). An-ti-EGFR antibody cetuximab enhances the cytolytic activity of natural killer cells toward osteosarcoma. Clinical cancer research, 18(2), 432-441.
25. Verhoeven, D. H., de Hooge, A. S., Mooiman, E. C., Santos,S. J., Monique, M., Gelderblom, H., ... & Lankester, A. C. (2008). NK cells recognize and lyse Ewing sarcoma cells through NKG2D and DNAM-1 receptor dependent path­ways. Molecular immunology, 45(15), 3917-3925.
26. Valipour, B., Velaei, K., Abedelahi, A., Karimipour, M., Da-rabi, M., & Charoudeh, H. N. (2019). NK cells: An attrac­tive candidate for cancer therapy. Journal of cellular physi­ology, 234(11), 19352-19365.
27. Xia, Y., Rao, L., Yao, H., Wang, Z., Ning, P., & Chen, X. (2020). Engineering macrophages for cancer immunothera­py and drug delivery. Advanced Materials, 32(40), 2002054.
28. Cersosimo, F., Lonardi, S., Bernardini, G., Telfer, B., Mandelli, G. E., Santucci, A., ... & Giurisato, E. (2020).Tumor-associated macrophages in osteosarcoma: from mechanisms to therapy. International journal of molecular sciences, 21(15), 5207.
29. Wang, X., Chen, K., & Zhao, Z. (2020). LncRNA OR3A4regulated the growth of osteosarcoma cells by modulating the miR-1207-5p/G6PD signaling. OncoTargets and thera­py, 13, 3117.
30. Salas, S., Jézéquel, P., Campion, L., Deville, J. L., Chibon, F., Bartoli, C., ... & Bouvier, C. (2009). Molecular charac­terization of the response to chemotherapy in convention­al osteosarcomas: predictive value of HSD17B10 and IF-ITM2. International journal of cancer, 125(4), 851-860.
31. Feng, W., Dean, D. C., Hornicek, F. J., Spentzos, D., Hoff-man, R. M., Shi, H., & Duan, Z. (2020). Myc is a prog­nostic biomarker and potential therapeutic target in osteo­sarcoma. Therapeutic advances in medical oncology, 12, 1758835920922055.
32. Shen, Y., Zhao, S., Wang, S., Pan, X., Zhang, Y., Xu, J.,... & Zhou, D. (2019). S1P/S1PR3 axis promotes aerobic glycolysis by YAP/c-MYC/PGAM1 axis in osteosarcoma. EBioMedicine, 40, 210-223.
33. Jiang, J., Liu, D., Xu, G., Liang, T., Yu, C., Liao, S., ... & Zhan, X. (2021). TRIM68, PIKFYVE, and DYNLL2: Thenovel autophagy-and immunity-associated gene biomark­ers for osteosarcoma prognosis. Frontiers in Oncology, 11, 1238.