Identification of Potential Therapeutic Targets and Immune Cell Infiltration Characteristics in Osteosarcoma Using Bioinformatics Strategy

Osteosarcoma is one of the most aggressive malignant bone tumors worldwide. Although great advancements have been made in its treatment owing to the advent of neoadjuvant chemotherapy, the problem of lung metastasis is a major obstacle in the improvement of survival outcomes. Thus, the aim of the present study is to screen novel and key biomarkers, which may act as potential prognostic markers and therapeutic targets in osteosarcoma. We utilized the robust rank aggregation (RRA) method to integrate three osteosarcoma microarray datasets downloaded from the Gene Expression Omnibus (GEO) database, and we identified the robust differentially expressed genes (DEGs) between primary and metastatic osteosarcoma tissues. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to explore the functions of robust DEGs. The results of enrichment analysis showed that the robust DEGs were closely associated with osteosarcoma development and progression. Immune cell infiltration analysis was also conducted by CIBERSORT algorithm, and we found that macrophages are the most principal infiltrating immune cells in osteosarcoma, especially macrophages M0 and M2. Then, the protein–protein interaction network and key modules were constructed by Cytoscape, and 10 hub genes were selected by plugin cytoHubba from the whole network. The survival analysis of hub genes was also carried out based on the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database. The integrated bioinformatics analysis was utilized to provide new insight into osteosarcoma development and metastasis and identified EGR1, CXCL10, MYC, and CXCR4 as potential biomarkers for prognosis of osteosarcoma.

Introduction

Osteosarcoma is one of the most aggressive malignant tumors in the bone (1), which derived from mesenchymal tissue and showed osteoblastic differentiation (2). Annually, the incidence rate of osteosarcoma is approximately four to five cases per million (3). In addition, osteosarcoma is most prevalent in children and adolescents (4), and 15–20% of osteosarcoma patients have lung metastasis at the initial diagnosis (5, 6). With the assistance of neoadjuvant chemotherapy, the treatment of osteosarcoma has been greatly improved, but the overall survival of patients with lung metastasis or relapse has not improved and remains low at approximately 20% (7, 8). Therefore, it is extremely necessary to seek novel prognostic factors and therapeutic targets for osteosarcoma.

In recent years, public databases including Gene Expression Omnibus (GEO), The Cancer Genome Atlas (TCGA), and Therapeutically Applicable Research to Generate Effective Treatments (TARGET) are widely used to explore diagnostic and prognostic biomarkers in osteosarcoma. In the previous studies, the limited number of samples and inappropriate analysis methods of multiple datasets led to deviation of the results. Zhang et al. (9) and Diao et al. (10) described the role of gene copy number alterations and methylation changes in the malignant progression of osteosarcoma using bioinformatics analyses, respectively. However, these studies were only based on single datasets and had a limited sample size, which may have biased the final results. To analyze more samples and avoid the sample heterogeneity of each independent experiment and the error caused by different technology platforms and different data processing methods, we used the robust rank aggregation (RRA) method to obtain robust differentially expressed genes (DEGs). RRA was used to compare the ranking of multiple gene lists. If a gene ranked the highest in all gene lists, then the smaller its calculated P-value is, the more likely it is to be a robust DEG (11). This method has been widely used in integrated analysis of multiple datasets, and it is robust to errors and noise (12–14). There have been no reports of the use of RRA in osteosarcoma.

In our study, microarray datasets GSE14827, GSE21257, and GSE32981 from GEO database were downloaded and analyzed by RRA method to identify robust DEGs between primary and metastatic osteosarcoma tissues. A total of 524 robust DEGs were determined, including 272 upregulated and 252 downregulated genes. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to explore the functions of robust DEGs. Immune cell infiltration analysis was also conducted by CIBERSORT algorithm. In addition, we also constructed the protein–protein interaction (PPI) network and key modules, and finally 10 hub genes were selected by plugin cytoHubba from the whole network. Survival analysis of hub genes was carried out using R packages. In conclusion, the integrated bioinformatics analysis was utilized to identify the significant robust DEGs and hub genes, which may act as novel and potential prognostic biomarkers in osteosarcoma.

Materials and Methods

Data Collection and Data Processing

We selected three gene chips of osteosarcoma from GEO database¹, including GSE14827, GSE21257, and GSE32981. The selection criteria were as follows: (1) inclusion of primary and metastatic osteosarcoma tissue samples; (2) expression profiling by array as the experiment type; (3) Homo sapiens; and (4) 20 samples as the minimal size. Among them, GSE14827 contained 18 primary samples and 9 metastatic osteosarcoma samples. GSE21257 contained 19 non-metastatic and 34 metastatic osteosarcoma samples. The GSE32981 dataset contained 23 samples. One sample was not available and was excluded. Twenty-two samples were selected for further study including four primary and 18 metastatic samples. The matrix files and platform annotation document of three microarray datasets were downloaded. The names of microarray probes were converted to the gene symbols by Perl. The DEGs were identified between primary and metastatic osteosarcoma samples in each dataset by limma package in R (15) with the cutoff criteria of |log2 fold change (FC)| > 0.585 and P-value < 0.05.

Robust Rank Aggregation Analysis

To integrate the three microarray datasets, RRA method was used to determine the robust DEGs (16), which is a standard method to minimize the bias and errors among several datasets. Before RRA analysis, the upregulated and downregulated genes were ranked by their FC in each dataset. Then, the RobustRankAggreg R package was performed to get robust DEGs on the basis of the ranked genes in the three datasets. Genes with FC > 1.3 and P-value < 0.05 were considered as the significant robust DEGs.

Gene Ontology Function and Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analyses

To identify the functional roles of the robust DEGs indicated above, GO enrichment results of biological process (BP), cellular component (CC), and molecular function (MF) were obtained using the R package “clusterprofiler.” The KEGG pathway analysis of robust DEGs was also conducted using the R package (17). P < 0.05 was considered statistically significant.

Immune Infiltration by CIBERSORT Analysis

The CIBERSORT algorithm is commonly used to predict the infiltration of 22 types of immune cells in each tissue sample (18). The 22 kinds of immune cells include seven types of T cells [CD8⁺ T cells, naïve CD4⁺ T cells, resting memory CD4⁺ T cells, activated memory CD4⁺ T cells, follicular helper T cells, regulatory T cells (Tregs), gamma delta T cells], three types of macrophages (M0, M1, and M2), naïve B cells, memory B cells, plasma cells, resting natural killer (NK) cells, activated NK cells, monocytes, resting dendritic cells, activated dendritic cells, resting mast cells, activated mast cells, eosinophils, and neutrophils. Normalized gene expression matrix was converted to 22 types of immune cell matrix by the CIBERSORT algorithm. The immune cell matrix was filtered according to the criteria of P < 0.05, and then the relative expression of 22 types of immune cells was identified between primary and metastatic osteosarcoma samples by R packages. The principal component analysis (PCA) was also performed to determine the difference between primary and metastatic samples.

Protein–Protein Interaction Network Construction and Module Analysis

We uploaded the robust DEGs to the STRING online database², and we chose confidence >0.9 as the screening criteria. The visualized PPI network was performed by Cytoscape (version 3.6.1) software³. Cytoscape plugin-MCODE was used to screen the significant modules in the PPI network.

Hub Gene Identification

cytoHubba, a plugin of Cytoscape, provides several topological analysis algorithms, including Degree, Edge Percolated Component (EPC), Maximum Neighborhood Component (MNC), Density of Maximum Neighborhood Component (DMNC), Maximal Clique Centrality (MCC), and six centralities that include BottleNeck, EcCentricity, Closeness, Radiality, Betweenness, and Stress. These algorithms can be used to identify hub genes (19).

Survival Analysis

The RNA-seq-FPKM data and prognostic information of osteosarcoma patients were downloaded from TARGET database⁴. TARGET is a database that only includes children’s tumors. Presently, TARGET database contains six kinds of tumors, including ALL (Acute Lymphoblastic Leukemia), AML (Acute Myeloid Leukemia), KT (Kidney Tumors), MDLS (Model Systems), NBL (Neuroblastoma), and OS (Osteosarcoma). The survival analyses of hub genes were conducted by R package survival and survminer. P < 0.05 was considered to be statistically significant.

Results

Identification of Differentially Expressed Genes in Each Dataset

In the present study, the biological characteristics of DEGs were identified by integrated bioinformatics analysis. The overall workflow of this study is showed in Figure 1. The osteosarcoma microarray data GSE14827, GSE21257, and GSE32981 were selected and analyzed using limma package in R. A total of 102 osteosarcoma samples including 41 primary and 61 metastatic tissues were identified in our study. According to the cutoff criteria of | log2 FC| > 0.585 and P < 0.05, there were 83 DEGs in GSE14827 including 23 upregulated and 60 downregulated genes. A total of 464 DEGs were screened from the GSE21257 dataset, including 247 upregulated and 217 downregulated genes. Additionally, a total of 650 DEGs were selected in GSE32981, including 299 upregulated and 351 downregulated genes. The distribution of DEGs is shown in volcano plots (Figures 2A–C), where the red and green dots represent the upregulated and downregulated genes, respectively.

FIGURE 1

Figure 1. Study workflow. GEO, Gene Expression Omnibus; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; RRA, robust rank aggregation; TARGET, Therapeutically Applicable Research to Generate Effective Treatments.

FIGURE 2

Figure 2. Identification of DEGs and robust DEGs. Volcano plots of the distribution of DEGs in GSE14827 (A), GSE21257 (B), and GSE32981 (C). Red and green dots represent the upregulated and downregulated genes, respectively. (D) The heatmap of top 20 upregulated and downregulated robust DEGs identified by RRA method. Red represents high expression robust DEGs, while blue represents low expression robust DEGs. DEG, differentially expressed gene; RRA, robust rank aggregation.

Identification of Robust Differentially Expressed Genes by Robust Rank Aggregation Method

To integrate the three datasets with minimal bias, the RRA method was used. A total of 524 robust DEGs were determined, including 272 upregulated and 252 downregulated genes (Supplementary Table S1). According to the P-value of robust DEGs, we assigned the top 20 upregulated and downregulated robust DEGs in the visualized heatmap (Figure 2D).

Functional Enrichment Analyses of Robust Differentially Expressed Genes

To explore the functions of robust DEGs, the GO and KEGG enrichment analyses were conducted by R packages. The results of GO analysis included three categories: BP, CC, and MF. For BP, the upregulated robust DEGs were mainly enriched in embryonic organ development, multicellular organismal homeostasis, and transmembrane receptor protein serine/threonine kinase signaling pathway. In the CC part, the upregulated genes were particularly enriched in lamellar body, cell–cell junction, and cell–cell adherens junction. The top three significantly enriched terms were DNA-binding transcription activator activity, growth factor receptor binding, and cell adhesion molecule binding in the MF group (Figure 3A and Supplementary Table S2). Moreover, the most significantly enriched GO BP terms of downregulated genes were extracellular matrix (ECM) organization, extracellular structure organization, and regulation of cellular response to growth factor stimulus. For CC category, the downregulated genes were enriched in collagen-containing ECM. In addition, the downregulated DEGs were mainly enriched in ECM structural constituent conferring compression resistance, ECM structural constituent, and glycosaminoglycan binding in the MF group (Figure 3B and Supplementary Table S3). The above results indicated that the robust DEGs were mostly associated with cancer-related functions.

FIGURE 3

Figure 3. Functional enrichment analysis of robust DEGs. GO enrichment analyses of upregulated robust DEGs (A) and downregulated robust DEGs (B) in three parts: BP, CC, and MF. (C) KEGG pathway enrichment analysis of robust DEGs. DEG, differentially expressed gene; GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

The result of KEGG pathway enrichment analysis is also shown in Figure 3C. TGF-beta signaling pathway, wnt signaling pathway, and IL-17 signaling pathway were highly associated with tumor progression.

Immune Cell Infiltration Analysis

With the use of CIBERSORT algorithm, the infiltration of 22 kinds of immune cells in 102 osteosarcoma tissues is shown in Figure 4A. There was no significant difference in the infiltration of immune cells between primary and metastatic osteosarcoma tissues. However, compared with other immune cells, such as T cells and B cells, macrophage infiltration dominated, whether in primary or metastatic osteosarcoma tissues (Figure 4B). The above results demonstrated that macrophages may play an important role in the development and progression of osteosarcoma. The visualized violin plot was also constructed to prove the above findings (Figure 4C). The PCA of Figure 4D showed nothing individual difference between primary and metastasis samples.

FIGURE 4

Figure 4. Immune cells infiltration analysis. (A) The distribution of 22 types of immune cells between primary and metastatic osteosarcoma tissues. (B) The difference of immune cells infiltration between primary and metastatic osteosarcoma tissues visualized by heatmap. (C) Violin plot visualizing the differentially infiltrated immune cells (P < 0.05). (D) PCA performed on all osteosarcoma tissues. The two principal components showed nothing significant variation. PCA, principal component analysis.

Protein–Protein Interaction Network Construction and Module Analysis

To further study the interaction of robust DEGs, we constructed the PPI network by STRING database. With the confidence >0.9 and hiding the disconnected nodes, a visualized PPI network was created by Cytoscape (Figure 5A). In the final network, there were 148 nodes and 302 edges, including 84 upregulated and 64 downregulated genes. By using MCODE plugin, three key modules were screened from the whole network (Figures 5B–D). The robust DEGs in module 1 were mainly enriched in type I interferon signaling pathway. BP of genes in module 2 was particularly enriched in chemokine-mediated signaling pathway. In addition, genes in module 3 were mainly enriched in cell–cell adhesion mediated by cadherin, cell–cell junction assembly, and adherens junction organization (Supplementary Table S4).

FIGURE 5

Figure 5. Construction of PPI network, analysis of key modules, and identification of hub genes. (A) The whole PPI network. Upregulated genes are marked in red, while the downregulated genes are marked in green. (B) PPI network of module 1. (C) PPI network of module 2. (D) PPI network of module 3. (E) Hub genes were identified by intersection of 50 genes from 10 algorithms including MCC, DMNC, MNC, Degree, EPC, BottleNeck, EcCentricity, Closeness, Radiality, and Betweenness. PPI, protein–protein interaction.

Hub Gene Identification

cytoHubba is a Cytoscape plugin that allows the use several topological analysis algorithms, including MCC, DMNC, MNC, Degree, EPC, BottleNeck, EcCentricity, Closeness, Radiality, Betweenness, and Stress. These approaches can be used to predict and explore important nodes in PPI networks. Scores from topological algorithms are assigned to each node in a PPI network. According to the gene score, the top ranked genes can be considered as the hub genes. In the present study, we used 10 topological analysis algorithms (MCC, DMNC, MNC, Degree, EPC, BottleNeck, EcCentricity, Closeness, Radiality, and Betweenness) to rank the top 50 genes of the whole network. The intersection of these 50 genes from the 10 algorithms revealed the 10 hub genes: POMC, EGR1, CXCL10, SERPINA1, OAS1, MYC, CXCR4, CXCL2, CHRDL1, and GNAI1 (Figure 5E). The description of the 10 hub genes is shown in Table 1, including full names, synonyms and primary functions.

TABLE 1

Table 1. Description of the 10 hub genes.

Survival Analysis

Association between 10 hub genes and the overall survival of osteosarcoma patient were analyzed using R package. Based on each hub gene’s best-separation cutoff value, osteosarcoma patients’ samples within the TARGET-osteosarcoma dataset were divided into two groups to get the Kaplan–Meier (K-M) survival curves. The results demonstrated that gene changes of CXCL10 (P = 0.044), GNAI1 (P = 0.048), MYC (P = 0.011), and OAS1 (P = 0.0091) were significantly correlated with the overall survival of osteosarcoma patients (Figure 6).

FIGURE 6

Figure 6. Survival analysis. Gene changes of CXCL10 (A), GNAI1 (B), MYC (C), and OAS1 (D) were significantly correlated with the overall survival of osteosarcoma patients (P < 0.05).

Discussion

More and more research based on public database such as GEO and TARGET database has been done to determine biomarkers in osteosarcoma. For example, Wang et al. (20) used the microarray data of 42 different age groups (<20- and >20-year-old) osteosarcoma samples from GSE39058 to find 2113 DEGs, including 1476 upregulated and 637 downregulated genes. Similarly, they also identified 15 differentially expressed miRNAs (DEMs) in GSE39040, and functional enrichment analysis showed that upregulated DEMs were mainly enriched in cell growth and response to growth factor, and downregulated DEMs were involved in cytokine receptor activity. Moreover, using GEO database, Dai et al. (21) screened candidate genes for predicting the response to chemoresistance in osteosarcoma by miRNA–mRNA network. However, the differentially selected genes in these studies are all based on a single dataset, and the small size of samples will cause the instability of results. We integrated three datasets using RRA method, which is standard and robust, compared with other studies on osteosarcoma.

In our study, a total of 524 robust DEGs were determined by RRA method, including 272 upregulated and 252 downregulated genes. The results of GO and KEGG pathway enrichment analyses indicated that the robust DEGs were significantly correlated with ECM organization, cell adhesion molecule binding, cell–cell adherens junction, collagen-containing ECM, and TGF-beta signaling pathway, which were associated with tumorigenesis and metastasis. Through immune cell infiltration analysis, we compared the infiltration of immune cells in primary and metastatic osteosarcoma specimens. We also constructed the PPI network by STRING database and module analysis; finally, we screened 10 hub genes by cytoHubba including POMC, EGR1, CXCL10, SERPINA1, OAS1, MYC, CXCR4, CXCL2, CHRDL1, and GNAI1. Survival analysis of hub genes based on the TARGET database was also performed in our study.

Based on the results of enrichment analyses, the GO terms and KEGG pathways were explored in osteosarcoma. A substantial body of studies indicated that epithelial-to-mesenchymal transition (EMT) is a process needed for metastasis, during which the loss of cell–cell junction such as adherens junctions, allows tumor cells dissociate from the primary site and acquire the motility to invade stroma (22–25). EMT is also reported to confer resistance to anoikis, which is necessary to survive in the circulation (26). In addition to the involvement of EMT process, the tumor-related ECM is also a key factor in tumor progression. In fact, ECM re-organization such as collagen deposition mediated by collagen-binding integrins may be a general cue for prognosis of tumors (27). Consistent with the above conclusions, the results of GO enrichment analysis, such as ECM organization, cell adhesion molecule binding, cell–cell adherens junction, and collagen-containing ECM, indicate their involvement in the progression of osteosarcoma. Additionally, enrichment of robust DEGs in some KEGG pathways, including TGF-beta pathway and wnt pathway, also demonstrates their relationship with osteosarcoma development. The overexpression of TGF-βs is related with the presence of lung metastasis (28) and is associated with high-grade osteosarcoma (29). Inhibition of wnt pathway can reduce osteosarcoma invasiveness by reversing the EMT (30). On the basis of the above results, we showed that the robust DEGs were highly associated with pathogenesis and progression of osteosarcoma. Furthermore, on the basis of the analysis of modules, we found that three key modules are mainly related to type I interferon signaling pathway, chemokine-mediated signaling pathway, and cell–cell adhesion functions. In the type I interferon pathway, interferon-α has been widely studied. Interferon-α reportedly enhanced the apoptosis of osteosarcoma cells mediated by etoposide (31) and doxorubicin (32). Whether the weakening of interferon signaling pathway plays an important role in the development and metastasis of osteosarcoma deserves further study. It has been reported that chemokine-mediated pathways induce the metastasis of primary tumors to distant target organs. In osteosarcoma, the interaction between chemokine CXCL12 and its receptor CXCR4 drives the metastasis of osteosarcoma cells to the lung (33). CXCR3 and its ligands have the same role in promoting lung metastasis of osteosarcoma (34). In addition, the CXCR7 receptor promotes osteosarcoma lung metastasis (35) and has been recognized as a second receptor with high affinity to CXCL12 (36). The interaction of CXCR7/CXCL12 in the progression of osteosarcoma needs to be further investigated.

In the past few decades, accumulating evidence has indicated that cancer initiation and progression are related with not only cancer itself but also tumor microenvironment (TME) (37, 38). TME is a complex including ECM, exosomes, and stromal cells (39). Among stromal cells, tumor-associated macrophages (TAMs), namely, the M2 type macrophages, have been reported to promote angiogenesis, matrix remodeling (40) and are closely associated with osteosarcoma progression and prognosis (41). In our study, we found that macrophages are the most principal infiltrating immune cells in osteosarcoma including undifferentiated macrophage M0 and macrophage M2; thus, the role of macrophages, especially M2 type macrophages, in the microenvironment of osteosarcoma needs to be further clarified.

Based on the PPI network construction, 10 hub genes were identified. Among these hub genes, six key genes were screened to explore their roles. Serpin peptidase inhibitor clade A member 1 (SERPINA1), a protease inhibitor, was reported to be a predictor in breast cancer (42) and colorectal cancer (43). However, the diagnostic and prognostic roles of SERPINA1 in osteosarcoma were still obscure. Chemokine CXCL2 was reported to be related with prognosis of bladder cancer (44), but its role in the progression of osteosarcoma was still unclear. Early growth response protein 1 (EGR1), a zinc-finger transcription factor, was reported to be involved in cell proliferation and migration (45, 46). And an increasing number of studies have shown that EGR1 is highly associated with cancer development and progression. Liu et al. reported that EGR1 was essential for HNF1A-AS1-mediated cell growth and invasion of gastric cells (47). EGR1 was also reported to promote prostate cancer bone and brain metastasis, as demonstrated by the reduction of blood vessel density in brain and bone caused by decreased EGR1 expression (48). In our study, the expression of EGR1 was higher in metastatic osteosarcoma tissues than primary tissues, but its role in osteosarcoma remained unclear. C-X-C motif ligand 10 (CXCL10) is a member of the CXC subfamily of chemokines and acts through CXC receptor 3 (CXCR3) (49, 50). The prognostic role of serum CXCL10 was proved by Yu et al. in colorectal cancer, and the authors also indicated that the high levels of serum CXCL10 were highly related with liver metastasis (51). Similarly, a study was also reported that high circulating levels of CXCL10 are a biomarker for worse survival in osteosarcoma (52). According to our result, high expression of CXCL10 in osteosarcoma tissues predicted a better survival. The differences in the above conclusions may be due to the different sources of CXCL10; the association between the CXCL10 expression in tissues and overall survival of osteosarcoma patients need to be further studied. Genetic mutations of tumor suppressor such as TP53 and RB1 are highly associated with osteosarcoma development (53, 54). In addition, the mis-regulated expression of oncogene MYC is often found in osteosarcoma patients (55). A study indicated that the overexpression of c-myc promoted osteosarcoma cells invasion through MEK-ERK pathway (56). Consistent with above result, the expression of MYC was higher in metastatic osteosarcoma tissues compared with the primary tissues, and the low expression of MYC predicted a better overall survival in our study. In a word, the MYC expression may act as a biomarker in osteosarcoma metastasis. Chemokine receptor 4 (CXCR4), a seven-transmembrane G protein, has been implicated to mediate the metastasis of several tumors and has become a potential target for tumor therapy (57). Several studies highlighted that the overexpression of CXCR4 potentiated osteosarcoma growth and lung metastasis (58–60). The CXCR4 antibody (61) and antagonist AMD3100 (62) were reported to suppress osteosarcoma cell invasion and lung metastasis, and it revealed that CXCR4 may act as a therapeutic agent to inhibit osteosarcoma progression.

Conclusion

In conclusion, by using integrated bioinformatics analysis such as RRA method, we identified the significant robust DEGs and gene modules in osteosarcoma. The enrichment analyses of DEGs showed that they were closely associated with osteosarcoma development and progression. We not only identified immune cell infiltration in osteosarcoma tissues, but we also screened 10 hub genes. After the above discussion, we found that genes EGR1, CXCL10, MYC, and CXCR4 may be considered as novel biomarkers of osteosarcoma, and more studies need to be done to illuminate their contribution in the diagnosis and prognosis of osteosarcoma.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

TY and JN conceived the study and thoroughly revised the manuscript. JN, WG, WW, ZZ, TR, and YH acquired and analyzed the data. JN, HZ, YY, and XL wrote the manuscript. All authors approved the final version of the manuscript.

Funding

This work was supported by a grant from the National Natural Science Foundation of China (#81572635 and #81972513).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2020.01628/full#supplementary-material

Footnotes

1. ^ https://www.ncbi.nlm.nih.gov/geo/
2. ^ http://www.string-db.org/
3. ^ https://cytoscape.org
4. ^ https://ocg.cancer.gov/programs/target

References

1. Xie BB, Li YR, Zhao RJ, Xu YZ, Wu YH, Wang J, et al. Identification of key genes and miRNAs in osteosarcoma patients with chemoresistance by bioinformatics analysis. Biomed Res Int. (2018) 10:4761064. doi: 10.1155/2018/4761064

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Luetke A, Meyers PA, Lewis I, Juergens H. Osteosarcoma. Treatment–Where do we stand? A state of the art review. Cancer Treat Rev. (2014) 40:523–32. doi: 10.1016/j.ctrv.2013.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Bielack SS, Kempf-Bielack B, Delling G, Exner GU, Flege S, Helmke K, et al. Prognostic factors in high-grade osteosarcoma of the extremities or trunk: an analysis of 1,702 patients treated on neoadjuvant cooperative osteosarcoma study group protocols. J Clin Oncol. (2002) 20:776–90. doi: 10.1200/jco.20.3.776

CrossRef Full Text | Google Scholar

4. Lindsey BA, Markel JE, Kleinerman ES. Osteosarcoma overview. Rheumatol Therapy. (2017) 4:25–43. doi: 10.1007/s40744-016-0050-2

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Damron TA, Ward WG, Stewart A. Osteosarcoma, chondrosarcoma, and Ewing’s sarcoma. Clin Orthopaed Relat Res. (2007) 459:40–7. doi: 10.1097/BLO.0b013e318059b8c9

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Miller BJ, Cram P, Lynch CF, Buckwalter JA. Risk factors for metastatic disease at presentation with osteosarcoma: an analysis of the SEER database. J Bone Joint Surgery Am. (2013) 95:e89–89. doi: 10.2106/jbjs.l.01189

PubMed Abstract | CrossRef Full Text | Google Scholar

7. He XJ, Gao ZZ, Xu HW, Zhang ZW, Fu P. A meta-analysis of randomized control trials of surgical methods with osteosarcoma outcomes. J Orthopaed Surgery Res. (2017) 12:6. doi: 10.1186/s13018-016-0500-0

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Meyers PA, Schwartz CL, Krailo MD, Healey JH, Bernstein ML, Betcher D, et al. Osteosarcoma: the addition of muramyl tripeptide to chemotherapy improves overall survival–a report from the children’s oncology group. J Clin Oncol. (2008) 26:633–8. doi: 10.1200/jco.2008.14.0095

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Zhang KF, Gao JW, Ni Y. Screening of candidate key genes associated with human osteosarcoma using bioinformatics analysis. Oncol Lett. (2017) 14:2887–93. doi: 10.3892/ol.2017.6519

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Diao CY, Xi Y, Xiao T. Identification and analysis of key genes in osteosarcoma using bioinformatics. Oncol Lett. (2018) 15:2789–94. doi: 10.3892/ol.2017.7649

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Yang X, Zhu SM, Li L, Zhang L, Xian S, Wang YQ, et al. Identification of differentially expressed genes and signaling pathways in ovarian cancer by integrated bioinformatics analysis. Oncotar Therapy. (2018) 11:1457–74. doi: 10.2147/ott.s152238

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Ni MW, Liu XK, Wu JR, Zhang D, Tian JH, Wang T, et al. Identification of candidate biomarkers correlated with the pathogenesis and prognosis of non-small cell lung cancer via integrated bioinformatics analysis. Front Genet. (2018) 9:14. doi: 10.3389/fgene.2018.00469

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Hou Q, Bing ZT, Hu C, Li MY, Yang KH, Mo Z, et al. RankProd combined with genetic algorithm optimized artificial neural network establishes a diagnostic and prognostic prediction model that revealed C1QTNF3 as a biomarker for prostate cancer. Ebiomedicine. (2018) 32:234–44. doi: 10.1016/j.ebiom.2018.05.010

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Yang Y, Lu Q, Shao XJ, Mo BH, Nie XQ, Liu W, et al. Development of a three-gene prognostic signature for hepatitis B virus associated hepatocellular carcinoma based on integrated transcriptomic analysis. J Cancer. (2018) 9:1989–2002. doi: 10.7150/jca.23762

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Ritchie ME, Phipson B, Wu D, Hu YF, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:13. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kolde R, Laur S, Adler P, Vilo J. Robust rank aggregation for gene list integration and meta-analysis. Bioinformatics. (2012) 28:573–80. doi: 10.1093/bioinformatics/btr709

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Yu GC, Wang LG, Han YY, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics J Integrat Biol. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol. (2018) 1711:243–59. doi: 10.1007/978-1-4939-7493-1_12

CrossRef Full Text | Google Scholar

19. Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Systems Biol. (2014) 8:7. doi: 10.1186/1752-0509-8-s4-s11

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wang J-S, Duan M-Y, Zhong Y-S, Li X-D, Du S-X, Xie P, et al. Investigating age-induced differentially expressed genes and potential molecular mechanisms in osteosarcoma based on integrated bioinformatics analysis. Mol Med Rep. (2019) 19:2729–39. doi: 10.3892/mmr.2019.9912

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Dai P, He Y, Luo G, Deng J, Jiang N, Fang T, et al. Screening candidate microRNA-mRNA network for predicting the response to chemoresistance in osteosarcoma by bioinformatics analysis. J Cell Biochem. (2019) 120:16798–810. doi: 10.1002/jcb.28938

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Martin TA. The role of tight junctions in cancer metastasis. Semin Cell Dev Biol. (2014) 36:224–31. doi: 10.1016/j.semcdb.2014.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Thiery JP. Epithelial-mesenchymal transitions in tumour progression. Nat Rev Cancer. (2002) 2:442–54. doi: 10.1038/nrc822

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. (2011) 144:646–74. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Yilmaz M, Christofori G. EMT, the cytoskeleton, and cancer cell invasion. Cancer Metastas Rev. (2009) 28:15–33. doi: 10.1007/s10555-008-9169-0

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Yu Y, Elble RC. Homeostatic signaling by cell-cell junctions and its dysregulation during cancer progression. J Clin Med. (2016) 5:26. doi: 10.3390/jcm5020026

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Hamidi H, Ivaska J. Every step of the way: integrins in cancer progression and metastasis. Nat Rev Cancer. (2018) 18:532–47. doi: 10.1038/s41568-018-0038-z

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Yang RS, Wu CT, Lin KH, Hong RL, Liu TK, Lin KS. Relation between histological intensity of transforming growth factor-beta isoforms in human osteosarcoma and the rate of lung metastasis. Tohoku J Exp Med. (1998) 184:133–42. doi: 10.1620/tjem.184.133

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Franchi A, Arganini L, Baroni G, Calzolari A, Capanna R, Campanacci D, et al. Expression of transforming growth factor beta isoforms in osteosarcoma variants: association of TGF beta 1 with high-grade osteosarcomas. J Pathol. (1998) 185:284–9. doi: 10.1002/(sici)1096-9896(199807)185:33.0.co;2-z

CrossRef Full Text | Google Scholar

30. Guo Y, Zi X, Koontz Z, Kim A, Xie J, Gorlick R, et al. Blocking Wnt/LRP5 signaling by a soluble receptor modulates the epithelial to mesenchymal transition and suppresses met and metalloproteinases in osteosarcoma Saos-2 cells. J Orthopaed Res. (2007) 25:964–71. doi: 10.1002/jor.20356

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Yuan XW, Zhu XF, Liang SG, Fan Q, Chen ZX, Huang XF, et al. Interferon alpha enhances etoposide-induced apoptosis in human osteosarcoma U20S cells by a p53-dependent pathway. Life Sci. (2008) 82:393–401. doi: 10.1016/j.lfs.2007.11.025

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Yuan XW, Zhu XF, Huang XF, Sheng PY, He AS, Yang ZB, et al. Interferon-alpha enhances sensitivity of human osteosarcoma U2OS cells to doxorubicin by p53-dependent apoptosis. Acta Pharmacol Sinica. (2007) 28:1835–41. doi: 10.1111/j.1745-7254.2007.00662.x

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Perissinotto E, Cavalloni G, Leone F, Fonsato V, Mitola S, Grignani G, et al. Involvement of chemokine receptor 4/stromal cell-derived factor 1 system during osteosarcoma tumor progression. Clin Cancer Res. (2005) 11:490–7.

Google Scholar

34. Pradelli E, Karimdjee-Soilihi B, Michiels JF, Ricci JE, Millet MA, Vandenbos F, et al. Antagonism of chemokine receptor CXCR3 inhibits osteosarcoma metastasis to lungs. Int J Cancer. (2009) 125:2586–94. doi: 10.1002/ijc.24665

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Goguet-Surmenian E, Richard-Fiardo P, Guillemot E, Benchetrit M, Gomez-Brouchet A, Buzzo P, et al. CXCR7-mediated progression of osteosarcoma in the lungs. Br J Cancer. (2013) 109:1579–85. doi: 10.1038/bjc.2013.482

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Burns JM, Summers BC, Wang Y, Melikian A, Berahovich R, Miao ZH, et al. A novel chemokine receptor for SDF-1 and I-TAC involved in cell survival, cell adhesion, and tumor development. J Exp Med. (2006) 203:2201–13. doi: 10.1084/jem.20052144

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Naito Y, Yoshioka Y, Yamamoto Y, Ochiya T. How cancer cells dictate their microenvironment: present roles of extracellular vesicles. Cell Mol Life Sci. (2017) 74:697–713. doi: 10.1007/s00018-016-2346-3

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Whiteside TL. Exosome and mesenchymal stem cell cross-talk in the tumor microenvironment. Semin Immunol. (2018) 35:69–79. doi: 10.1016/j.smim.2017.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Luo Z, Wang Q, Lau WB, Lau B, Xu L, Zhao L, et al. Tumor microenvironment: the culprit for ovarian cancer metastasis? Cancer Lett. (2016) 377:174–82. doi: 10.1016/j.canlet.2016.04.038

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Aras S, Zaidi MR. TAMeless traitors: macrophages in cancer progression and metastasis. Br J Cancer. (2017) 117:1583–91. doi: 10.1038/bjc.2017.356

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Segaliny AI, Mohamadi A, Dizier B, Lokajczyk A, Brion R, Lanel R, et al. Interleukin-34 promotes tumor progression and metastatic process in osteosarcoma through induction of angiogenesis and macrophage recruitment. Int J Cancer. (2015) 137:73–85. doi: 10.1002/ijc.29376

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Chan HJ, Li H, Liu Z, Yuan Y-C, Mortimer J, Chen S. SERPINA1 is a direct estrogen receptor target gene and a predictor of survival in breast cancer patients. Oncotarget. (2015) 6:25815–27. doi: 10.18632/oncotarget.4441

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Kwon CH, Park HJ, Choi JH, Lee JR, Kim HK, Jo HJ, et al. Snail and serpinA1 promote tumor progression and predict prognosis in colorectal cancer. Oncotarget. (2015) 6:20312–26. doi: 10.18632/oncotarget.3964

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Zhang H, Ye YL, Li MX, Ye SB, Huang WR, Cai TT, et al. CXCL2/MIF-CXCR2 signaling promotes the recruitment of myeloid-derived suppressor cells and is correlated with prognosis in bladder cancer. Oncogene. (2017) 36:2095–104. doi: 10.1038/onc.2016.367

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Wang D, Guan M-P, Zheng Z-J, Li W-Q, Lyv F-P, Pang R-Y, et al. Transcription factor Egr1 is involved in high glucose-induced proliferation and fibrosis in rat glomerular mesangial cells. Cell Physiol Biochem. (2015) 36:2093–107. doi: 10.1159/000430177

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Tarcic G, Avraham R, Pines G, Amit I, Shay T, Lu Y, et al. EGR1 and the ERK-ERF axis drive mammary cell migration in response to EGF. Faseb J. (2012) 26:1582–92. doi: 10.1096/fj.11-194654

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Liu H-T, Liu S, Liu L, Ma R-R, Gao P. EGR1-Mediated transcription of lncRNA-HNF1A-AS1 promotes cell-cycle progression in gastric cancer. Cancer Res. (2018) 78:5877–90. doi: 10.1158/0008-5472.can-18-1011

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Li L, Ameri AH, Wang S, Jansson KH, Casey OM, Yang Q, et al. EGR1 regulates angiogenic and osteoclastogenic factors in prostate cancer and promotes metastasis. Oncogene. (2019) 38:6241–55. doi: 10.1038/s41388-019-0873-8

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Singh UP, Singh S, Iqbal N, Weaver CT, McGhee JR, Lillard JW. IFN-gamma-inducible chemokines enhance adaptive immunity and colitis. J Interfer Cytokine Res. (2003) 23:591–600. doi: 10.1089/107999003322485099

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Neville LF, Mathiak G, Bagasra O. The immunobiology of interferon-gamma inducible protein 10kD (IP-10): a novel, pleiotropic member of the C-X-C chemokine superfamily. Cytokine Growth Fact Rev. (1997) 8:207–19. doi: 10.1016/s1359-6101(97)00015-4

CrossRef Full Text | Google Scholar

51. Toiyama Y, Fujikawa H, Kawamura M, Matsushita K, Saigusa S, Tanaka K, et al. Evaluation of CXCL10 as a novel serum marker for predicting liver metastasis and prognosis in colorectal cancer. Int J Oncol. (2012) 40:560–6. doi: 10.3892/ijo.2011.1247

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Flores RJ, Kelly AJ, Li Y, Nakka M, Barkauskas DA, Krailo M, et al. A Novel prognostic model for osteosarcoma using circulating CXCL10 and FLT3LG. Cancer. (2017) 123:144–54. doi: 10.1002/cncr.30272

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Kempf-Bielack B, Bielack SS, Jurgens H, Branscheid D, Berdel WE, Exner GU, et al. Osteosarcoma relapse after combined modality therapy: an analysis of unselected patients in the cooperative osteosarcoma study group (COSS). J Clin Oncol. (2005) 23:559–68. doi: 10.1200/jco.2005.04.063

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Kwiatkowski N, Zhang T, Rahl PB, Abraham BJ, Reddy J, Ficarro SB, et al. Targeting transcription regulation in cancer with a covalent CDK7 inhibitor. Nature. (2014) 511:616–20. doi: 10.1038/nature13393

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Gamberi G, Benassi MS, Bohling T, Ragazzini P, Molendini L, Sollazzo MR, et al. C-myc and c-fos in human osteosarcoma: Prognostic value of mRNA and protein expression. Oncology. (1998) 55:556–63. doi: 10.1159/000011912

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Han G, Wang Y, Bi W. c-Myc overexpression promotes osteosarcoma cell invasion via activation of MEK-ERK pathway. Oncol Res. (2012) 20:149–56. doi: 10.3727/096504012x13522227232237

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Wong D, Korz W. Translating an antagonist of chemokine receptor CXCR4: from bench to bedside. Clin Cancer Res. (2008) 14:7975–80. doi: 10.1158/1078-0432.ccr-07-4846

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Zhu Y, Tang L, Zhao S, Sun B, Cheng L, Tang Y, et al. CXCR4-mediated osteosarcoma growth and pulmonary metastasis is suppressed by MicroRNA-613. Cancer Sci. (2018) 109:2412–22. doi: 10.1111/cas.13653

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Ren Z, Liang S, Yang J, Han X, Shan L, Wang B, et al. Coexpression of CXCR4 and MMP9 predicts lung metastasis and poor prognosis in resected osteosarcoma. Tumor Biol. (2016) 37:5089–96. doi: 10.1007/s13277-015-4352-8

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Xi Y, Qi Z, Ma J, Chen Y. PTEN loss activates a functional AKT/CXCR4 signaling axis to potentiate tumor growth and lung metastasis in human osteosarcoma cells. Clin Exp Metastas. (2020) 37):173–85. doi: 10.1007/s10585-019-09998-7

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Brennecke P, Arlt MJE, Campanile C, Husmann K, Gvozdenovic A, Apuzzo T, et al. CXCR4 antibody treatment suppresses metastatic spread to the lung of intratibial human osteosarcoma xenografts in mice. Clin Exp Metastas. (2014) 31:339–49. doi: 10.1007/s10585-013-9632-3

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Fontanella R, Pelagalli A, Nardelli A, D’Alterio C, Ierano C, Cerchia L, et al. A novel antagonist of CXCR4 prevents bone marrow-derived mesenchymal stem cell-mediated osteosarcoma and hepatocellular carcinoma cell migration and invasion. Cancer Lett. (2016) 370:100–7. doi: 10.1016/j.canlet.2015.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar