WGCNA Identification of Genes and Pathways Involved in the Pathogenesis of Postmenopausal Osteoporosis.

PURPOSE: Postmenopausal osteoporosis (PMO) patients may suffer from chronic pain and increased fractures due to brittle bones that seriously affect their normal work and life. Exploring the pathogenesis of PMO can help clinicians construct individualized therapeutic targets.
METHODS: Differentially expressed genes (DEGs) were identified by analyzing the microarray assays of monocytes from 20 PMO and 20 control samples. Weighted correlation network analysis (WGCNA) and gene set enrichment analysis (GAEA) were performed. Genes associated with PMO were identified in the Comparative Toxicogenomics Database (CTD). miRNAs associated with osteoporosis were found in miRNet, and target genes were predicted. Hub genes and functional pathways associated with PMO were also identified. miRNA-mRNA networks were constructed. The association between hub genes and PMO was analyzed in the CTD.
RESULTS: A total of 1055 genes were up-regulated, and 694 genes were down-regulated in PMO samples (P<0.01). Five modules were identified by WGCNA. The blue module was significantly associated with PMO and selected for further analysis (P < 0.05). A total of 229 genes were significantly associated with PMO gene significance and module membership. Pathway variations were predominantly enriched in mRNA metabolic process, RNA splicing, Notch signaling pathway, apoptosis, cytokine-cytokine receptor interaction and so on. We identified 10 hub genes associated with PMO with different inference scores.
CONCLUSION: We identified genes, miRNAs, and pathways associated with PMO. These molecules may participate in the pathogenesis of PMO and serve as therapeutic targets.

Introduction

Many postmenopausal women suffer from postmenopausal osteoporosis (PMO). A survey of 3247 Italian postmenopausal women found that according to bone mineral density (BMD) diagnostic criteria, the prevalence of osteoporosis was 36.6%.1 PMO patients may suffer from chronic pain and fractures, and as a result their quality of life is seriously affected.2 With the aging of the population, osteoporosis will undoubtedly greatly increase the medical burden.3 The pathogenesis of osteoporosis includes inflammation, apoptosis, autophagy, and oxidative stress.4,5 Recent studies have shown that genes and miRNAs play an important role in the onset and development of PMO.6

Currently, bisphosphonates, calcium, and vitamin D are recommended for the prevention and treatment of osteoporosis. However, complications such as fractures may still occur. In addition, there are differences in therapeutic effects among different populations.7–9 Therefore, further exploration of the pathogenesis of PMO can help clinicians create individualized therapeutic targets.

Bioinformatics analysis is widely used to explore the molecular mechanisms and therapeutic targets of different diseases. Qiu et al found several genes and pathways involved in the pathogenesis of fibromyalgia, and provided an idea for the study of the molecular mechanism of fibromyalgia.10 Rodrigo identified several genes and miRNAs involved in the progression of osteoarthritis through cartilage sequencing analysis.11 Zhou identified several functional pathways involved in the regulation of BMD, such as the MAPK signal pathway, by sequencing and analyzing the peripheral blood monocytes of premenopausal and postmenopausal women with osteoporosis.12 Overall, there are still few reports on the use of bioinformatics analysis to explore the pathogenesis of PMO.

In this study, we identified differentially expressed genes (DEGs) by analyzing the microarray data of peripheral blood monocytes from PMO and normal postmenopausal women and then performed weighted correlation network analysis (WGCNA) and Gene set enrichment analysis (GSEA). We also identified hub genes and functional pathways related to the pathogenesis of PMO. Finally, we analyzed the association between hub genes and PMO in the Comparative Toxicogenomics Database (CTD).

Materials and Methods

Postmenopausal Osteoporosis Expression File

We downloaded dataset [GSE56815 (GPL96 platform)] from the GEO database (), an open-source platform for retrieving gene expression data.13 Gene expression of monocytes from 20 postmenopausal high bone mineral density (BMD) women and 20 postmenopausal low BMD women were then selected and analyzed.12

DEGs Identification and Principal Component Analysis

DEGs were identified by GEO2R, based on the limma packages of the R language. The cut-off was P-value<0.01. We visualized volcano diagrams with SangerBox software based on the R language as well () and performed principal component analysis (PCA) based on the expression level in ImageGP, an online tool also based on the R language ().

Weighted Gene Co-Expression Network Analysis

Weighted correlation network analysis (WGCNA) is a powerful method to identify significant clusters and module eigengene.14 WGCNA is widely used to find modules related to diseases and to explore their pathogenesis. To carry out this analysis, we constructed gene clustering of 20 PMO samples and 20 normal controls and identified the soft threshold according to scale independence and mean connectivity. Briefly, the soft threshold was 7. Next, we created the Topological Overlap Matrix (TOM) for clustering and constructed the module trait heat map according to the sample status. We randomly selected 400 genes to construct a TOM heat map and constructed an eigengene adjacency heatmap. A scatter plot of the turquoise module was created according to gene significance for PMO. The genes most associated with the blue module and PMO traits were selected for further analysis (cutoff: gene significance 0.2, module membership 0.5).

GO and KEGG Analysis

Gene Ontology (GO) covers biological process (BP), cellular component (CC) and molecular function (MF). KEGG (Kyoto Encyclopedia of Genes and Genomes) is widely used. To perform the GO and KEGG analysis of DEGs, we used the Database for Annotation, Visualization and Integrated Discovery (DAVID) () online tool. With a threshold for statistical significance of P<0.05. Pathways were visualized in ImageGP ().

Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) is a useful tool for interpreting genome-wide expression profiles.15 KEGG enrichment pathways of the GSE56815 were identified by GSEA analysis based on their expression levels. The parameters used were num: 100, scoring scheme: weighted, set max: 500, and set min: 15. Gene sets enriched at nominal P-value < 5% and FDR < 25% were regarded as significant.

Protein Network Construction and Hub Gene Selection

The Search Tool for the Retrieval of Interacting Genes (STRING) () was used to construct a (protein-protein interaction) PPI network of the identified DEGs and Cytoscape visualization software version 3.6.1 was used to present the network. Then, the most significant module in the PPI network was identified by Molecular Complex Detection tool (MCODE). The criteria were that scores >2, degree of cutoff=2, and node score cutoff=0.2.

The Comparative Toxicogenomics Database (CTD) can effectively predict the correlation between diseases, drugs and genes ().16 Hence, we found genes associated with PMO in the CTD. Then, we found osteoporosis associated miRNAs in miRNet (). At the same time, genes were predicted by osteoporosis associated miRNAs in miRNet. Common genes between MCODE genes, PMO associated genes in the CTD and the miRNAs predicted genes were then identified. The intersection of MCODE genes and postmenopausal osteoporosis associated genes in the CTD were identified as hub genes. Pearson clustering heatmap of hub genes were performed in ImageGP.

Hub Gene Interaction with Diseases

The inference score and reference count of hub genes associated with postmenopausal osteoporosis in the CTD were isolated. Interaction between hub genes and Musculoskeletal Diseases, Arthralgia, Cartilage Diseases, Osteoporosis, Metabolic Bone Diseases, Bone Resorption were analyzed in CTD. The inference score was visualized by the Radar chart.

miRNA-mRNA Network Construction

We constructed the miRNA-mRNA network was constructed in Cytoscape.

Results

DEGs and PCA

Our volcano diagrams show the DEGs (Figure 1A), and a total of 1749 DEGs were identified (P<0.01). We find that 1055 genes were up-regulated, and 694 genes were down-regulated in PMO samples. PCA results are shown in Figure 1B.

Figure 1

DEGs and PCA results. (A) DEGs between postmenopausal osteoporosis and healthy control samples; (B) PCA showed no significant batch effect.

WGCNA

Sample dendrogram and trait heatmap are shown in Figure 2A. The cluster was good. The scale independence and mean connectivity are shown in Figure 2B, with 7 being the perfect powerEstimate. We present the cluster dendrogram in Figure 2C, and find that 5 modules were identified. The network heatmap of 400 randomly selected genes is shown in Figure 2D.

Figure 2

Weighted gene co-expression network analysis results. (A) Sample dendrogram and trait heatmap; (B) Scale Independence and mean connectivity; (C) Cluster dendrogram; (D) Network heatmap of 400 randomly selected genes; (E) Module trait relationships; (F) Scatter plot genes in blue module.

The module trait relationships are shown in Figure 2E. The blue module was significantly associated with PMO and selected for further analysis (P<0.05). The scatter plot is shown in Figure 2F. A total of 229 genes were significantly associated with PMO gene significance and module membership.

DAVID Analysis and Gene Set Enrichment Analysis

We find that variations were predominantly enriched in mRNA metabolic process, RNA splicing, translational elongation, posttranscriptional regulation of gene expression, ribonucleoprotein complex, nucleotide binding, Insulin signaling pathway, Neurotrophin signaling pathway, Notch signaling pathway, Apoptosis and so on (Figure 3A–D). Detail information were shown in . GSEA showed that 67/166 gene sets were upregulated in phenotype Postmenopausal-osteoporosis. KEGG_PHENYLALANINE_METABOLISM, KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION, KEGG_NEUROACTIVE_LIGAND_RECEPTOR_INTERACTION, KEGG_HEDGEHOG_SIGNALING_PATHWAY, KEGG_STEROID_HORMONE_BIOSYNTHESIS, KEGG_MATURITY_ONSET_DIABETES_OF_THE_YOUNG were positive correlated with PMO (Figure 4A–F). Detail information for the results is shown in Table 1.

Figure 3

The enrichment analysis of the selected module genes. (A) Biological process; (B) Molecular function; (C) Cell component; (D) Kyoto Encyclopedia of Genes and Genomes (KEGG).

Figure 4

Gene Set Enrichment Analysis results. (A) KEGG_PHENYLALANINE_METABOLISM; (B) KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION; (C) KEGG_NEUROACTIVE_LIGAND_RECEPTOR_INTERACTION; (D) KEGG_HEDGEHOG_SIGNALING_PATHWAY; (E) KEGG_STEROID_HORMONE_BIOSYNTHESIS; (F) KEGG_MATURITY_ONSET_DIABETES_OF_THE_YOUNG.

Table 1

Detail of 20 Significant KEGG Pathway Enriched in Postmenopausal Osteoporosis Analyzed by GSEA

Term	Size	NES	Nominal P value
KEGG_SYSTEMIC_LUPUS_ERYTHEMATOSUS	103	1.8648449	<0.001
KEGG_PHENYLALANINE_METABOLISM	16	1.8480902	0.003656307
KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION	230	1.5717553	<0.001
KEGG_HEMATOPOIETIC_CELL_LINEAGE	82	1.5558474	0.014787431
KEGG_NEUROACTIVE_LIGAND_RECEPTOR_INTERACTION	245	1.5462313	<0.001
KEGG_HEDGEHOG_SIGNALING_PATHWAY	45	1.5296569	0.010928961
KEGG_STEROID_HORMONE_BIOSYNTHESIS	41	1.4420882	0.054104477
KEGG_PRIMARY_IMMUNODEFICIENCY	34	1.426645	0.06498195
KEGG_OLFACTORY_TRANSDUCTION	66	1.4145727	0.032315977
KEGG_GLYCOSPHINGOLIPID_BIOSYNTHESIS_LACTO_AND_NEOLACTO_SERIES	23	1.3720969	0.09310987
KEGG_ECM_RECEPTOR_INTERACTION	78	1.3505274	0.048543688
KEGG_ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM	25	1.3486422	0.10614525
KEGG_ETHER_LIPID_METABOLISM	25	1.3264432	0.1160221
KEGG_MATURITY_ONSET_DIABETES_OF_THE_YOUNG	22	1.3038279	0.13424124
KEGG_TASTE_TRANSDUCTION	32	1.2932026	0.12544803
KEGG_ARACHIDONIC_ACID_METABOLISM	45	1.2865875	0.10795455
KEGG_GLYCEROLIPID_METABOLISM	40	1.2472463	0.15458015
KEGG_CELL_ADHESION_MOLECULES_CAMS	116	1.236038	0.09878683
KEGG_HYPERTROPHIC_CARDIOMYOPATHY_HCM	73	1.2265385	0.14285715
KEGG_DRUG_METABOLISM_OTHER_ENZYMES	37	1.2201568	0.17333333

Abbreviations: KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, gene set enrichment analysis; NES, Normalized Enrichment Score.

Construction of a PPI network revealed 187 nodes and 426 edges, the average node degree was 4.56 (Figure 5A). A significant module was identified (Figure 5B). The CTD showed that 12,386 genes were associated with postmenopausal osteoporosis, and we find that 30 miRNAs were associated with osteoporosis in miRNet. Then, 5651 genes were predicted by the miRNAs. 10 hub genes (CCT5, SRP72, EIF5, RPS9, RPL27, CCT4, EXOSC2, RPL39, PSMA1, RPLP0) were identified (Figure 5C). Pearson correlation between hub genes was strong (Figure 5D). The information for the 10 hub genes in PMO is shown in Table 2.

Figure 5

Protein-protein interaction (PPI) network and hub gene identification. (A) PPI network; (B) MCODE genes; (C) Common genes between MCODE genes, postmenopausal osteoporosis associated genes in CTD and the miRNAs predicted genes; (D) Pearson clustering heatmap of hub genes.

Table 2

The Information of 10 Hub Genes in Postmenopausal Osteoporosis

Gene Symbol	Title	LogFC	P value
CCT5	Chaperonin containing TCP1 subunit 5	0.2178175	0.000584
SRP72	Signal recognition particle 72	0.321385	0.00197
EIF5	Eukaryotic translation initiation factor 5	0.335804	0.0000425
RPS9	Ribosomal protein S9	0.1325635	0.00798
RPL27	Ribosomal protein L27	−0.1227125	0.00035
CCT4	Chaperonin containing TCP1 subunit 4	0.1370405	0.00968
EXOSC2	Exosome component 2	0.15086	0.00102
RPL39	Ribosomal protein L39	−0.09098	0.0018
PSMA1	Proteasome subunit alpha 1	0.125466	0.00244
RPLP0	Ribosomal protein lateral stalk subunit P0	0.1207395	0.00649

Abbreviation: FC, Fold change.

CTD Analysis

As mentioned above, we find that 10 hub genes were associated with PMO with different inference scores and reference counts (Table 3). The genes CCT5, SRP72, EIF5, RPS9, RPL27, and CCT4 interacted with Musculoskeletal Diseases, Arthralgia, Cartilage Diseases, Osteoporosis, Metabolic Bone Diseases, Bone Resorption with different inference scores (Figure 6A–F).

Table 3

The Hub Genes Associated with Postmenopausal Osteoporosis in CTD

Disease	Gene Symbol	Inference Score	Reference Count
PMO	CCT5	2.97	1
PMO	SRP72	7.75	2
PMO	EIF5	15.8	16
PMO	RPS9	7.72	2
PMO	RPL27	2.85	1
PMO	CCT4	3.14	1
PMO	EXOSC2	12.53	3
PMO	RPL39	8	2
PMO	PSMA1	3.6	1
PMO	RPLP0	2.69	1

Abbreviations: CTD, Comparative Toxicogenomics Database; PMO, postmenopausal osteoporosis.

Figure 6

Interaction between hub genes and Musculoskeletal Diseases, Arthralgia, Cartilage Diseases, Osteoporosis, Metabolic Bone Diseases, Bone Resorption. (A) CCT5, (B) SRP72, (C) EIF5, (D) RPS9, (E) RPL27, (F) CCT4.

Part of the miRNA-mRNA network identified. Hsa-mir-139-5p/CCT5; hsa-mir-342-3p/RPS9; hsa-mir-320a/RPL27; hsa-mir-34c-3p/CCT4; hsa-mir-205-3p/SRP72; hsa-mir-21-5p/EIF5; hsa-mir-216a-3p/EIF5; hsa-mir-27a-3p/EIF5 is shown in Figure 7.

Figure 7

miRNA-mRNA network of hub genes and predicted miRNAs.

Discussion

PMO can make women prone to fractures due to brittle bones, which seriously impairs their normal work and life.17 The long-term effect of treatment after fracture in patients with osteoporosis is still far from satisfactory.18 Indeed, prevention of the progression of osteoporosis and complications such as fractures in the early stage is crucial in treating women with PMO. The current therapeutic drugs can improve osteoporosis to some extent, but the potential side effects of these drugs and the limited effect of fracture prevention are both shortcomings of current treatments.7,19

Currently much research has been focused on exploration of the pathogenesis of PMO in order to develop specific or individualized targeted therapeutic drugs.20 In this study, we found DEGs by analyzing the microarray data of PMO and normal control samples and identified genes associated with PMO in the CTD. miRNAs associated with osteoporosis were found in miRNet and target genes were predicted. 10 hub genes (CCT5, SRP72, EIF5, RPS9, RPL27, CCT4, EXOSC2, RPL39, PSMA1, RPLP0) and significant pathways were identified. miRNA-mRNA networks were constructed. EIF5, SRP72, CCT5 and RPS9 were selected for interpretation.

Eukaryotic translation initiation factor 5 (EIF5) is mainly involved in regulating translation initiation factor activity, eukaryotic initiation factor eIF2 binding, formation of translation preinitiation complex, translational initiation.21 EIF5 may participate in the regulation of oxidative stress and apoptosis by regulating miR-297.22 In addition, Nicole found that EIF5 can affect the proliferation and apoptosis of colorectal cancer cells and then affect the prognosis.23 Similarly, signal recognition particle 72 (SRP72) is involved in the regulation of RNA binding, signal recognition particle binding, ribosome binding, SRP-dependent cotranslational protein targeting to membrane and responses to drugs.21 Remko found that SRP72 may improve radiosensitivity by regulating apoptosis.24 SRP72 may also play an important role in regulating the expression of many proteins.25 We find that EIF5 and SRP72 were abnormally expressed between PMO and control groups. (P<0.05) At the same time, EIF5 and SRP72 interacted with Musculoskeletal Diseases, Arthralgia, Cartilage Diseases, Osteoporosis, Metabolic Bone Diseases, Bone Resorption with different inference scores. In addition, EIF5 and SRP72 were associated with PMO with different inference score and reference count analyzed in CTD. EIF5 was a target gene of osteoporosis associated miRNA hsa-mir-21-5p, hsa-mir-216a-3p and hsa-mir-27a-3p. SRP72 was a target gene of osteoporosis associated miRNA hsa-mir-205-3p. We speculate that EIF5 and SRP72 may be involved in the pathogenesis of PMO by regulating inflammation, oxidative stress, apoptosis, but these related molecular mechanisms warrant additional investigation.

Chaperonin containing TCP1 subunit 5 (CCT5) is involved in the regulation of ATP binding, G-protein beta-subunit binding, beta-tubulin binding, toxin transport.21 Mariana found that CCT5 may regulate autophagy and participate in the accumulation of neuroproteins, thus promoting the onset and development of Alzheimer’s disease.26 Wang demonstrated that CCT5 could regulate Muscovy duck reovirus (MDRV)-induced apoptosis and cell cycle.27 Additionally, ribosomal protein 9 (RPS9) is involved in the regulation of structural constituent of ribosome, translation regulator activity, positive regulation of cell proliferation and positive regulation of translational fidelity.21 Cheng demonstrated that RPS9 can participate in the onset and development of osteosarcoma by regulating MAPK signaling pathway.28 Yang found that RPS9 may be involved in the progression of non-small cell lung cancer by regulating apoptosis and proliferation. In addition, RPS9 can also regulate the JAK-STAT and calcium signal pathways.29 We find that CCT5 and RPS9 were abnormally expressed between PMO and control groups. (P<0.05) CCT5 and RPS9 interacted with Musculoskeletal Diseases, Arthralgia, Cartilage Diseases, Osteoporosis, Metabolic Bone Diseases, Bone Resorption with different inference scores. In addition, CCT5 and RPS9 were associated with PMO with different inference score and reference count. CCT5 was a target gene of osteoporosis associated miRNA hsa-mir-139-5p. RPS9 was a target gene of osteoporosis associated miRNA hsa-mir-342-3p. We speculate that CCT5 and RPS9 may be involved in the pathogenesis of PMO by regulating oxidative stress and apoptosis, and believe that the specific mechanism is worthy of further exploration.

Despite the rigorous analysis of this study, there are still some limitations. First, there is a lack of animal experiments. Second, the sample size we used was not large enough to consider population heterogeneity adequately. Finally, our predicted genes and functional pathways could benefit from further study.

Conclusion

We identified genes, miRNAs, and pathways associated with PMO. These molecules may participate in the pathogenesis of PMO and serve as therapeutic targets.