BioMed Research International
Volume 2018 (2018), Article ID 4761064, 10 pages
https://doi.org/10.1155/2018/4761064

Research Article

Identification of Key Genes and miRNAs in Osteosarcoma Patients with Chemoresistance by Bioinformatics Analysis

Binbin Xie,¹ Yiran Li,¹ Rongjie Zhao,¹ Yuzi Xu,² Yuhui Wu,¹ Ji Wang,3,4 Dongdong Xia,⁵ Weidong Han,¹ and Dake Chen⁶

¹Department of Medical Oncology, Sir Run Run Shaw Hospital, School of Medicine, Zhejiang University, Hangzhou, China
²Department of Sports Medicine, School of Medicine, Zhejiang University, Zhejiang, China
³Department of Surgical Oncology, Sir Run Run Shaw Hospital, College of Medicine, Zhejiang University, Hangzhou, Zhejiang 310016, China
⁴Biomedical Research Center and Key Laboratory of Biotherapy of Zhejiang Province, Hangzhou, Zhejiang 310016, China
⁵Orthopedic Department, Ningbo First Hospital, Ningbo 315000, China
⁶Department of Urology, Wenzhou People’s Hospital, Wenzhou, Zhejiang, China

Correspondence should be addressed to Dongdong Xia; xiadongdong@qq.com, Weidong Han; hanwd@zju.edu.cn, and Dake Chen; 61502485@qq.com

Received 23 November 2017; Revised 21 February 2018; Accepted 4 March 2018; Published 22 April 2018

Academic Editor: Jialiang Yang

Copyright © 2018 Binbin Xie et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Chemoresistance is a significant factor associated with poor outcomes of osteosarcoma patients. The present study aims to identify Chemoresistance-regulated gene signatures and microRNAs (miRNAs) in Gene Expression Omnibus (GEO) database. The results of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) included positive regulation of transcription, DNA-templated, tryptophan metabolism, and the like. Then differentially expressed genes (DEGs) were uploaded to Search Tool for the Retrieval of Interacting Genes (STRING) to construct protein-protein interaction (PPI) networks, and 9 hub genes were screened, such as fucosyltransferase 3 (Lewis blood group) (FUT3) whose expression in chemoresistant samples was high, but with a better prognosis in osteosarcoma patients. Furthermore, the connection between DEGs and differentially expressed miRNAs (DEMs) was explored. GEO2R was utilized to screen out DEGs and DEMs. A total of 668 DEGs and 5 DEMs were extracted from GSE7437 and GSE30934 differentiating samples of poor and good chemotherapy reaction patients. The Database for Annotation, Visualization, and Integrated Discovery (DAVID) was used to perform GO and KEGG pathway enrichment analysis to identify potential pathways and functional annotations linked with osteosarcoma chemoresistance. The present study may provide a deeper understanding about regulatory genes of osteosarcoma chemoresistance and identify potential therapeutic targets for osteosarcoma.

1. Introduction

Osteosarcoma is one of the most common primary malignant bone tumors in children and adolescents. The worldwide morbidity rates of osteosarcoma are approximately with an average incidence of 3.1 per million for each stage and 4.4 per million for groups <25 years old. Additionally, there is a bimodal age distribution: individuals aged 25–60 years and elderly individuals, respectively. In America the Annual age-standardized incidence of osteosarcoma has reached stabilization from 1976 to 2005 [1–5].

Many factors are associated with tumor genesis, including high birth weight [6], pubertal hormones [7], and germline genetic variants [8, 9]. The common subtypes are osteoblastic, chondroblastic, and fibroblastic osteosarcomas, which may account for 70–80% of total cases [10]. The standard therapy consists of neoadjuvant chemotherapy (NACT), surgical removal of the primary tumor, and adjuvant chemotherapy. Before the 1970s, no more than 20 percent of patients were alive after 5 years when excision was major therapeutic measure for osteosarcoma [11, 12], while it increased to 60–70 percent for children and young adults with localized disease after the chemotherapy was used as adjuvant therapy for surgical resection [13, 14]. The current predicament of osteosarcoma treatment is the five-year survival rate does not exceed 25% for patients aged 2–68 with poor initial response tending to have adverse outcomes [15]. Therefore, to improve and modify chemotherapy regimens, an increasing number of pharmacogenomics studies on osteosarcoma have been going on for some time, such as drug reactions and toxicity.

Multidrug resistance protein 1 (MDR1) that is encoded by gene ATP-binding cassette, subfamily B (MDR/TAP), member 1 (ABCB1) has been shown to serve as a plausible factor in doxorubicin resistance, which was validated to be linked with poor outcomes in many osteosarcoma studies [16, 17], but whether there would be more valuable biomarkers remained to be explored. In recent years, microarray technology has substantially promoted the advance of understanding the mechanisms underlying diseases. Additionally, the rapid development of bioinformatics enables us to comprehensively screen out the hub genes associated with chemoresistance by the process of high-throughput microarrays. MicroRNAs (miRNAs), a group of highly conserved short noncoding small RNAs including generally 18–25 nucleotides in length, can suppress the translation of mRNA and cleave it by the modality of base-pairing to the target genes’ 3′ untranslated region [18–20].

In the present study, we analyzed the data of GSE87437 and GSE30934 submitted by Serra Mand and Kobayashi E. et al., respectively, to get 668 differentially expressed genes (DEGs) and 5 differentially expressed miRNA (DEMs) between samples of poor and good chemotherapy reaction patients in GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/). To further understand the function of genes, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes pathway (KEGG), Protein-Protein Interaction (PPI) networks, and the connections among DEGs and DEMs were performed in sequence. We selected chemoresistance development related key genes and provided theoretical foundations for modifying and improving osteosarcoma treatment methods.

2. Materials and Methods

2.1. Microarray Data

The dataset of GSE87437, gene expression array, and GSE30934, miRNA expression array, included 10 and 8 samples from poor chemotherapy reaction patients and 11 and 16 samples from good ones, respectively. Moreover, the former dataset was based on GPL570 platform ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array) and the latter one was based on GPL10312 platform (3D-Gene Human miRNA Oligo chip v12-1.00).

2.2. Identification of DEGs

GEO2R, an R-associated web application, was applied to filtrate DEGs between good chemotherapy reaction samples and poor chemotherapy reaction samples. In total, 21 samples in GSE87437 and 24 samples in GSE30934 were divided into two groups, respectively, and the concrete grouping schemes were already shown in microarray data. The and ⩾ 1 were considered as cutoff criterion. All results of DEGs were downloaded in text format, hierarchical clustering analysis being conducted later in Morpheus (https://software.broadinstitute.org/morpheus/).

2.3. GO and Pathway Enrichment Analysis of DEGs

The online tool, Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) provided comprehensive information for list of genes by GO and KEGG pathway analyses. In addition, GO enrichment analysis included three different aspects: biological process (BP), molecular function (MF), and cellular component (CC) [21]. KEGG enrichment analysis was associated with genomic information’s functional interpretation and practical application [22]. The screened DEGs were uploaded to DAVID V6.8 to perform GO and KEGG pathway analysis with the criterion of , the results of which were downloaded in text format.

2.4. PPI Networks Construction and Module Analysis

To analyze the connection among proteins, DEGs were uploaded to Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/), a database covering 9,643,763 proteins from 2,031 organisms, and the result whose minimum interaction score was 0.4 was visualized in Cytoscape [23, 24]. Furthermore, the Molecular Complex Detection (MCODE) was used to screen out significant modules based on the constructed PPI networks with the criteria of degree cutoff = 2, node density cutoff = 0.1, node score cutoff = 0.2, -core = 2, and max. depth = 100 and hub genes were exported. The functional enrichment analysis of genes in each module was performed in DAVID. Besides, the genes in each module were uploaded to DAVID and KEGG pathway enrichment analysis was conducted with the condition of .

2.5. Survival Analysis of Hub Genes

The series matrix of GSE21257 that contained osteosarcoma patients’ prognostic information was downloaded from GEO database. The patients were split into two groups, high expression and low expression, according to the expression level of a specific hub gene. The data was processed by graphpad prism software and then exported the results.

2.6. Prediction of miRNA Targets

DEMs were acquired by the parallel method of DEGs mentioned above. miRWalk1.0 (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk/index.html), an integrated miRNA target prediction platform including 10 databases (DIANAmT, miRanda, miRDB, miRWalk, RNAhybrid, PICTAR4, PICTAR5, PITA, RNA22, and Targetscan), was utilized to explore the correlation between DEMs and DEGs. Besides, different colors were used to indicate the degrees of connections. For example, red color represented strong correlation.

3. Results

3.1. Identification of DEGs

A total of 668 DEGs were obtained from GSE87437 in the poor chemotherapy response samples and compared with those of good response with the criteria of and ⩾ 1.0, comprising 422 upregulated genes and 246 downregulated genes. The key DEGs are displayed in Table 1.

Table 1: Key differentially expressed genes (DEGs) obtained from GSE87374.

Gene symbol

value

ZNRD1
1.4
0.01504282

CDK1
1.31
0.03504221

MYH7B
1.19
0.03527194

GPR68
1.18
0.03724459

CAT
−1.15
0.03367398

FUT3
1.69
0.00007205

IMPG2
1.11
0.04785566

GPR180
−1.03
0.00511768

ANPEP
1.32
0.0354653

3.2. Hierarchical Clustering Analysis of DEGs

Hierarchical clustering analysis was conducted through Morpheus, a web-based online tool, with the series matrix data of the DEGs. The heat map is shown in Figure 1 (top 50 upregulated and 50 downregulated genes).

Figure 1: Heat map of the top 100 DEGs (50 upregulated genes and 50 downregulated genes). Red: up-regulation; Blue: down-regulation.

3.3. GO Term Enrichment Analysis

In order to understand the function of the identified DEGs deeply, GO and KEGG analyses were performed in DAVID, respectively. The result of GO analysis showed that DEGs were enriched in biological process (BP), including positive regulation of transcription and DNA-templated, positive regulation of sequence-specific DNA binding transcription factor activity, nitric oxide mediated signal transduction, positive regulation of transcription from RNA polymerase II promoter, and regulation of phosphatidylinositol 3-kinase signaling. As for molecular function (MF), the DEGs were enriched in estrogen response element binding, Rac guanyl-nucleotide exchange factor activity, calcium ion binding, zinc ion binding, and phosphatidylinositol-4,5-bisphosphate 3-kinase activity. Besides, Cellular Component (CC) analysis showed that the DEGs were enriched in proteinaceous extracellular matrix, cell surface, P granule, integral component of plasma membrane, and endocytic vesicle membrane, as shown in Figure 2.

(a)

(b)

(c)

(d)

Figure 2: GO and KEGG pathway analysis of DEGs associated with osteosarcoma. (a) Top 5 significantly enriched biological processes in DEGs. (b) Top 5 significantly enriched cell component in DEGs. (c) Top 5 significantly enriched molecular function in DEGs. (d) Top 5 significantly enriched KEGG pathway in DEGs.

3.4. KEGG Pathway Analysis

KEGG pathway analysis showed that DEGs were mainly involved in tryptophan metabolism, oxytocin signaling pathway, glyoxylate and dicarboxylate metabolism, cAMP signaling pathway, and dopaminergic synapse (Figure 2).

3.5. PPI Networks and Modules Selection

The PPI networks of DEGs were composed of 432 nodes and 428 edges (Figure 3). Then the networks were imported into Cytoscape software, analyzed by using plug-ins MCODE. Eventually, 3 significant modules were selected (Figure 4), and the KEGG pathway was mainly associated with ribosome biogenesis in eukaryotes, calcium signaling pathway, arachidonic acid metabolism, proteoglycans in cancer, and linoleic acid metabolism (Figure 4).

Figure 3: PPI network constructed by STRING database for the DEGs.

Figure 4: The top 3 modules from the PPI network. (a) module 1, (b) module 2, (c) module 3, (d) the enriched pathways of module 1, (e) the enriched pathways of module 2, and (f) the enriched pathways of module 3.

3.6. Hub Genes and Survival Analysis

9 hub genes were screened out, including zinc ribbon domain containing 1 (ZNRD1), myosin heavy chain 7B (MYH7B), G protein-coupled receptor 68 (GPR68), catalase (CAT), fucosyltransferase 3 (Lewis blood group) (FUT3), interphotoreceptor matrix proteoglycan 2 (IMPG2), G protein-coupled receptor 180 (GPR180), alanyl aminopeptidase, membrane (ANPEP), and cyclin dependent kinase 1 (CDK1) (Table 1). Next, survival analysis of these genes in GSE21257 which contained patients’ survival prognostic information showed that osteosarcoma patients with high mRNA expression of FUT3 meant a better overall survival (OS) despite its high expression in poor chemotherapy response samples (Figure 5). Additionally, the survival prognostic information of GPR180 and CDK1 was not included in GSE21257.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 5: The survival prognostic value of hub gene in osteosarcoma from the GSE21257.

3.7. MiRNA–DEG Pairs

After the differentially expressed analysis for the data of GSE30934, a total of 5 DEMs were obtained between the poor chemotherapy response samples compared with that of good response with the criteria of and ⩾ 1.0 (Table 2). Next, basing on miRWalk1.0 database (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk/index.html), the relationship between miRNAs and DEGs was acquired and different kinds of colors were on behalf of the number of miRNA–DEG pairs in different database which stand for the degrees of connection. For example, red color represented to a strong correlation (Figure 6). After comparing the targets with hub genes, we found that ZNRD1 was the potential target of hsa-miR-543, while CAT was the potential target of hsa-miR-518f. Both hsa-miR-543 and hsa-miR-518f matched the regulated gene in expression trends.

Table 2: Key differentially expressed genes (DEGs) obtained from GSE30934.

miRNA_ID

value

hsa-miR-543
−3.429933
0.00192

hsa-miR-409-5p
−2.70517
0.00729

hsa-miR-518f
1.448711
0.02332

hsa-miR-154
−2.619116
0.03838

ebv-miR-BART1-3p
−1.358071
0.04733

Figure 6: The network of miRNA–DEG pairs.

4. Discussion

In the present study, we observed whether there were more valuable genes like ABCB1 which could help improve and modify chemotherapy regimens in osteosarcoma. To find out the specific chemotherapy response-associated DEGs, we analyzed the osteosarcoma gene expression array of GSE87437 in GEO2R, where a total number of 668 DEGs were obtained between good and poor chemotherapy response samples. Besides, to further understand the potential biological functions, we conducted GO, KEGG, and STRING analyses. Subsequently, on the foundation of PPI networks, the selection of 9 hub genes and their survival prognosis were completed. In terms of the increasingly prominent role of miRNA in cancer, DEMs of osteosarcoma miRNA expression array of GSE30934 was screened out in the same way and criteria like DEGs, and DEGs-miRNA network was constructed to show relationship between them [25].

Our results showed that many genes and miRNAs may have functions in the development of chemoresistance in osteosarcoma and have the potential to become treatment targets. Here, we exclusively focused on 9 hub genes and two miRNAs. Firstly, 9 hub genes consisted of ZNRD1, MYH7B, GPR68, CAT, FUT3, IMPG2, GPR180, ANPEP, and CDK1. Our data showed that the expression of ZNRD1 was upregulated in chemoresistance osteosarcoma samples. Previous studies demonstrated that methotrexate-resistant, vincristine-resistant, multidrug resistant phenotypes of gastric cancer cells could be regulated by the inhibition of ZNRD1/Inosine monophosphate dehydrogenase 2 (IMPDH2), upregulated DARPP-32/downregulated ZNRD1, overexpressed miR-508-5p/ZNRD1/ABCB1 activities, respectively [26–28], but further researches of ZNRD1 in osteosarcoma chemoresistance remained to be conducted. Similar to ZNRD1, MYH7B was also found upregulated in chemoresistance osteosarcoma samples. At present, MYH7B was mainly involved in pathway of cardiomyocytes, such as mitochondrial apoptosis pathway [29] but studies about cancer were rare. GPR68, a kind of pH-sensing protein, was associated with tumor cell biology, such as tumor aggressiveness by triggering the intracellular signaling cascade to promote the development of microenvironment of extracellular acidification [30, 31]. As shown in Table 1, likewise, we found the expression of GPR68 was upregulated in chemoresistance samples. Daglioglu C validated that pH-responsive Fe3O4@SiO2(FITC)-BTN/QUR/DOX multifunctional nanoparticles could potentiate the chemotherapeutic efficacy of DOX against multidrug resistance as well as counteract the survival ability of chemoresistant lung carcinoma A549/DOX cell lines [32]. GPR68 has become an attractive target for drug development [33]. Several previous studies demonstrated that decreased CAT was highly associated with chemoresistance; for example, Xu et al. showed that intervention against miR-551b/CAT/reactive oxygen species (ROS)/Mucin-1 (MUC1) pathway might help overcome acquired chemoresistance [34]. Tumor microenvironment (TME) was characterized by hypoxia, acidosis, and dense extracellular matrix, providing tumors with resistance to various therapies, which could be effectively changed by the intravenous injection of human serum albumin (HAS)-chlorine e6 (Ce6)-CAT-paclitaxel (PTX) nanoparticles, enzyme-loaded therapeutic albumin nanoparticles. Meanwhile, H2O2 could relieve tumor hypoxia by generating oxygen within TME triggered by CAT of those nanoparticles, which made CAT a potential treatment target in various tumors [35]. Similarly, our data showed that CAT was downregulated in chemoresistance samples, which might exacerbate local microenvironment to strengthen tumor chemoresistance through the way of hypoxia, subsequent acidosis, and the like. High expression of FUT3 was proved to participate in the development of invasion, metastasis, and resistance to therapy by increased fucosylation activity in oral squamous cell carcinoma (OSCC), and the function could be blocked by inhibition of fucosylation [36]. In our study, FUT3 was also observed upregulated in chemoresistance samples, but it was found to be associated with a better survival prognosis (Figure 4). There were some reasons to explain that. In spite of its protumor action, death pathways were proved to be relied on fucosylation, and FUT3 was demonstrated to play an important role in natural killer-induced cytotoxicity after the recognition of sialyl Lewis X with the help of C-type lectin receptors [37, 38]. Therefore, the relationship between FUT3 and tumor was so complex that little was known about its function in osteosarcoma chemoresistance. IMPG-2, a gene mainly associated with retinal disease, was upregulated in chemoresistant osteosarcoma samples in our study, but more studies involving IMPG-2 and cancer needed to be conducted [39]. Furthermore, our data showed the decreased expression of GPR180 in poor chemotherapy response samples. Homoplastically, Honda et al. found that the methylation of GPR180 was probably to encode tumor suppressors and serves as a novel prognostic marker and therapeutic target in Hepatoblastoma [40]. Based on previous studies, whether the gene GPR180 could have function in osteosarcoma by producing tumor suppressors and its concrete role in chemoresistance remained to be explored. The gene ANPEP encoded aminopeptidase N (APN). A previous study [41] showed that ANPEP was downregulated in prostate cancer (PC). On the contrary, our study showed that ANPEP expression of good chemotherapy response samples was approximately two times that in chemoresistance osteosarcoma samples. However, the difference caused by the types of tumors or chemoresistance needed to be further studied. In urothelial carcinoma, APN could increase cytotoxicity of melphalan-flufenamide to play anticancer effect by amplifying the intracellular loading of melphalan [42]. The studies of chemoresistance associated with ANPEP have not been conducted so far, but Viktorsson et al. [42] offered researchers a new way for treatment, which made ANPEP a significant therapeutic target. Several researches in chemoresistance-associated fields have already demonstrated that CDK1 participated in the development of chemoresistance in pathways, such as GRP78/CDC2 [43]. Likewise, in our study, the expression of CDK1 was increased. Hayashi et al. found increased DNA repair activity in the G2-M transition promoted temozolomide (TMZ) resistance and CDC2 inhibitor flavopiridol (FP) treatment could resensitize TMZ-resistant clones in a p53-independent manner in glioma cells [44]. Besides, the combination of ERK inhibitor PD98059 and Taxol could improve the sensitivity of taxol-resistant tumor cells with the decreased CDC2 activity [45].

Compared with that of good response, 5 DEMs were acquired in GSE30934 in poor chemotherapy response samples (Table 2). Among them, hsa-miR-543 and hsa-miR-518f were found to have a relation to ZNRD1 and CAT, respectively. In our results, hsa-miR-543 was downregulated in chemoresistant samples. Previous study in this field was limited. In other aspects, decreased expression of it was involved in osteosarcoma angiogenesis which might be caused by connective tissue growth factor (CTGF) in phospholipase C (PLC)/protein kinase C (PKCδ) signaling pathway. Besides, hsa-miR-543 was also proved to be linked with tumor staging [46], cell proliferation, and the glycolytic pathway [47]. The studies between chemoresistance-promoted gene ZNRD1 and hsa-miR-543 have not been conducted yet, but biological functions mentioned above made hsa-miR-543 become an important therapeutic target. Moreover, the role of hsa-miR-518f in chemoresistance or the development of tumor was rarely known. Hsa-miR-518e and hsa-miR-518b, homologues of hsa-miR-518f, were demonstrated to be upregulated in hepatocellular carcinoma (HCC) [48]. Besides, a previous study showed that has-miR-518c-5p could regulate the growth and metastasis of oral cancer [49]. Consequently, further research to hsa-miR-518f was of great importance.

In summary, we identified 668 DEGs and 5 DEMs from GEO2R between good chemotherapy response samples and poor chemoresistance samples in osteosarcoma. And many of them, such as ZNRD1, GPR68, CAT, FUT3, ANPEP, CDK1, and hsa-miR-543, might be key genes related to osteosarcoma chemoresistance. These findings provided a series of promising treatment targets and enlightened us on the further investigations of the molecular mechanisms.

Conflicts of Interest

The authors declare that there are no conflicts of interest.

References

L. Mirabello, R. J. Troisi, and S. A. Savage, “International osteosarcoma incidence patterns in children and adolescents, middle ages and elderly persons,” International Journal of Cancer, vol. 125, no. 1, pp. 229–234, 2009.

L. Mirabello, R. J. Troisi, and S. A. Savage, “Osteosarcoma incidence and survival rates from 1973 to 2004: data from the surveillance, epidemiology, and end results program,” Cancer, vol. 115, no. 7, pp. 1531–1543, 2009.

H. D. Dorfman and B. Czerniak, “Bone cancers,” Cancer, vol. 75, no. 1, supplement, pp. 203–210, 1995.

K. P. Anfinsen, S. S. Devesa, F. Bray et al., “Age-period-cohort analysis of primary bone cancer incidence rates in the United States (1976-2005),” Cancer Epidemiology, Biomarkers & Prevention, vol. 20, no. 8, pp. 1770–1777, 2011.

P. C. Valery, M. Laversanne, and F. Bray, “Bone cancer incidence by morphological subtype: a global assessment,” Cancer Causes & Control, vol. 26, no. 8, pp. 1127–1139, 2015.

L. Mirabello, R. Pfeiffer, G. Murphy et al., “Height at diagnosis and birth-weight as risk factors for osteosarcoma,” Cancer Causes & Control, vol. 22, no. 6, pp. 899–908, 2011.

J. R. B. Musselman, T. L. Bergemann, J. A. Ross et al., “Case-parent analysis of variation in pubertal hormone genes and pediatric osteosarcoma: a children's oncology group (COG) study,” International Journal of Molecular Epidemiology and Genetics, vol. 3, no. 4, pp. 286–293, 2012.

M. Kansara, M. W. Teng, M. J. Smyth, and D. M. Thomas, “Translational biology of osteosarcoma,” Nature Reviews Cancer, vol. 14, no. 11, pp. 722–735, 2014.

L. Mirabello, K. Yu, S. I. Berndt et al., “A comprehensive candidate gene approach identifies genetic variation associated with osteosarcoma,” BMC Cancer, vol. 11, Article ID 209, 2011.

S. H. Orkin, D. E. Fisher, A. T. Look, S. Lux, D. Ginsburg, and D. G. Nathan, Oncology of Infancy and Childhood, Elsevier Health Sciences, 2009.

D. C. Dahlin and M. B. Coventry, “Osteogenic sarcoma, a study of six hundred cases,” The Journal of Bone & Joint Surgery, vol. 49, no. 1, pp. 101–110, 1967.

R. C. Marcove, V. Miké, J. V. Hajek, A. G. Levin, and R. V. Hutter, “Osteogenic sarcoma under the age of twenty-one,” The Journal of Bone & Joint Surgery, vol. 52, no. 3, pp. 411–423, 1970.

P. A. Meyers, G. Heller, J. Healey et al., “Chemotherapy for nonmetastatic osteogenic sarcoma: the memorial sloan-kettering experience,” Journal of Clinical Oncology, vol. 10, no. 1, pp. 5–15, 1992.

A. M. Goorin, D. J. Schwartzentruber, M. Devidas et al., “Presurgical chemotherapy compared with immediate surgery and adjuvant chemotherapy for nonmetastatic osteosarcoma: Pediatric Oncology Group Study POG-8651,” Journal of Clinical Oncology, vol. 21, no. 8, pp. 1574–1580, 2003.

B. Kempf-Bielack, S. S. Bielack, H. Jürgens et al., “Osteosarcoma relapse after combined modality therapy: an analysis of unselected patients in the Cooperative Osteosarcoma Study Group (COSS),” Journal of Clinical Oncology, vol. 23, no. 3, pp. 559–568, 2005.

H. I. Vos, M. J. H. Coenen, H.-J. Guchelaar, and D. M. W. M. te Loo, “The role of pharmacogenetics in the treatment of osteosarcoma,” Drug Discovery Therapy, vol. 21, no. 11, pp. 1775–1786, 2016.

M. Serra and C. M. Hattinger, “The pharmacogenomics of osteosarcoma,” The Pharmacogenomics Journal, vol. 17, no. 1, pp. 11–20, 2017.

V. N. Kim, J. Han, and M. C. Siomi, “Biogenesis of small RNAs in animals,” Nature Reviews Molecular Cell Biology, vol. 10, no. 2, pp. 126–139, 2009.

D. P. Bartel, “MicroRNAs: target recognition and regulatory functions,” Cell, vol. 136, no. 2, pp. 215–233, 2009.

M. A. Valencia-Sanchez, J. Liu, G. J. Hannon, and R. Parker, “Control of translation and mRNA degradation by miRNAs and siRNAs,” Genes & Development, vol. 20, no. 5, pp. 515–524, 2006.

P. Gaudet, N. Škunca, J. C. Hu, and C. Dessimoz, “Primer on the Gene Ontology,” in The Gene Ontology Handbook, vol. 1446 of Methods in Molecular Biology, pp. 25–37, Springer New York, New York, NY, USA, 2017.

M. Kanehisa, M. Furumichi, M. Tanabe, Y. Sato, and K. Morishima, “KEGG: new perspectives on genomes, pathways, diseases and drugs,” Nucleic Acids Research, vol. 45, no. 1, pp. D353–D361, 2017.

D. Szklarczyk, A. Franceschini, S. Wyder et al., “STRING v10: protein-protein interaction networks, integrated over the tree of life,” Nucleic Acids Research, vol. 43, pp. D447–D452, 2015.

P. Shannon, A. Markiel, O. Ozier et al., “Cytoscape: a software environment for integrated models of biomolecular interaction networks,” Genome Research, vol. 13, no. 11, pp. 2498–2504, 2003.

G. A. Calin, C. Sevignani, C. D. Dumitru et al., “Human microRNA genes are frequently located at fragile sites and genomic regions involved in cancers,” Proceedings of the National Acadamy of Sciences of the United States of America, vol. 101, no. 9, pp. 2999–3004, 2004.

L. Hong, T. Qiao, Y. Han et al., “ZNRD1 mediates resistance of gastric cancer cells to methotrexate by regulation of IMPDH2 and Bcl-2,” The International Journal of Biochemistry & Cell Biology, vol. 84, no. 2, pp. 199–206, 2006.

L. Hong, Y. Zhao, J. Wang et al., “Reversal of multidrug resistance of adriamycin-resistant gastric adenocarcinoma cells through the up-regulation of DARPP-32,” Digestive Diseases and Sciences, vol. 53, no. 1, pp. 101–107, 2008.

Y. Shang, B. Feng, L. Zhou et al., “The miR27b-CCNG1-P53-miR-508-5p axis regulates multidrug resistance of gastric cancer,” Oncotarget, vol. 7, no. 1, pp. 538–549, 2016.

J. Wang, Z. Jia, C. Zhang et al., “miR-499 protects cardiomyocytes from H₂O₂-induced apoptosis via its effects on Pdcd4 and Pacs2,” RNA Biology, vol. 11, no. 4, pp. 339–350, 2014.

W.-C. Huang, P. Swietach, R. D. Vaughan-Jones, O. Ansorge, and M. D. Glitsch, “Extracellular acidification elicits spatially and temporally distinct Ca2+ signals,” Current Biology, vol. 18, no. 10, pp. 781–785, 2008.

H. Saxena, D. A. Deshpande, B. C. Tiegs et al., “The GPCR OGR1 (GPR68) mediates diverse signalling and contraction of airway smooth muscle in response to small reductions in extracellular pH,” British Journal of Pharmacology, vol. 166, no. 3, pp. 981–990, 2012.

C. Daglioglu, “Enhancing tumor cell response to multidrug resistance with pH-sensitive quercetin and doxorubicin conjugated multifunctional nanoparticles,” Colloids and Surfaces B: Biointerfaces, vol. 156, pp. 175–185, 2017.

X.-P. Huang, J. Karpiak, W. K. Kroeze et al., “Allosteric ligands for the pharmacologically dark receptors GPR68 and GPR65,” Nature, vol. 527, no. 7579, pp. 477–483, 2015.

X. Xu, A. Wells, M. T. Padilla, K. Kato, K. C. H. Kim, and Y. Lin, “A signaling pathway consisting of miR-551b, catalase and MUC1 contributes to acquired apoptosis resistance and chemoresistance,” Carcinogenesis, vol. 35, no. 11, pp. 2457–2466, 2014.

Q. Chen, J. Chen, C. Liang et al., “Drug-induced co-assembly of albumin/catalase as smart nano-theranostics for deep intra-tumoral penetration, hypoxia relieve, and synergistic combination therapy,” Journal of Controlled Release, vol. 263, pp. 79–89, 2017.

V. Desiderio, P. Papagerakis, V. Tirino et al., “Increased fucosylation has a pivotal role in invasive and metastatic properties of head and neck cancer stem cells,” Oncotarget, vol. 6, no. 1, pp. 71–84, 2015.

K. Higai, A. Ichikawa, and K. Matsumoto, “Binding of sialyl Lewis X antigen to lectin-like receptors on NK cells induces cytotoxicity and tyrosine phosphorylation of a 17-kDa protein,” Biochimica et Biophysica Acta (BBA)—General Subjects, vol. 1760, no. 9, pp. 1355–1363, 2006.

C. Ohyama, S. Kanto, K. Kato et al., “Natural killer cells attack tumor cells expressing high levels of sialyl Lewis x oligosaccharides,” Proceedings of the National Acadamy of Sciences of the United States of America, vol. 99, no. 21, pp. 13789–13794, 2002.

D. Bandah-Rozenfeld, R. W. J. Collin, E. Banin et al., “Mutations in IMPG2, Encoding interphotoreceptor matrix proteoglycan 2, cause autosomal-recessive retinitis pigmentosa,” American Journal of Human Genetics, vol. 87, no. 2, pp. 199–208, 2010.

S. Honda, M. Minato, H. Suzuki et al., “Clinical prognostic value of DNA methylation in hepatoblastoma: four novel tumor suppressor candidates,” Cancer Science, vol. 107, no. 6, pp. 812–819, 2016.

K. D. Sørensen, M. O. Abildgaard, C. Haldrup et al., “Prognostic significance of aberrantly silenced ANPEP expression in prostate cancer,” British Journal of Cancer, vol. 108, no. 2, pp. 420–428, 2013.

K. Viktorsson, C.-H. Shah, T. Juntti et al., “Melphalan-flufenamide is cytotoxic and potentiates treatment with chemotherapy and the Src inhibitor dasatinib in urothelial carcinoma,” Molecular Oncology, vol. 10, no. 5, pp. 719–734, 2016.

W. Li, W. Wang, H. Dong et al., “Cisplatin-induced senescence in ovarian cancer cells is mediated by GRP78,” Oncology Reports, vol. 31, no. 6, pp. 2525–2534, 2014.

T. Hayashi, K. Adachi, S. Ohba, and Y. Hirose, “The Cdk inhibitor flavopiridol enhances temozolomide-induced cytotoxicity in human glioma cells,” Journal of Neuro-Oncology, vol. 115, no. 2, pp. 169–178, 2013.

X. Lin, Y. Liao, X. Chen, D. Long, T. Yu, and F. Shen, “Regulation of oncoprotein 18/Stathmin signaling by ERK concerns the resistance to Taxol in nonsmall cell lung cancer cells,” Cancer Biotherapy and Radiopharmaceuticals, vol. 31, no. 2, pp. 37–43, 2016.

L.-H. Wang, H.-C. Tsai, Y.-C. Cheng et al., “CTGF promotes osteosarcoma angiogenesis by regulating miR-543/angiopoietin 2 signaling,” Cancer Letters, vol. 391, pp. 28–37, 2017.

H. Zhang, X. Feng, T. Wang et al., “MiRNA-543 promotes osteosarcoma cell proliferation and glycolysis by partially suppressing PRMT9 and stabilizing HIF-1α protein,” Oncotarget, vol. 8, no. 2, pp. 2342–2355, 2017.

W. Wang, L.-J. Zhao, Y.-X. Tan, H. Ren, and Z.-X. Qi, “MiR-138 induces cell cycle arrest by targeting cyclin D3 in hepatocellular carcinoma,” Carcinogenesis, vol. 33, no. 5, pp. 1113–1120, 2012.

M. Kinouchi, D. Uchida, N. Kuribayashi et al., “Involvement of miR-518c-5p to growth and metastasis in oral cancer,” PLoS ONE, vol. 9, no. 12, Article ID e115936, 2014.