Original Research ARTICLE
Identification of Potential Therapeutic Targets and Immune Cell Infiltration Characteristics in Osteosarcoma Using Bioinformatics Strategy
- 1Musculoskeletal Tumor Center, Peking University People’s Hospital, Peking University, Beijing, China
- 2Beijing Key Laboratory of Musculoskeletal Tumor, Peking University People’s Hospital, Beijing, China
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.
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 database1, 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 database2, and we chose confidence >0.9 as the screening criteria. The visualized PPI network was performed by Cytoscape (version 3.6.1) software3. 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).
The RNA-seq-FPKM data and prognostic information of osteosarcoma patients were downloaded from TARGET database4. 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.
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. 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. 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. 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. 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. 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.
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. 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).
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.
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.
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.
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.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2020.01628/full#supplementary-material
- ^ https://www.ncbi.nlm.nih.gov/geo/
- ^ http://www.string-db.org/
- ^ https://cytoscape.org
- ^ https://ocg.cancer.gov/programs/target
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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.
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
Keywords: bioinformatics analysis, osteosarcoma, robust rank aggregation, immune cell infiltration, survival analysis
Citation: Niu J, Yan T, Guo W, Wang W, Zhao Z, Ren T, Huang Y, Zhang H, Yu Y and Liang X (2020) Identification of Potential Therapeutic Targets and Immune Cell Infiltration Characteristics in Osteosarcoma Using Bioinformatics Strategy. Front. Oncol. 10:1628. doi: 10.3389/fonc.2020.01628
Received: 14 June 2020; Accepted: 27 July 2020;
Published: 21 August 2020.
Edited by:Zhenyu Jia, University of California, Riverside, United States
Reviewed by:Md. Rakibul Islam, Daffodil International University, Bangladesh
Bikash Kumar Paul, Mawlana Bhashani Science and Technology University, Bangladesh
Copyright © 2020 Niu, Yan, Guo, Wang, Zhao, Ren, Huang, Zhang, Yu and Liang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Taiqiang Yan, email@example.com