Pathway and network analysis of genes related to osteoporosis.

As a common degenerative disease, osteoporosis (OS) is characterized by reduced bone mass and microarchitectural deterioration of bone tissue. Both genetic and environmental factors are involved in OS development. To date, ~300 genes have been confirmed to be involved in the pathogenesis of OS, a large majority of which have been independently investigated. As OS is a polygenetic disease, a comprehensive analysis focusing on the biological functions and interactions of OS‑related genes would provide valuable information. In this study, OS related research deposited in PubMed was retrieved and genes related to OS were catalogued. Pathways with an enriched biological function for these genes were extracted, and the crosstalk between the enriched pathways was analyzed. A comprehensive network was constructed, and a minimal network was extracted using the Steiner minimal network algorithm. In this study, a total of 294 genes in were retrieved from PubMed. Biological processes found to be enriched included those related to bone metabolism and the immune system. In total, 58 pathways were enriched. Furthermore, the comprehensive network consisting of 3,943 nodes and 7,976 edges was constructed, among which 631 nodes and 2,581 edges contributed to the OS‑specific molecular network. In this network, in excess of 300 potential genes associated with OS and two modules were identified. Thus, this study provides a mechanistic insight into OS and suggests more than 300 potential OS‑related genes for future research.

Introduction

As a common degenerative disease, osteoporosis (OS) is characterized by reduced bone mass and microarchitectural deterioration of bone tissue (1), which often means patients with OS suffer from back pain and bone fractures (2). It is estimated that ~50% of women and 30% of men >50 years old will suffer a fracture due to osteoporosis globally (3).

The development of OS can be influenced by genetic and environmental factors (4,5). Genetic factors play a greater role than environmental factors. It is estimated that the heritability of vertebral volumetric bone mineral density ranges from 0.27 to 0.51 (6). To date, animal models, gene expression, genome-wide association studies (GWAS) and systems biology approaches have identified hundreds of genes found to contribute to the process of OS (7–9). This suggests that OS is a polygenetic disease; that is, numerous genes exerts an effect on the development of OS, and multiple genes contribute to the pathophysiology of OS (10). Thus far, the function and role of these genes in the process of OS have only been investigated independently. However, a systematic strategy to collectively analyze the function and connection of these potential OS-related genes is preferable.

Recently, Qin et al (11) applied multiple computational approaches to analyze OS-associated single nucleotide polymorphisms (SNPs) and genes identified in GWAS. This study focused on gene regulatory networks, and retrieved some transcription factors (including NFATC2 and MEF2C) and microRNAS (miRNAs; including miR-3658 and miR-345-5p) that bind to loci of SNPs. However, the use of SNPs alone to identify genes involved in OS could lead to some genes known to be involved in OS being unaccounted for in such an analysis. To overcome this, genes known to be associated with OS were collected from PubMed and biological enrichment analyses was conducted to detect any significant functional themes within these genes. The biochemical pathways associated with these genes were analyzed for interactions among the enriched pathways. Finally, an OS-specific network based on the human protein-protein interaction network was constructed. The current study aimed to promote an understanding of the molecular mechanism of OS and to identify potential OS-related genes for future research.

Materials and methods

Identification of OS-related genes

OS-related gene candidates were curated by retrieving human genetic association studies deposited in PubMed (http://www.ncbi.nlm.nih.gov/pubmed/). Referring to published studies (12,13), reports related to OS were queried using the terms ‘osteoporosis, postmenopausal’ (MeSH Terms) OR ‘osteoporosis’ (MeSH Terms) AND ‘polymorphism, genetic’ (MeSH Terms) OR ‘genotype’ (MeSH Terms) OR ‘alleles’ (MeSH Terms) NOT ‘neoplasms’ (MeSH Terms). Up to June 1st 2018, a total of 1,335 publications were retrieved using these search terms. At this point, every reference was reviewed to guarantee that the conclusion was consistent with its contents. Studies reporting a significant association of gene(s) with OS were included in the present study. Thus, any gene involved in the pathogenesis or drug response to OS was incorporated, including research from GWAS and biochemical studies.

Functional enrichment analysis of OS-related genes

Pathway and process enrichment analysis was carried out using the following ontology sources: Kyoto Encyclopedia of Genes and Genomes (KEGG; Release 87.0; http://www.genome.jp/kegg), Gene Ontology (GO) Biological Processes (Version 2018-08-09; http://geneontology.org/), Reactome Gene Sets (Version 64; http://reactome.org), Molecular Signatures Database (Version 6.2; http://software.broadinstitute.org/gsea/msigdb/index.jsp) and CORUM (Version 3.0; http://mips.gsf.de/genre/proj/corum/index.html) (14). All genes in the genome were used as the enrichment background. In the gene set enrichment analysis, P-values were calculated based on cumulative hypergeometric distribution (15) and Q-values were calculated using the Benjamini-Hochberg procedure for multiple testing (16). The term was regarded as over-represented when P<0.05. To reduce redundancy in the ontology terms, the terms with P<0.01, a minimum count of 3 and an enrichment factor >2.0 were collected and grouped into clusters based on their membership similarities. In the process of hierarchical clustering of the enriched terms, k scores were used as the similarity metric, and sub-trees with similarity >0.3 were considered to be a cluster. The most statistically significant term within a cluster was chosen to represent the cluster. All of these results were obtained using the web-based tool, Metascape (Version 3.0, http://metascape.org) (14,17).

Functional enrichment analysis of OS-related genes and pathway crosstalk analysis

After enriching for biological process, the biochemical pathways of OS-related genes were further enriched using the web-based tool ToppGene (https://toppgene.cchmc.org/). In brief, the UniProt ID (https://www.uniprot.org/) list corresponding to the retrieved genes was uploaded into the server of ToppGene. Then, the list was automatically matched to canonical pathways using the KEGG and BioCarta (www.biocarta.com) pathway databases. A pathway was catalogued when it overlapped with one or more of the genes in the list. The overlap significance between the pathway and the input genes was analyzed using Fisher's exact test and the Bonferroni correction. The pathways were considered to be significantly enriched when the false discovery rate value was <0.05.

Pathway crosstalk was evaluated using the Jaccard coefficient (JC) and overlap coefficient (OC). Both of these describe the overlap between any given pair of pathways, and are defined as follows: where A and B are the number of genes included in the two tested pathways.

As pathways with too few genes may have insufficient biological information, only enriched pathways containing more than three candidate genes were included and pathway pairs with fewer than two overlapping genes were removed. All pathway pairs were ranked according to the average of the JC and OC coefficients. For clarity, crosstalk between pathways was also shown in the Cytoscape software (version 3.7.1; http://cytoscape.org).

Construction of the human interactome and an OS-specific subnetwork

To further investigate the interaction and correlation between these genes, a comprehensive human interactome was constructed based on the InnateDB database (http://www.innatedb.com). In this process, a novel web-based tool, OmicsNet (https://www.omicsnet.ca), was used. OmicsNet allows users to create molecular interaction networks and visually explore them in a three-dimensional space (18). The list of retrieved genes was uploaded to OmicsNet. Under the guidance of OmicsNet, an interactome was constructed based on the InnateDB database and the topological characteristics of potential molecular networks were analyzed using OmicsNet.

Next, an OS-specific network was extracted. The process was similar to the Steiner Tree problem, where the algorithm identifies a minimal sub-network containing all the terminal nodes from the complete network. OmicsNet was used for extraction, which implemented a heuristic approach that provides an approximate answer to this problem in order to reduce computation time.

To assess the non-randomness of the constructed network, 1000 random networks with the same number of vertices and interactions as the OS-specific network were generated using the Erdos-Renyi model (19) in the igraph R package (Version 0.7.1; http://igraph.org/) (20). The arithmetic average values of the shortest path distance and clustering coefficient were calculated. The number of random networks whose average shortest distance was less than that of the OS specific network was calculated and defined as ND. In a similar manner, the number of random networks whose average clustering coefficients was higher than that of the OS specific network was calculated and defined as NC. At last, the empirical P-values were calculated separately using ND/1000 and NC/1000. When both P-values were <0.05, the OS specific network was regarded as non-random.

Finally, potential modules were attained using the Infomap algorithm embedded in OmicsNet (https://www.omicsnet.ca) (18). Based on the principles of information theory, the Infomap algorithm transforms the problem of finding the potential modules into the problem of finding a description of minimum information for a random walk on the network graph (21). Using this strategy, the InfoMap algorithm joins neighboring nodes into modules.

Results

Identification of genes reported to be associated with OS

Up to June 1st 2018, a total of 1,335 publications were retrieved for OS. A total of 294 genes were retrieved (Table SI) and curated into a gene set (OS-related genes gene set, OSgset). Referring to a previously published article (9), this gene set can be divided into four subgroups: Calciotropic hormones and receptors (including CASR, CRHR and CTR), cytokines, growth factors and receptors (including BMP2, FGFR1 and IL-23), bone matrix proteins (including COL1A1, COL1A2 and ITGA1) and miscellaneous (including ADCY10, ALOX15 and ALOX5). The diversity of the genes significantly associated with OS clearly demonstrates that OS is a polygenetic disease.

Biological functions enriched in the OSgset

Functional enrichment analysis can reveal a more specific function of these genes. In total, >2,000 GO terms were significantly enriched in the genes analyzed. Considering the heavy overlap between ontology terms and the need to reduce redundancy, terms were collected with P-values <0.01, a minimum count of 3 and an enrichment factor >2.0, and grouped into clusters based on their membership similarity. The top 20 ranked clusters are shown as a heatmap (Fig. 1 and Table SII). Among these clusters, some biological processes can be discerned: These include the bone forming process (including ‘ossification’, ‘skeletal system development’, ‘tissue remodeling’, ‘tissue morphogenesis’, ‘heart development’, ‘negative regulation of cell differentiation’, ‘vasculature development’ and ‘regulation of secretion’), the bone regulation process (including ‘response to nutrient levels’, ‘response to steroid hormone’, ‘response to peptide’, ‘inflammatory response’, ‘cellular response to organic cyclic compound regulation of animal organ morphogenesis’ and ‘response to growth factor’), and signaling pathways in bone cells (including ‘regulation of MAPK cascade’, ‘GPCR ligand binding’, ‘transmembrane receptor protein tyrosine kinase signaling pathway’, ‘regulation of signaling receptor activity’ and ‘negative regulation of cell proliferation’). These results indicated that the candidate genes collected were reliable for follow-up bioinformatics analysis.

Figure 1.

Enriched biological processes. Biological enrichment was carried out with the following ontology resources: KEGG Pathway, GO Biological Processes, Reactome Gene Sets, Canonical Pathways and CORUM. Terms with P<0.01, a minimum count of 3 and an enrichment factor >2 were collected and grouped into clusters based on their membership similarities. Sub-trees with similarities >0.3 were considered to be a cluster. The most statistically significant term within a cluster was chosen to represent the cluster. Analysis was carried out Metascape. The x-axis denotes -log10(P) values based on the cumulative hypergeometric distribution. The colors denote the relative value of -log10(P): darker colors indicate a greater value of -log10(P).

Pathway enrichment analysis in OSgset

Enriching the biochemical pathways in which candidate genes are involved and analyzing their cross-talk can promote our understanding of the molecular mechanisms underlying OS. As shown in Table I (further details are provided in Table SIII), 58 significant enrichment pathways for OS were identified. Among them, the pathway ‘ensemble of genes encoding extracellular matrix and extracellular matrix-associated proteins’ was the most significantly enriched (P=3.60×10−18). In addition, some extracellular matrix-related pathways were also enriched, including ‘ensemble of genes encoding ECM-associated proteins including ECM-affiliated proteins, ECM regulators and secreted factors’, ‘genes encoding structural ECM glycoproteins’, ‘extracellular matrix organization’ and ‘ensemble of genes encoding core extracellular matrix including ECM glycoproteins, collagens and proteoglycans’. These findings confirmed a role for metabolism of the extracellular matrix in the development of OS.

Table I.

Pathways enriched in the OSgset.

Pathway	P-value[a]	pBH value[b]
Ensemble of genes encoding extracellular matrix and extracellular matrix-associated proteins	3.60×10−18	6.05×10−15
Wnt signaling pathway	3.71×10−13	6.23×10−10
Breast cancer	4.29×10−13	7.21×10−10
Genes encoding secreted soluble factors	4.63×10−13	7.78×10−10
Cytokine-cytokine receptor interaction	7.86×10−13	1.32×10−9
Ensemble of genes encoding ECM-associated proteins including ECM-affiliated proteins, ECM regulators and secreted factors	7.91×10−12	1.33×10−8
Canonical Wnt signaling	8.61×10−12	1.45×10−8
Class B/2 (Secretin family receptors)	1.10×10−11	1.85×10−8
GPCR ligand binding	1.14×10−11	1.91×10−8
Interleukin-4 and 13 signaling	3.51×10−11	5.90×10−8
Genes related to Wnt-mediated signal transduction	5.88×10−11	9.88×10−8
Proteoglycans in cancer	7.02×10−11	1.18×10−7
Cytokine Signaling in Immune system	1.52×10−10	2.56×10−7
Wnt signaling network	1.69×10−10	2.83×10−7
Cytokine Network	2.04×10−10	3.42×10−7
Negative regulation of TCF-dependent signaling by WNT ligand antagonists	2.35×10−10	3.95×10−7
Rheumatoid arthritis	7.22×10−10	1.21×10−6
Signaling by Wnt	1.07×10−9	1.80×10−6
Fluid shear stress and atherosclerosis	1.43×10−9	2.40×10−6
Basal cell carcinoma	1.49×10−9	2.51×10−6
Pathways in cancer	1.81×10−9	3.05×10−6
Alzheimer disease-presenilin pathway	1.83×10−8	3.08×10−5
Validated transcriptional targets of AP1 family members Fra1 and Fra2	2.84×10−8	4.78×10−5
Signaling pathways regulating pluripotency of stem cells	4.81×10−8	8.09×10−5
ALK in cardiac myocytes	6.37×10−8	1.07×10−4
Malaria	7.09×10−8	1.19×10−4
Cytokines and Inflammatory Response	1.21×10−7	2.03×10−4
Neuroactive ligand-receptor interaction	1.45×10−7	2.44×10−4
HTLV–I infection	1.54×10−7	2.59×10−4
Genes encoding structural ECM glycoproteins	2.30×10−7	3.86×10−4
Metabolism of steroid hormones	2.79×10−7	4.69×10−4
Signaling by Interleukins	4.17×10−7	7.02×10−4
Extracellular matrix organization	4.82×10−7	8.11×10−4
Ensemble of genes encoding core extracellular matrix including ECM glycoproteins, collagens and proteoglycans	5.09×10−7	8.56×10−4
Signaling by WNT in cancer	5.91×10−7	9.93×10−4
Wnt signaling pathway	7.16×10−7	1.20×10−3
Angiogenesis	8.07×10−7	1.36×10−3
Ovarian steroidogenesis	1.00×10−6	1.69×10−3
Hippo signaling pathway	1.13×10−6	1.90×10−3
Inflammatory bowel disease (IBD)	1.14×10−6	1.92×10−3
WNT ligand biogenesis and trafficking	1.22×10−6	2.05×10−3
AGE-RAGE signaling pathway in diabetic complications	1.36×10−6	2.28×10−3
G alpha (s) signaling events	3.12×10−6	5.24×10−3
Disassembly of the destruction complex and recruitment of AXIN to the membrane	4.21×10−6	7.07×10−3
Osteoclast differentiation	4.42×10−6	7.44×10−3
mTOR signaling pathway	4.64×10−6	7.80×10−3
PI3K-Akt signaling pathway	4.74×10−6	7.98×10−3
TCF dependent signaling in response to WNT	4.92×10−6	8.27×10−3
Wnt/beta-catenin Pathway	6.48×10−6	1.09×10−2
Interleukin-10 signaling	8.72×10−6	1.47×10−2
Glucocorticoid receptor regulatory network	9.88×10−6	1.66×10−2
HIF-1 signaling pathway	1.07×10−5	1.80×10−2
The IGF-1 Receptor and Longevity	1.12×10−5	1.88×10−2
IL23-mediated signaling events	1.17×10−5	1.97×10−2
Detoxification of Reactive Oxygen Species	1.41×10−5	2.37×10−2
Allograft rejection	1.41×10−5	2.37×10−2
Circadian Clock	2.01×10−5	3.38×10−2
RNF mutants show enhanced WNT signaling and proliferation	2.43×10−5	4.08×10−2

P-values were calculated by Fisher's exact test

pBH values were adjusted by the Benjamini and Hochberg method. OSgset, osteoporosis-related genes gene set.

The second enriched pathway was ‘Wnt signaling pathway’. Moreover, the number of Wnt-related pathways was the highest, including ‘Wnt signaling pathway’, ‘canonical Wnt signaling’, ‘genes related to Wnt-mediated signal transduction’, ‘Wnt signaling network’, ‘negative regulation of TCF-dependent signaling by WNT ligand antagonists’, ‘signaling by Wnt’, ‘signaling by WNT in cancer’, ‘WNT ligand biogenesis and trafficking’, ‘TCF dependent signaling in response to WNT’, ‘Wnt/beta-catenin Pathway’ and ‘RNF mutants show enhanced WNT signaling and proliferation’. These results suggested that Wnt-related pathways play an important role in the development of OS.

The third enriched pathway was ‘breast cancer’, which suggested that OS is closely related to breast cancer. In addition, immune-associated biological processes, including ‘cytokines and inflammatory response’ and ‘cytokine network’, were also significantly enriched, suggesting that the immune system may be involved in the etiology and pathological process of OS.

Crosstalk among significantly enriched pathways

To understand how these enriched pathways interact with each other, we performed a pathway crosstalk analysis for the 58 enriched pathways. All 58 pathways contained three or more members in OSgset. Among the pathways, 706 pairs of pathways met the criterion for crosstalk analysis, that is, each pathway pair shared at least two genes with one or more other pathways.

In Fig. 2, the crosstalk network is arranged based on the closeness centrality of every pathway. The maximum value of 0.84 was for the pathway ‘ensemble of genes encoding extracellular matrix and extracellular matrix-associated proteins’. The second and third largest were for ‘genes encoding secreted soluble factors’ and ‘ensemble of genes encoding ECM-associated proteins including ECM-affiliated proteins, ECM regulators and secreted factors’. The results suggested that these three pathways lie in the central loci of the crosstalk network, and that abnormalities in the extracellular matrix plays an important role in the development of OS.

Figure 2.

Crosstalk map of enriched pathways. Pathway crosstalk was evaluated using the Jaccard coefficient and overlap coefficient and was displayed in the Cytoscape software. The network map was displayed according the closeness centrality in an ascending order. The minimum value was from the ‘metabolism of steroid hormones’ pathway. The maximum value was from the ‘ensemble of genes encoding extracellular matrix and extracellular matrix-associated proteins’ pathway.

Among this network, the minimum value was from the pathway of ‘metabolism of steroid hormones’ (including CYP17A1, CYP19A1, SRD5A2, HSD11B1, HSD17B2, AKR1B1, LHB and STARD3NL). Moreover, this pathway had only one connection with the ‘ovarian steroidogenesis’ pathway, while the other nodes have more edges with other pathway nodes. These findings suggest that steroid hormones affect the formation of OS in an indirect manner.

Network topological characteristics of OSgset

Analyzing the topological properties of nodes and interactions between nodes via protein-protein interaction-based analysis can help to reveal any key biologically related mechanisms associated with the network. The interactome constructed in this present study comprised 3,943 genes/proteins and 7,976 interactions, with a mean degree of 8.23. As specified, 259 of the 294 genes were included in the analysis and accounted for 88.1% of the genes in OSgset, serving as seeds (Table II).

Table II.

List of genes in different functional groups.

A, Genes in OSgset serving as seeds
ADRA2A, JAG1, AGTR2, AKR1B1, ALOX5, ALOX15, ALPL, ANXA6, APC, ARHGAP1, ATP5E, AVP, B2M, BDNF, BGLAP, BLVRB, BMP2, BMP7, BMPR1B, CA8, CALCR, CASR, CAT, RUNX2, CCNE1, CD40, CD40LG, CFTR, CHD2, CLCN7, CNR1, COL1A1, COL6A1, COL9A1, COMT, CPB2, CRHR1, CSTA, CTNNB1, CTSZ, CYLD, CYP1A1, CYP1B1, CYP17A1, CYP19A1, CYP24A1, DBP, DLX5, DMP1, ESR1, ESR2, ESRRA, MECOM, FABP3, FDPS, FGFR1, FGFR2, FOXC2, FLNB, FLT1, FRZB, GC, MSTN, GHR, GIP, GIPR, GLA, NR3C1, GSR, GSTM3, GSTP1, HIF1A, HSD11B1, HSD17B2, HTR2B, IBSP, IRF8, ID4, IFNG, IGF1, IGFBP2, IL1A, IL6, IL6R, IL7, IL10, IL15, IL16, IL17A, ITGA1, JUND, LCT, LEPR, LHCGR, LRP1, LRP4, LRP6, LRP5, LTA, LTBP2, MARK3, MECP2, MEF2C, MGP, CIITA, MID1, MMP2, CD200, MSX1, MYO5A, NELL1, NF1, NFATC1, NFKB2, NOS3, NPY, ROR2, TNFRSF11B, P2RX1, P2RX4, P2RX7, PCSK6, SERPINE1, REG3A, PAPPA, PBX1, ENPP1, PGR, ABCB1, PLAUR, PLOD1, PLS3, PON1, PPARG, PRL, PTCH1, PTH, PTHLH, PTH1R, PTH2R, PTMA, PTN, RAC1, RAC2, RB1, BRD2, RPL29, SAA1, ATXN1, CCL2, SFRP1, SFRP2, SFRP4, SHBG, SLC6A4, SNCA, SOD1, SOD2, SOX4, SOX9, SPARC, SPP2, SPP1, SPTBN1, SREBF1, TERT, TGFBR3, TIMP1, TIMP2, TLR4, TNFRSF1B, TRPS1, TSHR, TWIST1, TYROBP, VDR, VEGFA, WNT1, WNT7B, WNT10B, PRDM2, CSDE1, CXCR4, SHFM1, HMGA2, GDF5, FZD1, FZD6, SPOP, PIR, TNFSF11, NCOA1, TNFRSF11A, PER3, NOG, MAP4K4, ATP6V1G1, TRAM2, ZBTB40, NR1I3, HDAC5, GPC6, SLC25A13, SPRY1, ZMPSTE24, AKR1A1, TUBA1B, FAM3C, GPNMB, CRTAP, NCOA2, NMU, AKAP13, AKAP11, MMRN1, DKK1, CLEC16A, SATB2, DICER1, DCAF13, FGF20, DKK2, TMEM14A, CNOT7, NOX4, ARHGEF3, SOST, SLMO2, ADIPOR1, WNT16, IL23A, ATP6V1H, SOX6, ITLN1, ADCY10, MEPE, RGMA, SQRDL, TRPV4, SMOC1, IFIH1, TMEM135, WNK4, VKORC1, FTO, WLS, CCDC170, WNT5B, STARD3NL, NLRC5, WNT3A, SP7, FAM210A, IL23R, VPS13B, ZNF384, THSD7A, NCF1, CCR2
B, Genes included in OS-specific network but not in OSgest
UBC, HNF4A, SP1, ELAVL1, EP300, SUMO2, CXR4, RELA, CREBBP, SUMO1, STAT3, HSP90AA1, TP53, SMAD3, SRC, HDAC1, APP, SIRT1, JUN, EGR1, COPS5, TRAF6, SMAD2, GRB2, PIK3R1, YWHAZ, SMARCA4, MYC, HSP90AB1, NFKB1, STAT1, CALM1, HDAC3, HSPA8, HDAC4, E2F1, STUB1, CEBPB, YWHAQ, FOS, CUL3, HDAC2, RXRA, CAV1, CREB1, NCOR1, CEBPA, MDM2, KAT2B, GSK3B, SMAD4, CEBPD, IRF1, POU2F1, YY1, FN1, HSPA4. ATF2, A2M, KPNB1, HNF1A, UBE2I, XRCC5, NOS2, MAPK1, SUMO4, NCOA6, ZBTB16, VCP, KIAA0101, NCK1, TCF3, PLCG1, FYN, H2AFX, HDAC6, EPAS1, ACTB, CANX, HSPA5, EEF1A1, TUBB, TUBA1A, FLNA, SPI1, ETS1, NEDD4, AKT1, CCND1, PRKDC, EZH2, XRCC6, TBP, PPP2CA, XPO1, BCL6, PIAS1, REL, CRK, HNRNPA1, DNMT1, UBE3A, CBL, TRAF3, PRKCA, IL7R, TFAP2A, TRAF2, PSMC5, USF1, TNF, PPP1CA, HIST1H4E, MYOD1, STAT5A, VDAC1, PSMC4, MED24, BIRC2, NPM1, STAT5B, UBQLN1, CDK1, ITGB1, SLC9A3R1, FBXO6, SMAD1, PLK1, EIF3F, BCL3, ONECUT1, CAND1, PSMD1, PSMA7, AHR, DHX9, BAG6, ABL1, SQSTM1, DDB1, PHB, THBS1, SMAD9, PAK2., STAT6, NME1, RPL5, NPEPPS, CRYAB, DDIT3, PPP1CC, P4HB, NR2F1, APEX1, HSP90B1, GNAI1, CD4, CSNK2A1, UBQLN4, PCMT1, DNMT3B, RBX1, USF2, FOSL1, MCM3, POLR2F, PSMD11, LDHA, SKIL, ITGA5, BMPR1A, PITX3, NR5A1, DNAJB11, YWHAH, PSMC3, CBX5, CCND3, BCL2, CGA, LRIF1, KHDRBS1, CACNA1A, CDC73, NFKBIA, SNTA1, ARHGDIA, SRPK1, ENG, CCL7, RAD21, NR2F6, MPG, HAX1, TAB2, CALR, PCNA, STAU1. GATA1, BCCIP, CLU, SGTA, SKP2, TF, RELB, NFYA, RORA, FTH1, PSMA6, VDAC2, POLR2A, ERBB2IP, MYD88, PEPD, JAK1, PPP2R4, ATF1, GATA2, CHD3, ARRB1, TBK1, ELF1, PLG, NR1H3, GNB1, CTCF, YBX1, MAP3K7, ACVR2B, ETS2, IGFBP1, FZD8, SORD, CHEK2, LNX1, ATP5F1, ZAP70, ACTA1, MAG, PNO1, TGFB1, ACTR2, PPIB, FIS1, ERG, TGFBR2, GIPC1, ATP5A1, TERF2IP, TOLLIP, ACTA2, TRAM1, SGK1, MMP14, PRKCB, PROCR, IGFBP5, MSN, CEBPE, CSNK1E, KDM5B, RANBP1, STK39, CTTN, BARD1, ARRB2, APPL1, FLOT1, JAK3, CD14, FASLG, WNT4, RIPK2, CLNS1A, SRPK2, PITX2, MCM2, EGFR, CSNK2B, ATF7IP, EIF2C3, APOA1, RBM23, GEMIN4, BATF3, SLC25A1, CCR5, TXN, LYN, GNAI2, ALB, RNF4, SOX2, CDC37, TCEB1, INSR, FTL, CSNK1D, SYK, TBCA, CD63, BABAM1, NR1H2, TIMP3, IGFBP3, UGGT1, COL2A1, PSAT1, EIF4A2, UCHL1, SPEN, STAM, IGF2R, SGSM2, GNAS, IRAK4, ECH1, IL1R1, CD9, KLRD1, LIG4, SERPINA1, MSMO1, LRP2, CAMP, IL4R, GPRASP2, SHARPIN, REST, COL1A2, PPP1R16A, TLN1, AHSG, RAB3A, CUBN, IGSF1, SFTPD, SDC2, DCN, GNAQ, HERPUD1, APOB, HOXA1, IGFBP7, ITSN2, TANK, DDX39B, BBS10, GATA4, IDE, ITGB3, SNCG, CDH5, CREM, BMI1, DSTN, BCL2A1, CDKN2A, CTSA, CD36, CASK, MMP3, EGR2, NOTCH2, PRKACB, SURF2, DES, MYBPC2, IKZF3, RBPJ, SH3GL2, GLI1, APCDD1
C, Genes in the main module extracted from the OS specific network
A2M, AHR, AKT1, APC, BIRC2, ARRB1, ATF1, B2M, CCND1, BCL3, BCL6, CALM1, CAV1, RUNX2, CCNE1, CD14, CD40, CDK1, CEBPA, CEBPB, CEBPD, CFTR, CLU, CCR5, CREB1, ATF2, CREBBP, CREM, CSNK2A1, CTNNB1, CYP19A1, CYP24A1, DDB1, DHX9, E2F1, EEF1A1, EGFR, EGR1, ELF1, EP300, EPAS1, NR2F6, ESR1, ESR2, ESRRA, ETS1, ETS2, MECOM, FLNA, FOS, FTH1, NR5A1, XRCC6, GATA2, NR3C1, GSK3B, HDAC1, HDAC2, HIF1A, HNRNPA1, HSPA4, HSPA5, HSPA8, HSP90AA1, HSP90AB1, IRF8, IFNG, RBPJ, IL6, IL10, IL16, IRF1, JUN, JUND, KPNB1, LRP4, LTA, SMAD2, SMAD3, SMAD4, MARK3, MCM3, MDM2, MECP2, MEF2C, CIITA, MID1, MMP2, CD200, MYC, MYO5A, MYOD1, NF1, NFATC1, NFKB1, NFKB2, NME1, NOS2, NOS3, NPM1, P2RX1, P4HB, SERPINE1, PBX1, PCNA, PGR, PHB, POLR2A, POU2F1, PPARG, PPP2CA, PRKDC, MAPK1, PRL, PSMA6, PSMA7, PSMC3, PSMC4, PSMC5, PSMD1, PSMD11, PTHLH, PTH1R, PTMA, RB1, REL, RELA, RNF4, RXRA, ATXN1, CCL2, CCL7, SFTPD, SGK1, SKP2, SMARCA4, SUMO2, SNCA, SOD1, SOD2, SOX4, SOX9, SP1, SPI1, SPTBN1, SRC, SREBF1, STAT1, STAT5A, STAT6, TBP, TCF3, TERT, TFAP2A, TLR4, TNF, TNFRSF1B, TP53, HSP90B1, TRPS1, TWIST1, UBE2I, UBE3A, SUMO1, USF1, USF2, VCP, VDAC1, VDAC2, VDR, VEGFA, XPO1, XRCC5, YY1, YWHAH, YWHAZ, ZBTB16, PRDM2, TUBA1A, CXCR4, SHFM1, FOSL1, HIST1H4E, PIAS1, NCOA1, EIF3F, RIPK2, HDAC3, KAT2B, SLC9A3R1, NCOR1, HDAC4, MED24, NR1I3, RBX1, HDAC6, HDAC5, NR1H3, STUB1, TUBA1B, HAX1, NCOA2, CTCF, YWHAQ, COPS5, NCOA6, TAB2, SIRT1, FGF20, TBK1, NOX4, WNT16. SOX6, CAND1, UGGT1, UBQLN4, SQRDL. WNK4, TUBB

OS, osteoporosis; OSgset, OS-related genes gene set.

OS-specific molecular network inference

An OS-specific minimal subnetwork was extracted. As shown in Fig. 3, this extracted minimal subnetwork shares 631 nodes and 2,581 edges with the 259 seeds. The average shortest path distance was 3.12, which was significantly smaller than the mean shortest-path distance (5.13; empirical P<0.01). Similarly, the average clustering coefficient of the OS specific network was significantly higher than that of the random networks (0.052 vs. 0.013; empirical P<0.001). These results indicate that the extracted OS-specific network is a non-random network. Therefore, it is proposed that the remaining 372 genes (listed in Table II) are un-identified genes that may be in involved in the OS process.

Figure 3.

Osteoporosis (OS)-specific network map. The OS-specific network was extracted from the listed comprehensive human interactome network by solving the Steiner minimal tree problem using OmicsNet with 631 nodes and 2,581 edges. The network map was displayed according to the closeness centrality in an ascending order. The minimum closeness value was that of NMU and the maximum closeness value was that of UBC. OS, osteoporosis; NMU, neuromedin U; UBC, ubiquitin C.

The functional modules in the minimal subnetwork were also investigated, and two modules were identified. The first shares 219 genes (listed in the Table II). For this module, the first three annotations were for the positive regulation of cellular metabolic processes (P=1.69×10−70), the positive regulation of metabolic processes (P= 6.23×10−68) and the positive regulation of transcription from the RNA polymerase II promoter (P=1.51×10−65). The second consisted of 9 genes (AGTR2, GIP, GIPR, HTR2B, SNTA1, TIMP3, NMU, ERBB2IP and LNX1). This third gene annotations were for the regulation of secretion (P=5.43×10−4), the positive regulation of secretion (P=6.15×10−4) and regulation of transport (P=9.79×10−4), which indicated that the modules functions in the regulation of secretion.

Discussion

In this study, genes associated with the development of OS were collected and their biochemical pathway functions enriched. Using this, a gene-interacting network was constructed. In doing this, a gene-interacting framework was developed and >300 potential OS-related genes were identified for further research.

Identifying the biochemical pathways that the OS-related genes are involved in will provide an insight to promote the understanding of OS. In the 58 pathways identified, the top three pathways are involved in extracellular matrix metabolism, the cytokine and cytokine receptor network, and in Wnt signaling. Further to this, it was found that all 58 pathways are involved in cross-talk and constitute a module. These findings suggest that pathways function in a concerted manner and not independently, which is supported by other studies (22–24). For instance, Zhang et al (25) reported that lipopolysaccharide activates the mitogen-activated protein kinase/Wnt/NF-κB pathway in cultured chondrocytes, and subsequently induces apoptosis. This process can be rescued by treatment with calcitonin.

When the networks with cross-talk were arranged based on the closeness centrality of every pathway, it was found that the maximum closeness value was from the ‘ensemble of genes encoding extracellular matrix and extracellular matrix-associated proteins’ pathway. Thus, it is suggested that the aforementioned pathway lies in the central locus of the network. On the contrary, the minimum value was from the ‘metabolism of steroid hormones’ pathway. Thus, it was proposed that abnormalities in the extracellular matrix plays an important role in the development of OS and that steroid hormones affect OS development in an indirect manner.

In the 58 pathways identified, the pathway enriched as ‘breast cancer’ is of interest as it shares 22 genes with the OSgset. It has been noticed that survivors of breast cancer often suffer from OS (26). One explanation for this is that many survivors of breast cancer experience a loss of ovarian function and a drop in estrogen levels due to chemotherapy or surgery, which promotes the development of OS (27). Another explanation is that breast cancer cells can secrete γ-secretase, cyclooxygenase-2 and interleukin-8, promoting the process of osteoclastogenesis (28,29). Based on the results presented here, it is suggested that OS and breast cancer may share a similar genetic background.

Notably, >300 potential genes were found in the OS-specific network presented here. It is predicted that these genes are potential candidates for further research into the molecular mechanism of OS. For instance, the receptor tyrosine-protein kinase ERBB2 interacting protein (ERBB2IP) binds to unphosphorylated ERBB2 protein and regulates ERBB2 function and localization (30). In addition, it can also disrupt the Ras-Raf interaction and affect the Ras signaling pathway (31). As the Ras-Raf signaling pathway participates in bone metabolism (32), it is possible that ERBB2IP is involved in the development of OS. Thus, the extracted network outlines a portrait of OS, and additionally identified a number of genes for future study.

A module is defined as having more connections within than the average number of connections across the entire network, which indicates that the genes involved in a module share the same biological function. In this study, two modules were extracted from the OS-specific network. The first consisted of 219 genes directly involved in bone metabolism. However, the second module, which includes nine genes and is annotated as ‘regulation of secretion’, is more attractive. Among these genes, neuromedin U (NMU), gastric inhibitory polypeptide (GIP), GIPR, AGTR2 and HTR2B have been confirmed to be involved in the process of OS. It has been reported that GIP (a gut hormone) promotes bone anabolism as an entero-osseous hormone by stimulating osteoblast differentiation and increasing osteoblast longevity (33). In addition, GIP attenuates the activity of osteoclastic cells, leading to a net increase in bone deposition and ultimately increasing bone mass (33). NMU is an anorexigenic neuropeptide that acts independently of leptin through poorly defined mechanisms. Physiological and cell-based assays indicate that NMU acts in the central nervous system, rather than directly on bone cells, to regulate bone remodeling (34). It has also been reported that NMU promotes small intestinal transit, and that NMU deficiency results in a lowered intestinal motility rate and diminishes the effect of serotonin-induced defecation and diarrhea (35). Based on these facts, a ‘neuro-entero-osseous’ model is proposed: The proteins α1-syntrophin, 5-TH2B and ANGII affect the activity of neurons in the gastro-interstitial tissue, alter the secretion of GIP, and influence bone metabolism. Due to a lack of experimental and clinical data supporting the role of other genes in this extracted module, a detailed mode of action cannot yet be given. Therefore, more experimental studies are required in the future.

It is difficult to diagnose or treat OS based on the analysis of all these potential genes. However, it is hypothesized that further studies examining these genes would be of valuable. It is further hypothesized that these genes may serve as a diagnostic test for OS and as drug targets in the future.

At least two studies have applied bioinformatic and computational methods to study OS. Qin et al (11) investigated OS-associated SNPs and genes identified by GWAS, and found a number of SNPs that may influence the binding affinity of transcription factors (NFATC2, MEF2C, SOX9, RUNX2, ESR2, FOXA1 and STAT3) and miRNAs. In 2018, Sheng et al (36) proposed a computational workflow to curate and evaluate OS related genes, and they emphasized the important role of OS-related genes, including TGFB1, IL6, IL1B, TNF, ESR2, IGF1, HIF1A, COL1A1 and IFENG. In comparison with these studies, the work presented here focuses on the relationship between known genes and suggests a number of unknown genes that may be involved in the process of OS.

In summary, this study revealed the pathways of all currently known genes in OS and the crosstalk between these pathways. Based on this analysis, a mechanism of OS for proposed. Additionally, in excess of 300 genes not currently associated with the development of OS have been identified. Therefore, this study increases the understanding of the contribution of genetic factors to OS and identifies genes for further investigation.