miRNA-Gene Interaction Network Construction Strategy to Discern Promising Traditional Chinese Medicine against Osteoporosis

Osteoporosis is a widespread bone disease that affects million cases annually. The underlying mechanisms behind the progress of osteoporosis remain enigmatic, which limits detections of biomarkers and therapeutic targets. Hence, this study was aimed at exploring hub molecules to better understand the mechanism of osteoporosis development and discover the traditional Chinese medicine potential drugs for osteoporosis. miRNA and gene expression profiles were downloaded from Gene Expression Omnibus (GEO). Weighted correlation network analysis (WGCNA) was used to identify the key modules for osteoporosis. DIANA Tools was applied to perform pathway enrichment. A miRNA-gene interaction network was constructed, and hub miRNAs and genes were distinguished using Cytoscape software. Receiver operating characteristic (ROC) curves of hub miRNAs and genes were plotted, and correlations with hub genes and osteoporosis-associated factors were evaluated. Potential drugs for osteoporosis in Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP) were screened, and molecular docking models between these drugs and target genes were showed by AutoDock tools. Two hub modules, 1 miRNA module and 1 gene module, were identified to be the most strongly correlated with osteoporosis by using WGCNA. Then, 3 KEGG pathways including focal adhesion, PI3K-Akt signaling pathway, and gap junction were shared pathways enriched with the miRNAs and genes screened out by WGCNA and differential expression analyses. Finally, after constructing a miRNA-gene interaction network, 6 hub miRNAs (hsa-miR-18b-3p, hsa-miR-361-3p, hsa-miR-484, hsa-miR-519e-5p, hsa-miR-940, and hsa-miR-1275) and 6 hub genes (THBS1, IFNAR2, ARHGAP5, TUBB2B, FLNC, and NTF3) were detected. ROC curves showed good performances of miRNAs and genes for osteoporosis. Correlations with hub genes and osteoporosis-associated factors suggested implicational roles of them for osteoporosis. Based on these hub genes, 3 natural compounds (kainic acid, uridine, and quercetin), which were the active ingredients of 192 herbs, were screened out, and a target-compound-herb network was extracted using TCMSP. Molecular docking models of kainic acid-NTF3, uridine-IFNAR2, and quercetin-THBS1 were exhibited with AutoDock tools. Our study sheds light on the pathogenesis of osteoporosis and provides promising therapeutic targets and traditional Chinese medicine drugs for osteoporosis.


Introduction
Osteoporosis is the most common bone disease with clinicopathological features of low bone mineral density (BMD) and deterioration of bone microstructure [1]. Osteoporosis affects up to 40% women, and more than 3 million cases in US suffer from osteoporosis per year [2]. Bone is a living tissue in which the bone formation by osteoblast and bone resorption by osteoclast are in critical and dynamic balance [1]. Despite the universality of osteoporosis, the underlying mechanisms related to its progress are still ambiguous, which severely limits detections of biomarkers and therapeutic targets.
MicroRNAs (miRNAs) are short noncoding RNAs able to regulate gene expression by interacting with the 3′ untranslated region of target mRNA. Multiple studies have indicated that various miRNAs are involved in the process of bone metabolism and remodeling. Numerous miRNAs are dysregulated in patients with osteoporosis, including miR-194, miR-16, miR-338, miR-2861, miR-223, miR-422, miR-133, miR-21, miR-186, and miR-3651 [3][4][5][6][7][8]. For example, miR-194 was upregulated in osteoporosis patients and may be a promising biomarker for osteoporosis [4]. In addition, miR-16 was suggested to inhibit osteogenesis by suppressing the expression of VEGFA [5]. Moreover, miRNAs could regulate several signaling pathways involved in osteoporosis. For instance, miR-133 and miR-203 promote osteogenic differentiation [9,10]. However, there is no study comprehensively analyzing hub miRNAs. In this study, we aimed to screen key miRNAs and their target genes by using comprehensive bioinformatics analysis.
Currently, high-throughput sequencing and bioinformatics technology has dramatically improved clinical, biological, and drug research. Weighted correlation network analysis (WGCNA) was a widely used computational tool for identification of hub genes [11]. In the present study, we constructed a coexpressed network with WGCNA and identified novel hub miRNAs and genes associated with osteoporosis. It is well established that traditional Chinese herbal medicine (TCHM) has been widely applied for thousands of years to treat diverse diseases, including osteoporosis [12]. However, multiple TCHMs have been reported to exert favourable antiosteoporosis actions through combining with Western medicine [13]. Herewith, it would be ponderable to distinguish the effective TCHM compounds used for treating osteoporosis. Based on it, this study established the target-compound-herb network to identify the significant target-compound.

Materials and Methods
2.1. Data Selection. miRNA profiles in GSE74209 and GSE93883 datasets were downloaded from Gene Expression Omnibus (GEO). Gene expression data were from GSE7158 and GSE35956 datasets. GSE74209 dataset contained 6 osteoporosis and 6 healthy cases. GSE93883 dataset contained 12 osteoporosis and 6 healthy cases. GSE7158 dataset contained 12 osteoporosis and 14 healthy cases. GSE35956 dataset contained 5 osteoporosis and 5 healthy cases. miRNA profiles from GSE74209 and GSE93883 datasets and mRNA profiles from GSE7158 and GSE35956 datasets were, respectively, merged and normalized using "normalization" package in R. Differentially expressed miRNAs (DEM) and genes (DEGs) were identified using "limma" package in R.

Weighted Correlation Network Analysis (WGCNA).
After detecting the outliers, "wgcna" package in R was applied to construct a WGCNA coexpression network. Firstly, soft thresholding power β was selected to ensure scale-free distribution (sale free R 2 = 0:9). Secondly, the coexpression network was constructed based on topological overlap matrix (TOM). And the cluster dendrogram and the module colors were plotted to display the modules. For miRNA, the dissimilarity degree of module Eigengenes (MEDissThres) was set to 0.25. And the minModuleSize was set to 30. Modules with these parameters were merged into a new module. Pearson's correlation coefficients were analyzed to exhibit the correlation between clinical traits and the module. P < 0:01 was seen as significant.

Functional Enrichment Analysis.
DEMs were subjected to DIANA Tools (http://www.microrna.gr/) to perform KEGG pathway analysis. "KEGG" package in R was used to enrich DGE-related KEGG pathways. The signaling pathways were seen to be significant when P < 0:05.
2.4. miRNA-mRNA Network Analysis. To examine the relationship between DEMs and DEGs, "network" package in R was downloaded. A miRNA-gene interaction network was constructed using Cytoscape software. The target genes of miRNAs were predicted by miRWalk.

Receiver Operating Characteristic (ROC) Analysis for Key miRNAs and Genes and Correlation with Key Genes and
Crucial Osteoporosis-Related Genes. ROC analysis was conducted to assess the predictive performance of hub miRNAs and genes, and the predictive values of the area under the curves (AUCs) of ROC curves specific to hub miRNAs and genes were calculated. In addition, correlations with these key genes and vital osteoporosis-associated factors were calculated to evaluate functional roles of these genes in osteoporosis.
2.6. Development of Molecule Docking Model. The 2D structures of chemicals were downloaded from PubChem database. ChemBio3D Ultra 14.0 software was applied to draw the 3D structures of the chemicals. The structure with minimum free energy conformation was saved as a mol2 file. The protein of gene was searched in UniPort database, and the 3D structure of the target protein was downloaded from the Protein Data Bank (PDB) archive and transformed into 2D structure with PyMOL software and saved as a PDB file. Then, the molecular docking model was developed using the AutoDock Vina software. The model with the minimum free energy conformation was visualized with PyMOL software.

Identification of Osteoporosis-Associated miRNAs.
After combination and normalization of miRNAs from GSE74209 and GSE93883 datasets, 1,810 miRNAs in 18 osteoporosis and 12 healthy cases were obtained (Table 1). To screen out miRNAs correlated with osteoporosis, WGCNA was performed to construct a miRNA coexpression network. No outliers were detected in the 30 samples (Supplementary Figure 1A). The soft thresholding power β was selected as 12 for scale-free network (scale-free R 2 = 0:9 with relatively high mean connectivity) (Supplementary Figure 1B). Then, 5 modules were obtained after similar modules were merged with the settings MEDissThres = 0:25 and minModuleSize = 30 (Figure 1(a)). The correlation between the 5 modules and osteoporosis was shown in Figure 1(b), and the strongest negative correlation was observed between the turquoise module and osteoporosis (R = −0:68, P = 3e − 05). Figure 1(c) displays high correlation between module membership and gene significance in the turquoise module, suggesting that the miRNAs in this module were strongly associated with osteoporosis. Therefore, further analyses were performed with 55 miRNAs in turquoise module.
Further, to identify osteoporosis-associated miRNAs, we screened out the DEMs in 18 osteoporosis and 12 healthy cases. 132 DEMs including 56 upregulated and 76 downregulated miRNAs were identified in response to osteoporosis with the cutoff criteria adjusted P value < 0.01 and jlog 2FCj ≥ 1 (Figure 1(d)). The overlapped genes between 132 DEMs and 55 miRNAs in the turquoise module were selected as miRNAs associated with osteoporosis ( Figure 1(e)), and 17 miRNAs were identified.

Identification of Osteoporosis-Associated Genes.
In GSE7158 and GSE35956 dataset, 23,520 mRNAs were identified in 17 osteoporosis and 19 healthy samples ( Table 2). The 36 samples were then subjected to WGCNA analysis. The soft thresholding power β was selected as 15, and 9 modules were obtained (Figures 2(a)-2(c)). The darkolivegreen module contained 3,344 mRNAs that were the most strongly correlated with osteoporosis ( Figure 2(c)). 1,572 DEGs including 807 upregulated and 765 downregulated genes were identified in response to osteoporosis (Figure 2(d)). The 222 osteoporosis-associated shared genes were obtained after comparing the 3,344 genes and 1,572 DEGs (Figure 2(e)).

KEGG Pathway
Analysis and Construction of miRNA-Gene Interaction Network. Next, we explored the biological functions of the 17 osteoporosis-associated miRNAs with DIANA Tools. The 38 KEGG pathways were significantly enriched such as ErbB signaling pathway, Wnt signaling pathway, AMPK signaling pathway, TGF-beta signaling pathway, and MAPK signaling pathway (Supplementary Table 1). The 222 osteoporosis-associated genes were also subjected to KEGG pathway analysis with "KEGG" package in R. And 14 KEGG pathways were significantly enriched such as cysteine and methionine metabolism, Salmonella infection, focal adhesion, gap junction, pathogenic Escherichia coli infection, glycerolipid metabolism, PI3K-Akt signaling pathway, and p53 signaling pathway (Figure 3(a)).
To determine the predictive performance of hub miR-NAs and genes, the AUCs of ROC curves about hub miR-NAs and core genes were calculated as 0.907 and 0.827 (Figures 4(a) and 4(b)), suggesting that these key miRNAs and genes had predictive actions for osteoporosis. In addition, correlation coefficients with these key genes and vital osteoporosis-associated factors including RUNX2, CALCA, and BMP2 were presented to be a trend of correlation with osteoporosis (Figures 4(c)-4(e)); complementally, these correlations coefficients were not significant, caused by insufficient sample size, the presence of missing values, and the existence of outliers. Overall, these data show a demonstrative effect of these key miRNAs and genes for osteoporosis.
3.4. Screening of Potential Drugs. Next, we screened potential drugs for osteoporosis in TCMSP (Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform) based on the 6 hub genes. And 3 genes   BioMed Research International (NTF3, IFNAR2, and THBS1) were targeted by 3 natural compounds (kainic acid, uridine, and quercetin), which were the active ingredients of 192 herbs. Then, we developed a target-compound-herb network to visualize the relationship among the 3 genes, 3 compounds, and 192 herbs ( Figure 5). For example, kainic acid, the active ingredient in herb Abri Herba, could affect osteoporosis by targeting NTF3; uridine, the active ingredient in herbs Atractylodes, Cordyceps, and Semiaquilegia, could affect osteoporosis by targeting IFNAR2; quercetin, the active ingredient in herbs Ardisia japonica Herba, Folium Artemisiae Argyi, Anisi stellati fructus, and Ginkgo Semen could affect osteoporosis by targeting THBS1. To investi-gate the relationship between compounds and the targets, molecular docking was performed. Figure 6 displays the 3 best molecular docking models with minimum free energy. The 3 compounds were, respectively, buried in the pockets of the 3 proteins, suggesting that the 3 compounds could be potential drugs of osteoporosis targeting the 3 genes.

Discussion
WGCNA is widely used to detect the coexpressed genes correlated with certain phenotype by developing a correlation network. Currently, accumulating studies have performed   [14]. Using WGCNA, Hu et al. discovered 7 genes that affected the progression of osteoporosis [15]. In this study, we detected osteoporosis-associated miRNAs and genes using WGCNA and differential expression analyses. We identified 5 modules for miRNAs, of which the turquoise module was the most strongly correlated with osteoporosis. For genes, 9 modules were grouped, and the darkolivegreen module was the most relevant to osteoporosis. Subsequently, 17 miRNAs and 222 genes were obtained associated with osteoporosis after comparing WGCNA and differentially expression results. For miRNAs, KEGG pathway analysis indicated that the 17 miRNAs were mainly enriched in 38 signaling pathways. Meanwhile, the 222 genes were mainly enriched in 14 pathways. And 3 shared KEGG pathways including focal adhesion, PI3K-Akt signaling pathway, and gap junction were discovered. Multiple studies have demonstrated the roles of the 3 signaling pathways in osteoporosis. Focal adhesions (FAs) are macromolecules that mediate the interactions between cells and extracellular matrix [16]. Multiple studies have reported that FAs participated in cell differentiation, mobility, and angiogenesis. H. Xie et al. demonstrated that FA signaling ameliorates the progression of osteoporosis by promoting the formation of H-type vessels in bone mesenchymal stem cells [16]. Besides, Hu et al. displayed that enhanced vessels in bones were positively associated with the level of FA kinases [17]. Phosphoinositide 3-kinase-(PI3K-) protein kinase B (AKT) signaling pathway, involved in cell survival, apoptosis, and growth, has been indicated to play critical roles in regulating osteoporosis development [18]. For example, the activation of PI3K-Akt signaling pathway elevated the expression of osteogenic differentiation genes such as ALP and BMP2 to promote the growth and differentiation of osteoblast [19]. However, the inactivation of PI3K-Akt signaling attenuated bone resorption of osteoclasts [20]. Gap junction proteins were mainly identified in osteocytes, osteoblasts, and osteoclasts. Many investigations have indicated that gap junction proteins are critical to the conservation of bone structure and formation. Gap junctions were decreased in osteocytes in mouse osteoporosis model [21]. Connexin 43, a gap junction protein, is essential for the function of osteocyte and reduced in high glucoseinduced osteoporosis [22]. These investigations all indicated the inhibitory roles of these signaling pathways in osteoporosis progression.
Currently, estrogen and bisphosphonate have been widely used in the treatment of osteoporosis. However, their side effects greatly limited their application in osteoporosis. For example, the use of estrogen elevated the incidence rate of breast cancer and heart failure [28]. Meanwhile, bisphosphonate compounds led to the irritation in gastrointestinal tract [28]. Therefore, given rare adverse effects, TCHMs,   [29]. For uridine, Li et al. demonstrated that uridine triphosphate, a kind of nucleotides, could inhibit osteogenic differentiation and promote adipogenic differentiation [30]. There was no study about therapeutic potential of kainic acid and uridine in osteoporosis.   Therefore, further investigations are needed for elucidate the role of the two compounds in the treatment for osteoporosis.

Conclusion
In conclusion, our study detected 6 hub miRNAs and 6 hub genes associated with osteoporosis using comprehensive bioinformatics analysis including WCGNA, differential expression, and KEGG pathway analyses. Then, the miRNA-gene interaction network was constructed. Predictive performances of hub miRNAs and genes were evaluated, and correlations with key genes and vital osteoporosis-associated factors were determined. And 3 natural compounds (kainic acid, uridine, and quercetin) were screened out based on the hub genes.
Our study sheds light on the pathogenesis of osteoporosis and provides promising therapeutic targets and drugs for osteoporosis. Some limitations of this study including correlations between targets and osteoporosis-related factors and specific mechanism of target-compound will be explored in the future work by larger samples and silencing/overexpressing the targets combined with compound administration.

Data Availability
The data that support the findings of this study are available in GEO at (https://www.ncbi.nlm.nih.gov/geo/).

Conflicts of Interest
There is no conflict of interests in this study to declare.

Authors' Contributions
Ran Han conceived and wrote the manuscript; Lubing Li designed the study and processed the data; Xiahatai Ayiding downloaded the data and processed the data, tables, and figures. All authors read and approved the manuscript.

Supplementary Materials
Supplementary Figure