Identication of Biomarkers and Study of Mechanisms Related to Metastatic of Osteosarcoma Based on Integrated Bioinformatics Analyses

Osteosarcoma (OS) is a serious threat to public health. Because of high morbidity and fairly complicated pathogenesis. The study aim to identify candidate biomarkers and research the molecular mechanisms correlated of patients with metastatic OS. The GSE21257 was downloaded from Gene Expression Omnibus(GEO) database, and the differentially expressed RNAs (DERs) were identied and functional enriched analysis by statistical soft-ware in R. Subsequently, the co-expression modules and its clinical characteristics of OS were identied by weighted gene co-expression network analysis (WGCNA) Following, the KEGG pathways directly related to metastatic OS was to researched by the Comparative Toxicogenomics Database 2019 update (CTD). Finally, the “survival” package in R was used to survival analysis and the DERs were veried using another independent proling GSE14827. A total of 1,464 DERs were classied including 702 up-regulated and 762 down-regulated. In addition, a total of 1248 DERs were obtained by WGCNA analysis, the blue modules is the highest negative correlation (P=0) and the turquoise modules is highest positive correlation (P=3E-196) among all correlations with OS metastatic. The lncRNA-mRNA co-expression network including 4 lncRNAs and 507 mRNAs, and the cytokine-cytokine receptor interaction and JAK-STAT signaling pathway were found signicantly correlation with metastatic. Finally, the increased expression levels of IFNGR1, lower DLEU1 and DLEU2 related to better prognosis. Which were signicantly consistent in the another independent proling GSE14827. each each for A among and a genes adjacency between (35) was used to identify the prognosis-related lncRNAs and mRNAs by COX regression univariate analysis. The Kaplan-Meier survival analysis was used to study the associations between the expression levels and the survival prognosis. The survival curves of samples with low expression of the RNAs and high expression of the RNAs were compared using the log-rank test. P < 0.05 indicated statistically signicant differences.

Screening of signi cant differentially expressed RNAs (DERs) The Limma package (Version 3.4.1, https://bioconductor.org/packages/release/bioc/html/limma.html) (20) in R was used to performed DERs and calculated the FDR values and Fold Change values between metastasis and non-metastasis OS samples. The DERs was identi ed with FDR (false discovery rate) < 0.05 and |log 2 fold change (FC)|> 0.263 were as the signi cance cut-off criteria. The hierarchical clustering analyses of DERs were performed using pheatmap package (Version 3.4.1, https://cran.r-project.org/package=pheatmap) (21) in R, and were presented by two-way hierarchical clustering heatmaps. A P < 0.05 as the signi cance cut-off criteria (22,23).

Functional enrichment analyses
The online tool DAVID (Version 6.8; https://david.ncifcrf.gov/) (24,25) was employed to identify the Gene Ontology (GO) biological processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways that signi cantly associated with the mRNAs in OS samples. A P < 0.05 was considered to achieve signi cant enrichment and have statistical signi cance.
Weighed gene co-expression network analysis (WGCNA) To identify modules associated with OS and its clinical characteristics by using the WGCNA (Version1.61, https://cran.r-project.org/web/packages/WGCNA/) (26,27) package in R (28,29). In this study, all RNAs were detected and performed network construction and division. Based on the clinical information of OS samples to calculated the correlation between each partitioning module and each clinical indicator, with the minimal module size of 50 and the merge cut height of 0.99. Functional annotation was conducted for each stable module using the WGCNA package (30), with a Fold enrichment > 1. Additionally, DERs was performed for each module with a P-value and false discovery rate (FDR) of < 0.05. A cluster dendrogram among modules and a genes adjacency heatmap between modules were generated.

LncRNA-mRNA co-expression network construction
The co-expression analysis of lncRNAs and mRNAs was conducted based on the Pearson correlation coe cient (PCC) of the Cor function in the R language (Version 3.4.1, http://77.66.12.57/R-help/cor.test.html). Their expression levels was conduct a network visualization display via Cytoscape (Version 3.6.1, http://www.cytoscape.org/) (31). Then, Gene Set Enrichment Analysis (GSEA, http://software.broadinstitute.org/gsea/index.jsp) (32) in R was used to performe KEGG pathway enrichment analysis on optimal mRNAs in the lncRNA-mRNA network. The enrichment score (ES), normalised enrichment score (NES) and nominal P value were used in this analysis, and the NES absolute value increases and the P value decrease, suggesting a higher degree of enrichment and a higher signi cance of the result (33). A P < 0.05 was considered to screen KEGG pathways that were signi cantly enriched in the relevant mRNAs OS relevant KEGG pathway network construction The Comparative Toxicogenomics Database 2019 update (CTD, http://ctd.mdibl.org/)(34) was used to research the KEGG pathways directly related to OS with the keywords of "osteosarcoma". Obtained the overlapping pathways by compared with the RNAs of signi cantly participated in pathways from the coexpression network, and to constructed a OS relevant pathway network.
Survival analysis and DERs veri cation.
The R package "survival" (http://bioconductor.org/packages/survivalr/) (35) was used to identify the prognosis-related lncRNAs and mRNAs by COX regression univariate analysis. The Kaplan-Meier survival analysis was used to study the associations between the expression levels and the survival prognosis. The survival curves of samples with low expression of the RNAs and high expression of the RNAs were compared using the log-rank test. P < 0.05 indicated statistically signi cant differences.
In addtion, to verify the reliability of the important RNAs of the metastasic and non-metastasic OS samples. Another independent pro ling of GSE14827 containing 27 OS samples (18 non-metastases and 9 metastases) was downloaded from GEO (https://www.ncbi.nlm.nih.gov/) using the information supported by the platform GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array.

Survival analysis and DERs veri cation
According to the expression levels of 13 RNAs, the OS samples were classi ed into high-expression and low-expression group. The Kaplan-Meier survival analysis showed that samples with high expression levels of IFNGR1 had better prognosis (P = 3.823e-02), while reduced the expression levels of DLEU1 (p = 4.345e-02) and DLEU2 (P = 2.933e-02) related to better prognosis (Fig. 5).
The expression level of IFNGR1, DLEU1 and DLEU2 of the GSE21257 was signi cantly consistent in the metastatic and non-metastatic OS samples of another independent pro ling GSE14827 (Fig. 6).

Discussion
OS is a highly aggressive bone tumor, which early systemic metastasis in young people, and results in poor survival for patients with OS (36). Therefore, it is very important to develop new therapeutic strategies and to treat the metastatic OS patients. However, the molecular mechanisms underlying the metastasis of OS are still very urgently to understand. This study aimed to investigate the metastatic processes of OS and to identify potential biomarkers of metastatic OS.
The gene expression pro le of GSE21257 using bioinformatics analysis and shown that there were 762 up-regulated and 702 down-regulated DERs which may be potential biomarkers of metastatic OS. Furthermore, The results of GO term and KEGG enricment analysis were demonstrated that the DERs previously been associated with the cell adhesion molecules (CAMs), antigen processing and presentation, lysosome (2). The functional enrichment analysis suggested that aberrant regulation of RNAs may contribute to metastasis of OS. The DERs identi ed in this study may be associated with metastatic processes of OS.
They can be used as a criterion for malignancy phenotype of bone cancer and prediction of metastatic potentiality of OS (37). Furthermore, the WCGNA results showed the blue module (1 lncRNAs and 495 mRNAs) and the turquoise module (5 lncRNAs and 529 mRNAs) have sni cantly relevance with metastatic OS. The highly stable blue and turquoise modules were mainly involved in cellular immune responses, cell adhesion and the cell cycle, and indicating their likely association with metastasis OS pathogenesis. In addition, 13 RNAs with cytokine-cytokine receptor interaction and JAK-STAT signaling pathway RNAs-pathway co-expression network was constructed. The Kaplan-Meier curve analysis was performed to explore the associations between overall survival and these 13 essential RNAs in patients with OS. Consequently, IFNGR1, DLEU1 and DLEU2 were found to be signi cantly associated with metastatic OS (P < 0.05). The three prognostic RNAs have all been reported to be associated with tumor, suggesting the reliability of the methods we used in this study.
To our surprise, IFNGR1 expression level was increased in metastasis OS group ( Figure. 5). Since IFNGR1 have been identi ed in 1996, IFNGR1 plays a importent role in the immune response against mycobacteria (38). IFNGR1 was also found to be mainly mutated in non-small-cell lung cancer and found that interferon associated IL2-STAT5 pathways were signi cantly enriched (39). Importantly, IFNGR1 plays an important role in the cytokine-cytokine receptor interaction and JAK-STAT signaling pathway simultaneously, which are reportedly involved in up-regulated to against OS. So it may be involved immune response for OS and as a protective factor in the metastasis of OS (40).
LncRNA DLEU1 was acted as an oncogene in various types of cancer. DLEU1 was overexpressed and promoted the cell proliferation, migration and invasion in DLEU1 was up-regulated in samples with metastasis OS. Thus, the OS genes may identi ed as signi cantly gene in metastasis OS (42). DLEU2 is the host gene of miR-15a, the study was indicated the transcription of DLEU2 was repressed by hypoxia. These ndings indicated that miR-15a may be a valuable target for the treatment of OS (43), suggesting the carcinogenic role of DLEU2 in OS. DLEU1 and DLEU2 elicited OS cell metastasis by inhibiting the activation of the JAK-STAT signaling pathway, thus, they are considered as potential target genes for cancer treatment.

Conclusion
In conclusion, the increased expression levels of IFNGR1 may be as a protective factor in metastastic OS. Meanwhile, high DLEU1 and DLEU2 expression predicts poor overall survival of OS patients with metastasis, and they were signi cantly consistent in the metastatic and non-metastatic OS samples of another independent pro ling GSE14827. The IFNGR1, DLEU1 and DLEU2 all played important roles in OS with metastasis by in uencing cytokine-cytokine receptor interaction and JAK-STAT signaling pathway. Therefore, our study revealed that the IFNGR1, DLEU1 and DLEU2 might provide a potential new therapeutic strategy for metastatic OS treatment. The raw data were collected and analyzed by the Authors, and are not ready to share their data because the data have not been published.

Abbreviations
Competing interests: The authors declare no con ict of interest.