Construction of an miRNA-Regulated Pathway Network Reveals Candidate Biomarkers for Postmenopausal Osteoporosis

We aimed to identify risk pathways for postmenopausal osteoporosis (PMOP) via establishing an microRNAs- (miRNA-) regulated pathway network (MRPN). Firstly, we identified differential pathways through calculating gene- and pathway-level statistics based on the accumulated normal samples using the individual pathway aberrance score (iPAS). Significant pathways based on differentially expressed genes (DEGs) using DAVID were extracted, followed by identifying the common pathways between iPAS and DAVID methods. Next, miRNAs prediction was implemented via calculating TargetScore values with precomputed input (log fold change (FC), TargetScan context score (TSCS), and probabilities of conserved targeting (PCT)). An MRPN construction was constructed using the common genes in the common pathways and the predicted miRNAs. Using false discovery rate (FDR) < 0.05, 279 differential pathways were identified. Using the criteria of FDR < 0.05 and |log⁡FC| ≥ 2, 39 DEGs were retrieved, and these DEGs were enriched in 64 significant pathways identified by DAVID. Overall, 27 pathways were the common ones between two methods. Importantly, MAPK signaling pathway and PI3K-Akt signaling pathway were the first and second significantly enriched ones, respectively. These 27 common pathways separated PMOP from controls with the accuracy of 0.912. MAPK signaling pathway and PI3K/Akt signaling pathway might play crucial roles in PMOP.


Introduction
Estrogen plays a crucial role in keeping healthy bones. Estrogen deficiency has significantly negative effect on bone cell functions [1]. Postmenopausal osteoporosis (PMOP) is believed to result directly from the decreased endogenous estrogen in menopausal women [2,3]. PMOP is characterized by the severe loss of bone mass in the vertebrae and long bones. Globally, thousands of women over 50 years of age are influenced by PMOP which causes an increase of economic and societal burden recently. Therefore, studying PMOPrelated molecular mechanisms is urgently needed.
In recent years, previous studies have demonstrated that there are powerful genetic effects on the etiology of osteoporosis (OP) in postmenopausal women [4][5][6]. For example, polymorphisms in a set of genes have been found to be related to PMOP women in China, including ESR2 [7], HLA-B [8], and OPG [9]. Nevertheless, the etiology of PMOP still remains limited. MicroRNAs (miRNAs), as a kind of small noncoding RNAs of 22-25 nucleotides, have been implicated to exert indispensable regulatory functions in gene expression in animal kingdoms [10]. miRNAs have also been reported to play important roles in the regulation of biological processes, such as cell apoptosis, proliferation, and differentiation [11,12]. Growing evidence has demonstrated that miRNAs have an important role in regulating bone mass. In recent years, the role of miRNAs in the progress of OP has gained interest. Numerous studies have demonstrated that the dysregulation of miRNAs is associated with OP. miRNA-2861 has been found to affect osteoblast differentiation, contributing to OP in the femur via its effects on osteoblasts [13]. Moreover, miRNA-34c has been implicated to enhance the osteoclastogenesis through regulating the genes involved in the Notch signaling pathway [14]. Moreover, Chen et al. [15] have demonstrated that miRNA-125b is significantly overexpressed in PMOP, and miRNA-125b may be a potentially  noninvasive biosignature for PMOP. Similarly, Xiao et al. [16] have suggested the role of miRNA-129-5p in osteoblast differentiation and bone homeostasis. These findings suggest that miRNAs have potential effects on the pathways in the progression of OP and can thereby affect the treatment of OP. However, most of the studies have focused on only one or a few miRNAs/pathways using limited samples. Of note, little is known about specific pathways regulated by miRNAs during the development of PMOP.
Hence, we aimed to detect the risk pathways regulated by PMOP-related miRNAs through establishing an miRNAregulated pathway network (MRPN). Then, we analyzed the topological characteristics of the MRPN, and we also determined how the miRNAs regulated PMOP risk pathways. The results of our study may provide a novel viewpoint with respect to the mechanism and treatment of PMOP.

Methods
The current analysis was comprised of the following stepsdata selection (microarray profile, pathway data): pathway identification through calculating gene-level statistics and pathway-level statistics based on the accumulated normal samples using individual pathway aberrance score (iPAS), finding differentially expressed genes (DEGs) and significant pathways based on DEGs using DAVID software, identification of the common pathways of these two methods, miRNAs prediction based on the common genes via calculating TargetScores, MRPN construction, and the topological properties for MRPN. The detailed information of each step was described in Figure 1.

Data Availability.
To reveal the molecular mechanisms of PMOP, microarray analyses of monocytes were performed using Affymetrix HG-133A arrays in 80 females, including 40 pre-and 40 postmenopausal subjects. These data were deposited in the E-GEOD-56815 of the ArrayExpress database [14]. Probe IDs owning concentrated expression level were transformed into human gene symbols. Duplicated genes in matrix were eliminated. Overall, we obtained 8450 genes.

Collection of Pathway Data.
All biological pathways of humans were retrieved from REACTOME database (http://www.reactome.org/) [17]. Then, several pathways were filtered out when the gene set size of these pathways was larger than 100. In addition, the pathways with the intersection of 0 between genes enriched in the pathways and 8450 genes identified in the present analysis were also discarded. Ultimately, 1078 pathways were reserved for further analysis.

Gene-Level Statistics.
Herein, accumulated premenopausal samples were used as a reference and the expression level of genes was computed by comparing one postmenopause sample with many accumulated premenopause samples according to robust multichip average [18]. In detail, we firstly normalized the genes in all the premenopause subjects for the reference one by one to the reference, and then we calculated the mean value as well as standard deviation (SD) of the expression level. For a single postmenopause sample, quantile normalization was conducted after integrating the individual postmenopause sample with all reference samples. Assuming the genes having several probes, gene expression value was computed via averaging probe expression value. Gene-level statistics of each gene in an individual postmenopause subject was standardized Computational and Mathematical Methods in Medicine 3 according to the mean value and SD of the reference genes.

Pathway-Level Statistics.
We used the average method to calculate the pathway activity. Specifically, the gene-level statistics of all genes in each pathway were measured and summarized, following by transforming the mean value of gene-level into the pathway statistics of this pathway. The pathway statistics was calculated using the following formula: In this formula, iPAS is the expression status of a pathway, represents the gene-level statistics of th gene, and denotes the count of genes enriched in this given pathway.

Identification of Differential Pathways.
With the goal of analyzing the changes of pathways in PMOP and control samples, Gitools [19] was utilized to construct the cluster heatmap of pathways. -test was used to measure the pathway statistics of each pathway in all samples of the two groups. All the premenopause samples for the reference were compared one by one with reference to output the null distribution of pathway statistics. On the basis of the comparison between an individual postmenopause sample and the null distribution, we obtained the raw value. Subsequently, the values were corrected for multiple hypothesis testing based on false discovery rate (FDR) control using Benjamini and Hochberg (BH). Pathways with a FDR value of less than 0.05 were considered significantly differential.
2.6. Identification of DEGs. In our study, LIMMA package [20] of R language was used to calculate the expression level of each gene between two groups to further identify DEGs. Next, all expression scores were converted into foldchanges (FC) with log 2 base (herein, this means log FC). The log FC for each gene was determined as log(postmenopausal) − log(premenopausal), and the distribution of the log FC value of each gene expression was investigated. log FC was on behalf of the differential expression degree. Then, we used multiple test to adjust the original values by means of Benjamini et al. [21] method using FDR. DEGs were finally extracted when the threshold was set at FDR < 0.05 and | log FC| ≥ 2.

Pathway Analyses for DEGs.
To further reveal the biological roles of DEGs, DAVID [22] was used to perform the pathway analysis using the Expression Analysis Systematic Explorer (EASE) test [23]. The EASE score was utilized to extract the significant categories. In our work, the pathways having < 0.05 were regarded to be significant pathways.

Comparison of iPAS and DAVID.
To evaluate whether the iPAS method was feasible, the method proposed in the current study was compared with traditional DAVID software. The common pathways between the two methods were identified. In an attempt to evaluate the classification capacity of the common pathways, the hierarchical clustering analysis was implemented for the common pathways using Gitools [19].

TargetScore Calculation and Prediction of Potential miRNAs.
In detail, TargetScan context score (TSCS) is a sequence-based score for the individual target site, which was computed using TargetScan [24]. Probabilities of conserved targeting (PCT) are the probability of conserved targeting for single target site [25]. TargetScore is a flexible package which takes log FC values and sequence scores as inputs to calculate the probability of each gene being the target of a given miRNA. The function "TargetScore" offers a convenient way to acquire TargetScore with precomputed input (log FC, TSCS, and PCT), and the range of TargetScore is from 0 to 1 [26]. The greater the TargetScore is, the higher the accuracy in detecting known targets is. Therefore, relying on the log FC, TSCS, and PCT values, we counted the TargetScore for the genes enriched in the significant pathways to further identify the potential miRNAs by means of Variational Bayesian-Gaussian Mixture Model (VB-GMM) [27]. The TargetScore distribution for validated and nonvalidated targets of all miRNA-mRNA interactions each having at least 1 validated targets was evaluated. In the current study, the predefined was determined as the threshold when little overlap between the two distributions when TargetScore = existed. Furthermore, the miRNA-regulated pathways were achieved using their target gene sets. Next, we constructed an MRPN consisting of miRNAs, genes, and their common pathways. Cytoscape 2.8.3 was used to visualize the MRPN, and the topological properties of the MRPN were analyzed based on the Network Analysis plugin [28].

Identification of Differential Pathways.
After the genelevel statistics of all genes was transformed into the pathway-level statistics value of each pathway, -test was used for measuring the pathway statistics of each pathway in control and PMOP samples. Moreover, FDR approach was applied to correct significance levels ( values) for multiple hypothesis testing. Using the criteria of FDR < 0.05, a total of 279 pathways were extracted between the two groups. Cluster analysis was conducted to explore the changes of these 279 pathways. The heatmap of the 279 pathways was shown in Figure 2.

Identifying DEGs.
Before pathway analysis, DEGs between the two groups were firstly identified. Using the criteria of FDR < 0.05 and | log FC| ≥ 2, a total of 39 DEGs were retrieved. Specific information about DEGs was exhibited in Table 1. The most 5 significant DEGs were PTGES2, HAB1, SCT, ASNSD1, and TUT1.

Comparison of iPAS and DAVID.
To evaluate whether the iPAS method was feasible, the method proposed in Computational and Mathematical Methods in Medicine 5 the current study was compared with traditional DAVID software. It was observed that a total of 27 pathways were the common ones obtained from the iPAS approach and DAVID software. These common pathways were shown in Table 2. More importantly, hsa0410 (MAPK signaling pathway) was the first significantly enriched pathway. Of note, the pathway quantity selected from iPAS was increased, relative to the DAVID method (a total of 279 differential pathways were obtained by means of iPAS method, while only 64 significant pathways were extracted based on the traditional DAVID). Thus, our study demonstrated that the iPAS method showed more efficiency to identify significant pathways, compared with DAVID. The clustering analysis was carried out to evaluate the classification capacity for microarray data samples using the common 27 pathways detected by iPAS and DAVID method (Figure 3). The classification efficiency for microarray data samples was computed based on accuracy. These pathways separated PMOP from controls with the accuracy of 0.912. In light of this, we further infer that these pathways could be used to be marker pathways for PMOP with a high accuracy.

Construction of MRPN for the Common Pathways Based on the Common Genes of miRNAs and Pathways.
Tar-getScores were firstly calculated based on the combination of log FC values, and TSCS as well as PCT scores. Then, the distribution of TargetScores was obtained. According to the distribution of TargetScores, little overlap between the two distributions at TargetScore = 0.4 was found. Hence, TargetScore > 0.4 was selected as the threshold to filter out less confidence miRNAs. Using the TargetScore value > 0.4, a total of 39 genes, 50 miRNAs, and 145 miRNA-gene interactions were identified, and specific relationship was displayed in Figure 4.
We then constructed an MRPN for PMOP, which was shown in Figure 5. In this network, the most significant common pathway (MAPK signaling pathway) was regulated by has-miRNA-9. Furthermore, has-miRNA-34b-5p and has-miRNA-530f indirectly regulated 10 significantly differential pathways; among these 10 pathways, hsa-04151 (PI3K-Akt signaling pathway) was the second most significant common pathway.

Topological Features of the MRPN.
We determined the topological characteristics of the network, including the degree distribution. The degree distribution followed a power law distribution using all of the nodes in the MRPN. Based on the degree distribution, we selected 5 nodes as hub nodes, including NXT2, has-miR-1285, GMFB, has-miR-530f, and RPRD1A.

Discussion
In our study, a novel method of extracting pathway-based biomarkers in PMOP was developed, by investigating the interactions among genes, miRNAs, and pathways associated with pathogenesis through constructing an MRPN. The results of our work identified 279 differential pathways between PMOP and normal samples, including MAPK signaling pathway and PI3K/Akt signaling pathway.
MAPKs play crucial roles in cellular response to growth factors or cytokines [29]. The musculoskeletal system is frequently observed in several of the Ras/MAPK pathway syndromes, suggesting that activation of the Ras/MAPK pathway impacts cells regulating bone development and homeostasis [30]. OP is the systemic bone remodeling disease caused by imbalance between bone formation and bone resorption. Significantly, MAPK has been implicated to trigger the induction of c-Fos [31], and c-Fos is a key transcription factor of osteoclastogenesis [32]. Moreover, osteoclastogenesis exerts important functions in the development of OP [33]. MAPK signaling pathway has been demonstrated to be differentially activated in MSCs derived from osteoporotic postmenopausal women [34,35]. Therefore, we infer that MAPK is important for bone formation and may represent a target to prevent bone loss associated with PMOP. Of note, this pathway was regulated by has-miRNA-9. Has-miRNA-9 was a proinflammatory miRNA [36]. A former study has demonstrated that has-miRNA-9 plays important roles in suppressing tumor growth via mediating the downstream p38 MAPK pathway [37]. Accordingly, has-miRNA-9 might play important roles in the PMOP progression, partially through regulating the MAPK signaling pathway.
Previous studies have suggested that PI3K-Akt signaling pathway is implicated in the regulation of cell growth, proliferation, and survival [38]. Osteoblastic bone formation is stimulated through activating the PI3K-Akt signaling pathway [39]. Bone remodeling is an essential physiological process regulating bone mass and strength; an imbalance in bone resorption and formation results in bone loss during aging and in OP. Moreover, PI3K/Akt pathway is very important   Computational and Mathematical Methods in Medicine in osteoclastogenesis and bone mass in mice, which may offer a potential target in OP [40]. Remarkably, osteoclastogenesis exerts important functions in OP progression [33]. Has-miRNA-530f as a hub node in the MDPN regulated this pathway of PI3K/Akt signaling pathway in our study. According to the literatures, no studies have reported the association of has-miRNA-530f and PMOP. However, based on our results, we speculate that PMOP might be related to has-miRNA-530f-regulated PI3K/Akt signaling pathway.

Conclusion
In conclusion, we comprehensively collected the PMOP risk pathways and the miRNAs of PMOP to characterize the molecular mechanism and treatment strategy for PMOP. We then constructed an MRPN to better reveal the relationship among genes, pathways, and miRNAs and to develop promising drug candidates for PMOP. Based on the results, differential pathways, including MAPK signaling pathway and PI3K/Akt signaling pathway, were successfully identified. Moreover, these pathways-related miRNAs may be involved in the pathogenic process of PMOP. The findings of our study may be applied clinically for the diagnosis and treatment of PMOP patients.

Conflicts of Interest
The author declares no conflicts of interest.