Identification of potential micro-messenger RNAs (miRNA–mRNA) interaction network of osteosarcoma

ABSTRACT Osteosarcoma (OS) is the most common primary malignant tumor in children and adolescents. Numerous studies have reported the importance of miRNA in OS. The purpose of this study is to predict potential biomarkers and new therapeutic targets for OS diagnosis and prognosis by analyzing miRNAs of OS plasma samples from the Gene Expression Omnibus (GEO) database. Data-sets were downloaded from the GEO and analyzed using R software. Different expressions of miRNAs (DE-miRNAs) in plasma and mRNAs (DE-mRNAs) in OS patients were identified. Funrich was used to predict the transcription factors and target genes of miRNAs. By comparing the target mRNAs and DE-mRNAs, the intersection mRNAs were identified. The intersection mRNAs were imported to perform Gene Ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. MiRNA-mRNA regulatory network and a protein-protein interaction (PPI) network were constructed by using Cytoscape. Finally, a total of 164 DE-miRNAs, 256 DE-mRNAs, and 76 intersection mRNAs were identified. The top 10 TF of up- and down-regulated DE-miRNAs were also predicted. In addition, GO and KEGG analyses further revealed the intersection mRNAs. By constructing the miRNA–mRNA networks, we found miR-30d-5p, miR-17-5p, miR-98-5p, miR-301a-3p, and miR-30e-5p were the central hubs. COL1A1, COL1A2, MMP2, CDH11, COL4A1 etc. were predicted to be the key mRNA by constructing the PPI networks. Through a comprehensive bioinformatics analysis of miRNAs and mRNAs in OS, we explored the potential effective biomarkers and novel therapeutic targets for the diagnosis and prognosis of OS.


Introduction
Osteosarcoma (OS) is the most common primary malignant bone tumor in children and adolescents, which has high mortality and disability rate [1][2][3]. With the progression of surgery, neoadjuvant chemotherapy, etc., the 5-year survival rate of patients with localized osteosarcoma has increased from less than 20% to 70%-80% [2,4]. However, the prognosis is still poor for patients with lung metastasis or recurrence. The 5-year survival rate of these patients is only 20% [4,5]. Therefore, it is necessary to further explore its potential pathogenesis.
Bioinformatics is a new field of biological research, processing, and analysis of biological data by using mathematical, statistical, and computational methods. Due to the massive amounts of data generated by new technologies, such as genomic sequencing and microarray chips, the traditional gene-by-gene approaches are insufficient to meet the growth and demand of biological research. Therefore, bioinformatics is a valuable way to achieve biological understanding and therapeutic progress [13].
MiRNAs play an important role in OS development and metastasis. Therefore, the aim of this study was to evaluate and summarize the evidence relating to the miRNA in OS plasma to determine the most effective diagnosis, treatment, and prognosis evaluation strategy.

Data processing of DE-miRNAs
We selected datasets from GEO (http://www. ncbi.nlm.nih.gov/geo), which is a publicly available database of gene and microarray profiles. The search strategy ('osteosarcoma' [MeSH Terms] and miRNA [All Fields]) AND ('Homo sapiens'[Organism] AND 'Expression profiling by RT-PCR'[Filter]) were adopted until August 2020. Inclusion criteria were as follows: plasma of miRNA from OS patients or healthy people. In the end, the data-set GSE65071 based on the platform of GPL19631 (Exiqon human V3 microRNA PCR panel I+ II) was chosen for analysis. Background corrections and normalizations were processed using the R package 'affy' from the Bioconductor project. The DE-miRNAs analysis was conducted using the R package 'limma'. The P value <0.05 and |log2FC| >1 was set as the threshold for DE-miRNAs.

Identification of TF and targeted genes of DE-miRNAs
The FunRich (http://www.funrich.org) provides annotations on pathways, TF, biological processes (BP), cellular components (CC), molecular functions (MF), etc. [14]. In the current study, up-regulated and down-regulated miRNAs were separately subjected to FunRich, top 10 TF of DE-miRNAs and targeted genes were separately identified.

Data processing of DE-mRNAs
We selected GEO and the search strategy ('osteosarcoma' [MeSH Terms] and mRNA [All Fields]) AND ('Homo sapiens'[Organism] AND 'Expression profiling by array'[Filter]) were adopted until August 2020. Inclusion criteria were as follows: mRNAs from patients with OS or healthy people. In the end, the data-set GSE16088 based on the platform of GPL96 (Affymetrix Human Genome U133A Array) was chosen for analysis. Background correction, normalization and log2 transformation were processed by the R package 'affy' from the Bioconductor project. The differential analysis was conducted using the R package 'limma'. The P value <0.05 and |log2FC| >2 was set as the criterion for DE-mRNAs.

Intersection genes of target mRNAs and DE-mRNAs
The target genes of DE-miRNAs predicted by FunRich were compared with DE-mRNAs predicted by GEO. The intersection genes were presented as Venn diagrams and the overlap mRNAs were identified. It is widely acknowledged that there is an inverse relationship between miRNA and the expression of target genes. As a result, upregulated mRNA and downregulated mRNA were defined.

Functional annotation and pathway enrichment analysis
GO functional annotation is a widely used bioinformatics tool for analyzing functional relationships between gene products, including three categories: CC, BP, and MF. KEGG analysis is used to study the enrichment pathway of cross genes, so as to further understand gene function. GO functional annotation and KEGG pathway enrichment analysis were performed by the R package 'cluster Profiler' and 'enrich plot' [14]. Fisher's exact test was used to classify the GO category and select the significant KEGG pathway, P value < 0.05 and q value < 0.05 were considered a significant difference.

Construction of miRNA-mRNA network and Protein-protein interaction network
The miRNA-mRNA network of the mRNAs was constructed by Cytoscape v3.8.0 [15,16]. The PPI network for dysregulated mRNAs was established via the Search Tool for the Retrieval of Interacting Genes database (STRING, http:// string-db.org). Then, the PPI network was constructed by Cytoscape. Based on the combined score>0.4 for PPI pairs of DEGs in our study, we further built a PPI network model using Cytoscape software.

Identification of potential Hub mRNAs and hub miRNAs
We calculated the node number of each miRNA to identify hub regulators in the above miRNA-mRNA network. MiRNAs with more than 5 nodes in the network were selected as hub-miRNAs. Higher degree nodes play a vital role in maintaining the stability of the whole network.

Results
The goal of this study was to evaluate and summarize the evidence relating to the miRNA in OS plasma to determine the most effective diagnosis, treatment, and prognosis evaluation strategy. We predicted TF of miRNAs, pathways hub miRNAs, and hub mRNAs and Constructed of miRNA-mRNA network and PPI network.  Table 1A and 1B.      down-regulated mRNAs were also present to perform GO (Figure 7a, Figure 7b) and KEGG pathway enrichment analysis. But this list of down-regulated mRNAs had no outcome of the KEGG pathway.

Discussion
Previous studies have reported that miRNAs play an important role in the occurrence and development of OS, but the interaction between miRNAs and mRNA is still unclear [5]. The miRNA-mRNA regulatory network plays an important role in the occurrence and development of cancer. Based on the GSE65071 profile dataset, 15 primary normal plasmas, 20 OS plasmas, were enrolled and analyzed. A total of 256 DE-miRNAs (191 up-and 65 down-regulated DE-miRNAs) were identified. Five miRNAs including miR-30d-5p, miR-17-5p, miR-98-5p, miR-301a-3p, and miR-30e-5p were identified as the hub miRNAs in our network. Liao et al. found miR-98-5p was sponged with LncRNA SNHG16 to regulate cellular processes in OS [17]. Many studies have reported that microRNA-301a promoted cell proliferation in OS and other sarcomas [18][19][20][21]. The normal expression of miR-30d-5p could inhibit cancer proliferation, migration, and invasion of non-small cell lung cancer [22]. It was reported that a low expression of miR-30e-5p was associated with many kinds of cancers, such as nasopharyngeal carcinoma, esophageal cancer, and bladder cancer [8,23,24]. However, a previous study reported miR-17-5p is down-regulated in the OS patients, which is contrary to our results [7]. It may be related to different types of specimens.
TFs and miRNAs regulate mRNA gene expression. Additionally, miRNAs and TF could alter the expression of each other [25]. We predicted the TF of DE-miRNA using FunRich software such as EGR1, Sp1, SP4, POU2F1. Early growth response protein-1 (EGR1) is a zinc-finger protein, belonging to the EGR family [9]. Cell differentiation and mitogenesis require the participation of target gene products activated by the EGR family [26]. EGR1 is reported to be a cancer suppressor gene and the downregulating of EGR1 in OS samples was usually associated with poor prognosis [9]. Fiorillo et al. found EGR1 could regulate MG-specific miRNAs, miR-21-5p and miR-30e-5p [27]. Sp1 is related to many cellular processes, including cell growth, cell differentiation, immune responses, apoptosis, chromatin remodeling, and response to DNA damage [28--28-30]. Xing et al. reported that SP1 promoted OS progression via the miR-655/SOX18 axis [28]. Gao et al. found that Sp1-mediated up-regulation of lnc00152 promotes invasion and metastasis of retinoblastoma cells via the miR-30d/SOX9/ZEB2 pathway [31]. POU2F1 is also known as OCT-1, located on chromosome 1q24 [32]. The POU2F1 transcription factor belongs to the POU   transcription factor family that contains the POU domain with a necessary amino acid region for DNA binding to the octameric sequence ATGCAAAT [33]. Li et al. reported that LncRNA SND1-IT1 via sponging miRNA-665 upregulated POU2F1 to accelerate the proliferation and migration of OS [32].
Using target genes for DE-miRNAs overlapped with DE-mRNAs from GSE16088 profile dataset, 55 up-, and 21 downregulated target mRNAs were identified. According to GO analysis, the 76 mRNAs enriched in osteoblast differentiation, focal adhesion, and integrin binding. Bone remodeling is an important process, which is the balancing activities of the bone-forming osteoblasts and the bone-resorbing osteoclasts [34]. Impairments in these balanced activities may result in osteo-condensing bone pathologies, such as osteoporosis, osteopetrosis, and can also be the origin of bone cancers. The suspected initiated cell in OS is the osteoblast [35]. Tanaka et al. reported that miR-138 could inhibit Ewing's sarcoma cells metastatically by targeting focal adhesion kinase [36]. Liu et al. found that miR-128 inhibited epithelial-mesenchymal transition of human OS cells by directly targeting integrin α2 [37]. KEGG analysis found that these genes were mainly involved in ECM−receptor interaction, protein digestion, and absorption, PI3K−Akt signaling pathway, Relaxin signaling pathway, and Focal adhesion. The occurrence, development, invasion, and metastasis of malignant tumors are often related to the abnormal expression of ECM and ECM−receptor [38,39]. The tumor cells activate or secrete protein-degrading enzymes to degrade the matrix after adhering to various components of the ECM through their surface receptors, thereby forming a local lysis zone, which constitutes a channel for tumor cell metastasis. Chen et al. reported miR-191-5p promoted the development of OS via targeting EGR1 and activating the PI3K/AKT signaling pathway [9]. These different pathways may be the potential mechanism of OS.   In order to construct a PPI network of target genes, we screened out the top 20 genes. There were COL1A1, COL1A2, MMP2, CDH11, etc. COL1A1 has a triple helix structure, comprising two alpha1 chains and one alpha2 chain, which encodes the pro-alpha1 chains of type I collagen. It is associated with a particular type of skin tumor, and is promoted through overproduction of platelet-derived growth factor (PDGF), such as dermatofibrosarcoma [40]. Bi et al. reported that MicroRNA-98 inhibited the cell proliferation of human hypertrophic scar fibroblasts by targeting Col1A1 [41]. In clinical melanoma samples, the upregulation of COL1A1 was negatively correlated with disease-free survival [42]. MMP2 is a member of the MMP gene family, which has the ability to cleave ECM components and participate in signal transduction. MMPs and their endogenous inhibitors (tissue inhibitors of metalloproteinases, TIMPs) are the most important substances in the regulation of the ECM degradation process [39]. The imbalance of MMP-TIMP results in associating tumor invasion by extracellular matrix proteolysis. The previous study has reported overexpression of MMP-2 and MMP-3 promoted OS cell migration [43]. Gong et al. found that up-regulated miR-17 in ovarian cancer cells could decrease the production of activated MMP-2 and inhibit ovarian cancer cell peritoneal metastasis [44].In addition, Jia et al. reported downregulation of TWIST1 decreased cell viability, inhibited migration, and promoted apoptosis of OS cells [45], Sun et al. found miR-143-3p inhibited the proliferation, migration, and invasion of OS cells by targeting FOSL2 [46]. Croset et almiR-30 could inhabit tumor cell invasion by silencing CDH11 or ITGA5 in ER-/PR-negative breast cancer [47]. By analyzing the miRNA of the blood of patients with OS, we can predict potentially effective biomarkers associated with the invasion and metastasis of OS, which provides a new perspective on early diagnosis and finding novel therapeutic targets for OS patients. However, the present study has several limitations. First of all, this study lacks more convincing evidence such as qPCR and immunohistochemistry results. Secondly, there are only sarcoma data   in the gepia database, while there are only two normal samples in the normal group so we cannot provide the expression level of hub genes between OS and normal samples in the public database. Third, the present data was obtained from the GEO database, which cannot evaluate the reliability and quality of statistical data. Fourth, this study lacks prognosis information from public databases.

Conclusions
Through a comprehensive bioinformatics analysis of miRNAs and mRNAs in OS, we explored the potential effective biomarkers and novel therapeutic targets for the diagnosis and prognosis of OS.

Highlights
(1) We predicted underlying TF and target genes of DE-miRNAs from OS plasma samples.
(2) Based on the blood samples of OS patients, GO functional annotation and KEGG pathway enrichment analysis were performed.
(3) Based on the blood samples of OS patients, MiRNA-mRNA network, and PPI network were constructed.

Availability of data and material
All data are fully available without restriction.

Disclosure statement
The authors declare that they have no competing interests.

Funding
This study was supported by the National Natural Science