Objectives In order to screen the altered gene expression profile in peripheral blood mononuclear cells of patients with osteoporosis, we performed an integrated analysis of the online microarray studies of osteoporosis.
Methods We searched the Gene Expression Omnibus (GEO) database for microarray studies of peripheral blood mononuclear cells in patients with osteoporosis. Subsequently, we integrated gene expression data sets from multiple microarray studies to obtain differentially expressed genes (DEGs) between patients with osteoporosis and normal controls. Gene function analysis was performed to uncover the functions of identified DEGs.
Results A total of three microarray studies were selected for integrated analysis. In all, 1125 genes were found to be signiﬁcantly differentially expressed between osteoporosis patients and normal controls, with 373 upregulated and 752 downregulated genes. Positive regulation of the cellular amino metabolic process (gene ontology (GO): 0033240, false discovery rate (FDR) = 1.00E + 00) was significantly enriched under the GO category for biological processes, while for molecular functions, flavin adenine dinucleotide binding (GO: 0050660, FDR = 3.66E-01) and androgen receptor binding (GO: 0050681, FDR = 6.35E-01) were significantly enriched. DEGs were enriched in many osteoporosis-related signalling pathways, including those of mitogen-activated protein kinase (MAPK) and calcium. Protein-protein interaction (PPI) network analysis showed that the significant hub proteins contained ubiquitin specific peptidase 9, X-linked (Degree = 99), ubiquitin specific peptidase 19 (Degree = 57) and ubiquitin conjugating enzyme E2 B (Degree = 57).
Conclusion Analysis of gene function of identified differentially expressed genes may expand our understanding of fundamental mechanisms leading to osteoporosis. Moreover, significantly enriched pathways, such as MAPK and calcium, may involve in osteoporosis through osteoblastic differentiation and bone formation.
Cite this article: J. J. Li, B. Q. Wang, Q. Fei, Y. Yang, D. Li. Identification of candidate genes in osteoporosis by integrated microarray analysis. Bone Joint Res 2016;5:594–601. DOI: 10.1302/2046-3758.512.BJR-2016-0073.R1.
- Integrated analysis
- Peripheral blood mononuclear cells
- Differentially expressed genes
To screen the altered gene expression profile in peripheral blood mononuclear cells of patients with osteoporosis.
To interpret the biological roles of the identified differentially expressed genes in patients with osteoporosis.
To identify the potential candidate genes in osteoporosis.
A total of 1125 genes were found to be signiﬁcantly differentially expressed in peripheral blood monocytes (PBMs) between osteoporosis patients and normal controls, with 373 upregulated and 752 downregulated genes.
Some genes, such as nuclear receptor interacting protein 1, interleukin 1 receptor associated kinase 3, and ubiquitin specific peptidase 9 X-linked, were highly correlated with the development of osteoporosis.
The significantly enriched pathways, such as the mitogen-activated protein kinase and the calcium signalling pathways, may involve in osteoporosis through osteoblast differentiation and bone formation.
Strengths and limitations
This study revealed the gene expression profiles in PBMs between osteoporosis patients and normal controls with large statistical power.
The gene function analysis of identified differentially expressed genes may expand our understanding of fundamental mechanisms leading to osteoporosis.
However, the differentially expressed genes in our study are predicted without experimental evidence, such as reverse transcription polymerase chain reaction and Western blot. Therefore, further research is needed to uncover gene functions in osteoporosis pathogenesis.
Osteoporosis is one of the most common and serious metabolic bone diseases, which affects 200 million people worldwide, among which 80% are women aged 60 years or older.1 Osteoporosis manifests clinically by reduced bone mass, a decreased amount of normally mineralised bone and microarchitectural deterioration of bone tissues, accounting for an increase in susceptibility to fracture.2 Genetic predisposition, combined with advanced age, gender, immobilisation, and other risk factors, contribute to the development of osteoporosis.3 Genome-wide association studies have been performed to search for an association between bone mineral density (BMD) and osteoporosis in twin and family studies.4,5 A meta-analysis of genome-wide association studies identified 56 low BMD-associated loci, and 14 loci which were associated with risk of fracture.6 The association of vitamin D (1,25- dihydroxyvitamin D3) receptor, oestrogen receptor 1 (ER) and LDL receptor related protein 5 with osteoporosis, has been most widely studied as candidate genetic genes.
The imbalance between bone resorption (by osteoclasts) and bone formation (by osteoblasts) is a key pathophysiological mechanism of osteoporosis as a consequence of pathologically prolonging the life span of osteoclasts and early apoptosis of osteoblasts and osteocytes.7-9 Circulating monocytes play important roles in this process, as they differentiate into osteoclasts to be involved in bone resorption,10,11 and produce a wide variety of factors for osteoclastogenesis and bone resorption.12 Osteoclasts from circulating mononuclear cells of osteoporotic patients demonstrated increased bone resorption activity compared with normal controls.13 Primary osteoporosis is related either to the postmenopausal state or advanced age (senile osteoporosis) due to the cellular changes associated with oestrogen deﬁciency and ageing. For postmenopausal osteoporosis, oestrogen deficiency leads to upregulation of receptor activator of nuclear factor-kappa ligand (RANKL) and tumour necrosis factor-α, resulting in excessive bone resorption and rapid bone loss.14 For senile osteoporosis, there is also a signiﬁcant reduction in bone formation due to osteoblastic dysfunction. Osteoblasts and adipocytes both differentiate from mesenchymal stem cells (MSCs), and there is a shift from osteoblastogenesis to adipogenesis in the bone marrow with increasing age, causing reduction in the life span of osteoblasts.15,16
Recently, a number of genes and pathophysiological mechanisms have been identified for osteoporosis through microarray studies in peripheral blood monocytes (PBMs) of osteoporotic patients.17,18 However, there are inconsistencies among these studies due to the use of different sample sources, platforms and analysis techniques.19 Towards this end, we integrated gene expression data sets from multiple microarray studies of osteoporosis to overcome these limitations of individual studies. The identification of differentially expressed genes (DEGs) between osteoporosis and normal controls and their biological function would expand our understanding of fundamental mechanisms leading to osteoporosis, and may be meaningful for further diagnosis and therapy for this disease.
Materials and Methods
Eligible gene expression profiles of PBMs in osteoporosis
The Gene Expression Omnibus (GEO) database was searched to obtain gene expression profiling studies of osteoporosis subjects in PBMs.20 The following key search terms were used: “osteoporosis”, “peripheral blood monocytes”, “gene expression”, “microarray”. We included original microarray studies which analysed the differential gene expression profiles of PBMs between patients with osteoporosis and normal controls (Fig. 1).
The existence of heterogeneity among multiple microarray studies arising from different microarray platforms, gene nomenclature and clinical samples, makes it infeasible to compare the gene expression data directly. Therefore, normalisation is necessary to minimise the heterogeneity. Consequently, we preprocessed the raw microarray data of each study by quantile normalisation and log2 transformation to obtain the z-score for the purpose of global normalisation.
The MATrix LABoratory (MATLAB) software (MathWorks, Natick, Massachusetts) was used to identify the differently expressed probe sets in PBMs between osteoporosis patients and normal controls. Gene specific t-tests were carried out, and p-values were calculated. Multiple testing adjustment was performed. The genes with a p-value < 0.05 were selected as DEGs.21 Heat map analysis was performed using the “heatmap.2” function of the R/Bioconductor package “gplots” (R Foundation for Statistical Computing, Vienna, Austria).22
Functional annotation of DEGs
Gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed to detect the biological functions and potential pathway of DEGs. GO-rilla,23 a web-accessible programme was used for GO enrichment analysis. The online-based software GENECODIS was used24 in pathway enrichment analysis based on the KEGG pathway database which is a recognised and comprehensive database including an extensive array of biochemistry pathways.25
Protein-protein interaction (PPI) network construction
The PPIs analysis is expected to explore the functions of proteins at the molecular level. PPIs may provide novel insights into the regulatory mechanisms of cellular growth, development, metabolism, differentiation and apoptosis.26 The Biological General Repository for Interaction Datasets (BioGRID),27 which is an online interaction repository, was used to construct the PPI network. The PPI network was visualised using Cytoscape software (Institute of Systems Biology, Seattle, Washington).28 Nodes denote proteins and edges denote interactions between two proteins.
Gene expression changes in PBMs between osteoporosis and NC
In this study, a total of three microarray studies (GSE7429, GSE2208 and GSE7158) were available by searching the GEO database, all of which analysed the differential gene expression profiles of PBMs between osteoporosis patients and normal controls. The included blood samples were from unrelated pre- and postmenopausal women. Patients with a hip T-score ⩽ -2.5 were deﬁned as having osteoporosis according to the World Health Organization criteria.29 The details of the included microarray studies are presented in Table I.17,30,31 According to the selection criteria (p-value < 0.05), 1125 genes were found to be signiﬁcantly differentially expressed in PBMs between osteoporosis patients and normal controls, with 373 upregulated and 752 downregulated genes. The top ten most signiﬁcantly up- or downregulated genes are listed in Table II. The upregulated gene with the lowest p-value was intraflagellar transport 52 (IFT52), which is a core component of the IFT particle. IFT is responsible for the formation and maintenance of flagella and cilia, and is involved in bone development by regulating cell differentiation in osteoblasts, osteocytes, chondrocytes, and MSCs.32 The downregulated gene with the lowest p-value was glucose-fructose oxidoreductase domain containing 2 (GFOD2), which has been found to be associated with blood lipid levels.33 The full list of the differentially genes is provided as supplementary material. The top 100 most significant DEGs are displayed in a heat map across different microarray studies of osteoporosis (Fig. 2).
For the identified DEGs, we performed enrichment analysis in GO categories including biological process and molecular function. Positive regulation of the cellular amine metabolic process (GO: 0033240, false discovery rate (FDR) = 1.00E + 00) and neural tube formation (GO: 0050681, FDR = 1.00E + 00) were significantly enriched under the GO category of biological processes (Fig. 3a), while for molecular functions, flavin adenine dinucleotide binding (GO: 0050660, FDR = 3.66E-01) and androgen receptor binding (GO: 0050681, FDR = 6.35E-01) were significantly enriched (Fig. 3b).
Furthermore, KEGG pathway enrichment analysis was also conducted. Hypergeometric tests with p-value < 0.05 were identified as the criteria for significantly enriched pathways. The most signiﬁcantly enriched pathway was measles (FDR = 6.53E-06). Furthermore, neuroactive ligand-receptor interaction (FDR = 5.08E-05) and the mitogen-activated protein kinase (MAPK) signalling pathway (FDR = 5.78E-05) are also highly enriched. Interestingly, the calcium signalling pathway (FDR = 5.83E-05) was found to be signiﬁcantly enriched.
PPI network construction
PPI networks of the top ten upregulated and downregulated DEGs were established by Cytoscape software, consisting of 398 nodes and 433 edges. In the PPI network, the nodes with ‘high degree’ are defined as ‘hub proteins’, and degrees are defined by the number of neighbours directly connected to the node. The significant hub proteins contained USP9X (Degree = 99), USP19 (Degree = 57) and UBE2B (Degree = 57) (Fig. 4).
Recent studies have advanced our understanding of pathophysiological mechanisms for osteoporosis. Several key pathways such as the RANK-RANKL-osteoprotegerin pathway, MSC differentiation, endochondral ossification and the Wnt signalling pathways, have been found to be involved in the regulation of bone mineral homeostasis.34 Based on the published evidence, novel therapeutic targets for osteoporosis have been identified and applied in a clinical setting, such as anti-RANKL receptor activator of NF-κB ligand.35 In this study, we integrated three microarray studies, of which gene expression profiles of PBMs between osteoporosis patients and normal controls were compared with the aim of identifying new candidates to target for osteoporosis drug discovery.
We identified 1125 DEGs in PBMs between osteoporosis patients and normal controls, which included 373 upregulated and 752 downregulated genes. Direct connections have been found between some of these genes and osteoporosis, such as nuclear receptor interacting protein 1(NRIP1), interleukin 1 receptor associated kinase 3 (IRAK3), bone morphogenetic protein 7 (BMP7), SMAD family member 1 (SMAD1), and MAPK3. The single-nucleotide polymorphisms (SNPs) of the NRIP1 gene which modulates transcriptional activity of the oestrogen receptor, were correlative with postmenopausal osteoporosis.36,37 Chen and Xia38 showed that there was a close relationship between NRIP1 and BMD in B cells in a microarray study. The expression of IRAK3 is predominantly in peripheral blood leukocytes and induced strongly upon the maturation of macrophages.39 Li et al40 reported that the IRAK3 plays a critical role as a negative regulator of toll-like receptor/interleukin (IL)-1 receptor (TLR/IL-1R) signalling in osteoclast differentiation and activation and, concomitantly, in the development of osteoporosis.
In this study, we also found several other genes which have been implicated in bone and skeletal muscle development, such as IFT52, GATA binding protein 6 (GATA6), tripartite motif containing 25 (TRIM25) and ubiquitin specific peptidase 19 (USP19), but the correlation with osteoporosis has never been reported. IFT52 is a core component of the IFT particle, which is involved in bone development by regulating cell differentiation in osteoblasts, osteocytes, chondrocytes, and MSCs.33 GATA6, as a direct target of micro-RNA miR-181a, plays a key role during BMP-induced osteoblastic differentiation of C2C12 (a mouse myoblast cell line) and MC3T3 cells (a clonal osteoblast-like mouse calvarial cell line).41 TRIM25, as one of many oestrogen-responsive genes, functions in the general maintenance of healthy bone.42
Furthermore, the results from PPI network analysis of the top ten upregulated and downregulated DEGs indicated that the most significant hub protein was USP9X. USP9X, located on the X-chromosome, plays a pivotal role in development and cancer by regulating transforming growth factor (TGF) signalling, Notch, epidermal growth factor, Wnt and mechanistic target of rapamycin signalling pathways. USP9X has been shown to regulate the apoptotic signalling networks.43 Osteoporosis is associated with dysregulated apoptosis based on the evidence that osteoporosis occurs mostly as a result of oestrogen deficiency which accelerates osteoblast apoptosis and susceptibility to osteoporotic fractures.44
KEGG pathway enrichment analysis showed that DEGs were enriched in many osteoporosis-related signalling pathways including the MAPK signalling pathway and the calcium signalling pathway. MAPK signalling pathways contribute greatly to osteoblast differentiation and bone formation via TGF-β and BMP signalling pathways.45,46 Vitamin D deficiency and reduced calcium absorption in the elderly population may be a risk factor of osteoporosis. An increase in calcium intake is the most commonly recommended preventive measure for osteoporosis.47
However, the DEGs in our study are predicted without experimental evidence such as reverse transcription polymerase chain reaction and Western blot. Therefore, further research is needed to uncover gene functions in osteoporosis pathogenesis. Overall, in our study a number of DEGs were identiﬁed and some of them may play important roles in the pathogenesis of osteoporosis. These ﬁndings may be useful in identifying the molecular mechanisms of osteoporosis and advancing therapy development.
The present study has some limitations. First, although global normalisation was performed, there may exist heterogeneity among the individual microarray studies arising from different platforms and clinical samples (such asseverity, gender, or geographical regions). Secondly, results from different software algorithms may be slightly different. In the process of data analysis, we selected software algorithms (e.g. GO-rilla, GENECODIS, BioGRID) that are commonly used based on relevant references.23,24,27
In conclusion, we used an integrated analysis approach to integrate three gene expression data sets of osteoporosis, and identified DEGs in PBMs between osteoporosis patients and normal controls. The gene function analysis showed that these DEGs were highly correlated with the development of osteoporosis, expanding our understanding of fundamental mechanisms leading to osteoporosis, and may be meaningful for further diagnosis and therapy for this disease.
Supplementary material Table showing the top 15 most significantly enriched KEGG pathway of DEGs is available alongside the online version of this article at www.bjr.boneandjoint.org.uk
Funding Statement This study was supported by a grant from Beijing High-level health technical personnel training plan (No: 2015-3-009).
ICMJE conflict of interest None declared
- Received August 3, 2016.
- Accepted October 5, 2016.
- © 2016 Fei et al.
This is an open-access article distributed under the terms of the Creative Commons Attributions licence (CC-BY-NC), which permits unrestricted use, distribution, and reproduction in any medium, but not for commercial gain, provided the original author and source are credited.