Screening and Analysis of Biomarkers in the miRNA-mRNA Regulatory Network of Osteosarcoma

Osteosarcoma is a malignant disease, and few effective strategies can completely overcome the prognosis of these patients. This study attempted to reveal the key factors and related molecular mechanisms of osteosarcoma via excavating public microarray datasets. The data were obtained from the Gene Expression Omnibus (GEO) database; the differentially expressed miRNAs and differentially expressed genes were obtained in GSE69470 and GSE12685l, respectively; the target of miRNAs were predicted with the miRDIP database; the functions of the factors were analyzed and visualized by the David database and R language, respectively. Moreover, the protein-protein interaction network and miRNA-mRNA network were performed with the STRING database and Cytoscape software to identify the hub nodes in GSE69470 and GSE12685. The results showed that 834 DEGs were found in GSE12685 and 37 miRNAs were found in GSE69470. Moreover, the target of 37 miRNAs were enriched in PI3K/AKT, P53, Wnt/β-catenin, and TGF-β pathways and related with skeletal system development and cell growth. Besides, the miRNAs including miR-22-3p, miR-154-5p, miR-34a-5p, miR-485-3p, miR-93-5p, and miR-9-5p and the genes including LEF1, RUNX2, CSF1R, CDKN1A, and FBN1 were identified as the hub nodes via network analysis. In conclusion, this study suggested that the miRNAs including miR-22-3p, miR-154-5p, miR-34a-5p, miR-485-3p, miR-93-5p, and miR-9-5p and the genes including LEF1, RUNX2, CSF1R, CDKN1A, and FBN1 act as key factors in the progression of osteosarcoma.


Introduction
Osteosarcoma is a serious bone malignancy with a high tendency of aggressiveness and metastases, which has a high incidence in childhood [1,2]. Statistically, more than 60% of patients with osteosarcoma are under 25 years old [3]. At present, surgery combined with neoadjuvant chemotherapy is the pillar for osteosarcoma treatment [4]. In recent ten years, the application of multiple radiotherapy and chemotherapy techniques and mature modern surgical skills have effectively changed the prognosis of patients [5]. With the current therapy strategies, a 5-year survival rate of patients has increased to 70% [6]. However, many patients still find it difficult to escape the bad outcomes in the long term. Although accumulating research has continually figured out the potential reasons of osteosarcoma development, the detailed mechanisms of osteosarcoma remain unclear, and deep research is urgently necessary.
MicroRNAs play important roles in the metabolic activities of cells via regulating the expression of proteins [7,8]. Several reporters have indicated that the miRNA profile of osteosarcoma exhibits a high difference with the normal tissues, and some miRNAs have been identified as the hub nodes in regulating the malignant behavior of the tumor cells, such as rapid proliferation and high invasion [9]. For osteosarcoma, the miRNA profiles of patients also exhibit significant differences compared with normal people. Microarray analysis has been frequently used for revealing RNA profiles of diseases, which is useful for analyzing the related pathological mechanisms [10]. Some research has revealed the related molecular mechanism of some tumors via microarray analysis [11,12].
is study aimed to analyze the miRNA profile and identify the key miRNAs of osteosarcoma via the GEO database, reveal the functions and potential regulation mechanisms of the miRNAs, and finally provide some novel references for osteosarcoma treatment.

Data Source.
e microarray datasets of GSE12685 and GSE69470 were obtained from the GEO database (https:// www.ncbi.nlm.nih.gov). e raw data of the datasets were downloaded with the packages of the R language including GEOquery, dplyr, and Limma.

Identification of Differentially Expressed Genes.
All datasets were analyzed with the GEO2R tools of GEO, and the expressed matrix files of the datasets were obtained. For GSE12685, 11 clinical pathological samples of osteosarcoma tumor and 2 differential normal human osteoblasts were used to identify the DEGs. For GSE69470, the DEGs were analyzed in human osteosarcoma cell lines and four normal human cell lines including chondrocyte cells, mesenchymal stem cells, osteoblasts, and skeletal muscle cells. e genes with |logFC| > 2 and P < 0.05 were considered as the differentially expressed genes (DEGs).

Screening of Potential Targets of the miRNAs.
e common genes in the differentially expressed miRNAs in the matrix files were selected and then used for the next analysis. For GSE69470, the common differentially expressed miRNAs in the matrix files are the osteosarcoma tumor cells and differentially normal human osteoblast cells, osteosarcoma tumor cells, and mesenchymal stem cells (MSCs). e potential targets of common miRNAs were predicted with the miRDIP database (https://ophid.utoronto.ca/mirDIP/). e genes with a high minimum score (top 5%) were selected for the targets of the common miRNAs. Moreover, the common genes of the matrices were screed and visualized with VENNY (https://bioinfogp.cnb.csic.es/tools/venny/ index.html).

Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis.
A KEGG enrichment analysis was performed for investigating the related pathways of osteosarcoma. For GES12685, the DEGs were analyzed with KEGG enrichment. For GSE69470, the potential targets of differentially expressed miRNAs were used for KEGG enrichment analysis. e genes were annotated and analyzed by the David database (https://david.ncifcrf.gov/), and the related pathways in the results with a P value <0.05 were visualized with the R language.

Gene Ontology (GO) Enrichment
Analysis. GO enrichment analysis was used to identify the molecular functions of genes. For GES12685, the DEGs were annotated by the David database, and then the related ENTREZID of genes were analyzed with the packages of the R language including clusterProfiler and org.Hs.eg.db. For the miRNAs, the targets of them were analyzed with GO enrichment. Finally, the related functions module with a P value <0.05 were visualized with the R language.

Network Analysis.
e PPI network analysis was performed to identify the hub nodes of osteosarcoma. In brief, the common genes of the targets of the miRNAs in GSE69470 and the DEGs in GSE69470 were analyzed with the STRING database (https://cn.string-db.org/), and then visualized with Cytoscape software. For the miRNA-mRNA network, the connection of the differentially expressed miRNAs in GSE69470 and the DEGs in GES12685 were structured with Cytoscape. For the miRNA-mRNA network, the differentially expressed miRNAs in GSE69470 and their targets were visualized by Cytoscape software (3.4.0).

Identification of miRNAs of the Datasets.
e datasets were used to analyze the differentially expressed miRNAs. For GSE12685, 11 clinical pathological samples of osteosarcoma tumor and 2 healthy samples of normal human osteoblast were used to identify the DEGs, and 424 downregulated genes and 410 upregulated genes were found in osteosarcoma (Figures 1(c) and 2(c)). For GSE69470, the differentially expressed miRNAs were analyzed in human osteosarcoma cell lines and four normal human cell lines including chondrocyte cells, mesenchymal stem cells, osteoblasts, and skeletal muscle cells. Compared with the expression profiles of chondrocyte, osteoblasts, and skeletal muscle cells, 45 downregulated genes and 10 upregulated genes were found in human osteosarcoma cell line (Figure 1). Compared with the expression profile of human mesenchymal stem cells, 43 downregulated genes and 16 upregulated genes were found in human osteosarcoma cell lines ( Figure 1). Moreover, 37 common miRNAs were found in the matrices of GSE69470.

Identification of the Targets of the Common miRNAs.
To identify the key miRNAs in the development of osteosarcoma, the 37 common miRNAs including 30 downregulated miRNAs and 7 upregulated miRNAs were screened in the differentially expressed miRNAs of osteosarcoma cell lines compared with normal differential cells or mesenchymal stem cells (Figure 2(c)). 37 common miRNAs were screened in GSE69470, and 4037 potential targets were predicted by the miRDIP database. Moreover, 114 common genes were screened from the DEGs of GSE12685 and the targets of miRNAs in in GSE69470 (Figures 3(c) and 3(d)).

Functional Modules Analysis.
To investigate the pathological mechanism of osteosarcoma, the targets of the screened miRNAs in GSE69470, the DEGs in GSE12685, and the common genes of them were analyzed with the KEGG enrichment. e DEGs in GSE12685 were associated with the extracellular matrix organization, regulation of vasculature development, skeletal system development, and so on (Figure 3(a)). For GSE69470, it was found that the targets of the screened miRNAs were involved in the skeletal system development, covalent chromatin modification, cell growth, and so on (Figure 3(b)). Moreover, it was found that the common genes were related with skeletal system development, regulation of cellular response to growth factor stimulus, regulation of cell growth, and so on (Figure 3(e)).

Pathways Analysis.
To reveal the regulation mechanism of osteosarcoma, the related pathways of the factors in GSE12685 and GSE69470 were investigated with KEGG enrichment. e results showed that the DEGs in GSE12685 were related with the PI3K-AKT pathway, the Hippo pathway, the MAPK pathway, and microRNAs in cancer, P53 pathways, the Rap1 pathway, the TGF-β pathway, and so on (Figure 4(a)). e targets of the screened miRNAs in GSE69470 were connected with pathways in cancer, the PI3K-AKT pathway, the MAPK pathway, and microRNAs in cancer, the Rap1 pathway, the P53 pathways, the TGF-β pathway, and so on (Figure 4(b)). Moreover, it was also found that the common genes of the targets in GSE69470 and DEGs in GSE12685 were also involved in the PI3K/AKT pathway, the Hippo pathway, the MAPK pathway, and microRNAs in cancer, P53 pathways, the TGF-β pathway, and the Rap1 pathway (Figure 4(c)).
Besides, miR-485-3p has been also proved to be involved in the cancer development via regulating TGF-β pathways, and RUNX2 has been widely confirmed as the downstream factors of TGF-β. In this study, RUNX2 was identified as the potential target of miR-485-3p.
Data Availability e data to support the findings of this study are available on reasonable request from the corresponding author.