Skip to main content

ORIGINAL RESEARCH article

Front. Cell Dev. Biol., 27 August 2020
Sec. Molecular and Cellular Pathology
This article is part of the Research Topic Omics Data Integration towards Mining of Phenotype Specific Biomarkers in Cancers and Diseases View all 67 articles

Identification and Analysis of Genes Underlying Bone Mineral Density by Integrating Microarray Data of Osteoporosis

\r\nHaihong ZhangHaihong Zhang1Jinghui FengJinghui Feng2Zhiguo LinZhiguo Lin1Shuya WangShuya Wang1Yan WangYan Wang1Siming DaiSiming Dai1Weisi KongWeisi Kong1Yanli WangYanli Wang1Zhiyi Zhang*Zhiyi Zhang1*
  • 1Department of Rheumatology and Immunology, The First Affiliated Hospital of Harbin Medical University, Harbin, China
  • 2Department of Gerontology, The First Affiliated Hospital of Harbin Medical University, Harbin, China

Osteoporosis is a kind of brittle bone disease, which is characterized by a reduction in bone mineral density (BMD). In recent years, a number of genes and pathophysiological mechanisms have been identified for osteoporosis. However, the genes associated with BMD remain to be explored. Toward this end, we integrated multiple osteoporosis microarray datasets to identify and systematically characterize BMD-related genes. By integrating the differentially expressed genes from three osteoporosis microarray datasets, 152 genes show differentially expressed between high and low BMD osteoporosis samples in at least two of the three datasets. Among them, 88 were up-regulated in high BMD samples and 64 were up-regulated in low BMD samples. The expression of ZFP36, JUNB and TMEM8A were increased at high BMD samples in all three datasets. Hub genes were further identified by co-expression network analysis. Functional enrichment analysis showed that the gene up-regulated in high BMD were enriched in immune-related functions, suggesting that the immune system plays an important role in osteoporosis. Our study explored BMD-related genes based on the integration of osteoporosis microarray data, providing guidance to other researchers from a new perspective.

Introduction

Osteoporosis, a common systemic bone disease, can lead to weak bones and increase the risk of fractures (Ensrud and Crandall, 2017; Wilson et al., 2017; Khosla et al., 2018). A third of women and a fifth of men over the age of 50 have broken bones due to osteoporosis (Brown, 2017). Osteoporosis is characterized by low bone mass, microstructure degeneration and reduced bone strength. Patients often have a reduction in bone mineral density (BMD) (Coughlan and Dockery, 2014; Black and Rosen, 2016). Bone homeostasis depends on osteoclast absorption and osteoblast formation. The imbalance of this tightly coupled process can lead to the development of osteoporosis (Chen et al., 2018b). This process involves changes in a variety of signaling pathways, including MAPK signaling pathway, NF-κB pathway, Notch signaling pathway, etc. (Chen et al., 2018a, 2019; Lee and Long, 2018; Wang et al., 2019). Notch signaling pathway regulates the differentiation and function of osteoblasts and osteoclasts and participates in the process of bone reconstruction, activation of notch signaling pathway can inhibit glucose metabolism and osteoblast differentiation of bone marrow mesenchymal progenitor cells (Lee and Long, 2018).

In recent years, some genes and pathophysiological mechanisms have been identified by microarray analysis for patients with osteoporosis (Liu et al., 2005; Lei et al., 2009; Li et al., 2016). The study based on peripheral blood monocyte cells (PBMCs) microarray data of osteoporosis patients revealed the pathophysiological mechanism of osteoporosis, which is characterized by increased recruitment of monocytes into bone and then differentiating into osteoclasts (Liu et al., 2005). Zhou et al. (2018b, 2019) predicted osteoporosis related transcription factors (TF) and long non-coding RNA (lncRNA) via exon arrays. Another study explored osteoporosis-related pathways based on microarray data (Zhou et al., 2018a). The studies mentioned above were based on a limited set of data, and have certain limitations. In order to obtain more robust results, it is critical to integrate multiple datasets for obtaining new insights.

In this study, we collected three microarray datasets of osteoporosis. Differential expression analysis was conducted to identify differentially expressed genes (DEGs) between high and low BMD samples. Then DEGs were integrated to obtain uniformly expressed BMD-related genes. Analysis of gene co-expression networks revealed key genes of different BMD conditions. And enrichment analysis showed that they were enriched in immune-related functions and biological pathways, suggesting a potential role of the immune system in osteoporosis.

Materials and Methods

GEO Datasets

Microarray datasets of high and low BMD osteoporosis samples were downloaded from Gene Expression Omnibus (GEO) database1, including GSE2208, GSE56814 and GSE56815. A total of 172 samples were collected from the three datasets, of which 92 were high BMD samples and 80 were low BMD samples (Table 1).

TABLE 1
www.frontiersin.org

Table 1. The overview of three GEO datasets.

Identification of Differentially Expressed Genes (DEGs)

The up-regulated genes in high or low BMD samples of each dataset were identified using the R package named limma with a threshold of |log2FoldChange|>0 and P < 0.05. Genes up-regulated in at least two datasets were identified as uniformly expressed BMD-related genes and used for subsequent analysis.

The Construction of Gene Co-expression Network

We calculated the Pearson Correlation Coefficient (PCC) of integrated up-regulated genes in high or low BMD samples (PCC > 0.4 and P < 0.05). For co-expressed gene pairs that appeared in multiple datasets, the averaged PCC was calculated. Then Cytoscape (V.3.8.0) was used for network visualization, and hub genes were identified by cytoHubba plugin (Shannon et al., 2003; Chin et al., 2014).

Enrichment Analysis

We performed GO and KEGG enrichment analysis of genes included in co-expression network using R package “clusterProfiler” (Yu et al., 2012). The significance threshold is 0.05.

Results

Differentially Expressed Genes Between High and Low BMD Osteoporosis Samples

To explore the key genes associated with BMD, we downloaded three osteoporosis microarray datasets from NCBI GEO database. In total, 19718 genes and 172 samples were involved. The number of high and low BMD samples in each dataset was comparable (Table 1).

We performed the analysis of differentially expressed genes between the two types of samples (Figure 1). There were 12.59% (GSE2208), 3.30% (GSE56814), and 18.49% (GSE56815) genes identified as differentially expressed, respectively. In GSE56814, compared to the low BMD samples, there were more genes up-regulated in high BMD samples, while GSE56814 had more genes up-regulated in low BMD samples. In GSE2208, the two types of genes were almost the same in number (Figure 1D).

FIGURE 1
www.frontiersin.org

Figure 1. Differential expression analysis. (A–C) Volcano plot of DEGs in each dataset, red nodes represent upregulation in high BMD samples and blue nodes represent upregulation in low BMD samples. (D) Statistics of two types of DEGs.

Integration of BMD-Related Genes

We integrated the DEGs of three datasets to obtain robust different BMD-related genes. Based on the criterion that the same gene was up-regulated in at least two datasets, we identified 88 and 64 uniformly expressed genes with high and low BMD by UpSetR, respectively (Figures 2A,B and Supplementary Tables 1, 2; Conway et al., 2017). In all datasets, there were three genes, ZFP36, JUNB and TMEM8A shown up-regulated in high BMD (Figure 2C and Table 2). JUNB is a member of c-Jun protein family, which interacts with c-Fos protein family to form transcription factor AP-1. AP-1 can activate osteoclast specific genes (Wagner and Eferl, 2005; Asagiri and Takayanagi, 2007; Hamamura et al., 2015). Only one gene, CACNA2D3, had elevated expression with low BMD in all datasets (Figure 2C and Table 2).

FIGURE 2
www.frontiersin.org

Figure 2. Integration of BMD-related DEGs. (A) UpSetR plot of genes up-regulated in high BMD samples, the number of genes present in at least two datasets is marked (red). (B) UpSetR plot of genes up-regulated in low BMD samples, the number of genes present in at least two datasets is marked (red). (C) Boxplot of genes that expressed consistently across all three datasets.

TABLE 2
www.frontiersin.org

Table 2. Genes expressed consistently across all three sets of data.

Co-expression Network and Functional Analysis of BMD-Related Genes

Through the similarity of gene expression, the possible interaction between gene products can be analyzed, so as to understand the interaction between genes and find out the core genes. Therefore, we performed a co-expression analysis of integrated up-regulated genes in high and low BMD samples, respectively. We obtained 58 co-expressed gene pairs of high BMD. The top5 hub genes were identified by cytoHubba, a plug-in of Cytoscape, including KDM2A, APH1A, DNPEP, NFKBIB, and TMEM8A. KDM2A was co-expressed with the largest number of genes (Figure 3A). KDM2A can regulate Mesenchymal stem cells (MSCs) osteo/dentinogenic differentiation and cell proliferation (Dong et al., 2013; Gao et al., 2013). Co-expression network of genes that were up-regulated in high BMD samples contained 48 gene nodes. We did functional enrichment analysis on them. In addition to osteoporosis-related functions and pathways, Notch signaling pathway and osteoclast differentiation, enrichment analysis results showed that they were mainly enriched in immune-related functions and pathways (Figure 3B; Boyle et al., 2003; Lee and Long, 2018). Studies have shown that the immune system plays an important role in osteoporosis (Faienza et al., 2013). For example, postmenopausal osteoporosis patients had higher T-cell activity and increased TNFα and RANKL production, which can promote osteoclast differentiation (D’Amelio et al., 2008; Mirza et al., 2010; Kim et al., 2012). In addition, other studies have shown that B lymphocytes were closely related to bone metabolism (Takayanagi et al., 2005; Boyce and Xing, 2008; Breuil et al., 2010).

FIGURE 3
www.frontiersin.org

Figure 3. Co-expression network and function enrichment analysis of genes that are consistently expressed in high BMD samples. (A) The gene co-expression network. Node size represents how many other genes interact with it. Node color indicates in which data sets it is up-regulated in high BMD samples. The darker the edge color, the greater the correlation coefficient. (B) GO and KEGG results of genes in network.

In the co-expression network of genes that were up-regulated in low BMD samples, there were 44 gene pairs, involving 38 genes (Figure 4A). The top5 hub genes were CCT7, DGUOK, MPHOSPH10, RARS, and DMTF1. The results of enrichment analysis showed that low BMD-related genes were mainly enriched in basic biological processes and pathways, including ribosome biogenesis, rRNA processing and translation elongation (Figure 4B).

FIGURE 4
www.frontiersin.org

Figure 4. Co-expression network and function enrichment analysis of genes that are consistently expressed in low BMD samples. (A) The gene co-expression network. Node size represents how many other genes interact with it. Node color indicates in which datasets it is up-regulated in low BMD samples. The darker the edge color, the greater the correlation coefficient. (B) GO and KEGG results of genes in network.

Discussion

In this study, we identified and analyzed BMD-related genes through three osteoporosis microarray datasets. Differential expression analysis identified DEGs between high and low BMD in each dataset. By integration, we screened out 152 uniformly differentially expressed genes. Gene co-expression network analysis further identified key top 5 hub genes. Co-expression genes that had elevated expression in high BMD were enriched in functions of immune systems, suggesting its important potential role in osteoporosis. In addition, they were also enriched in osteoclast differentiation pathway, indicating that high BMD is developing toward low BMD.

In conclusion, we used bioinformatic methods to systematically characterize genes underlying BMD levels based on osteoporosis microarray data. We analyzed osteoporosis from a new perspective and provided new guidance for its diagnosis and treatment.

Data Availability Statement

All datasets presented in this study are included in the article/Supplementary Material.

Author Contributions

ZZ designed the study. HZ, JF, ZL, SW, and YW analyzed the data. SD, WK, and YLW wrote the manuscript. All authors read and approved the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2020.00798/full#supplementary-material

TABLE S1 | The uniformly up-regulated DEGs in high BMD samples.

TABLE S2 | The uniformly up-regulated DEGs in low BMD samples.

Footnotes

  1. ^ https://www.ncbi.nlm.nih.gov/geo/

References

Asagiri, M., and Takayanagi, H. (2007). The molecular understanding of osteoclast differentiation. Bone 40, 251–264. doi: 10.1016/j.bone.2006.09.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Black, D. M., and Rosen, C. J. (2016). Clinical practice. postmenopausal osteoporosis. N. Engl. J. Med. 374, 254–262.

Google Scholar

Boyce, B. F., and Xing, L. (2008). Bruton and tec: new links in osteoimmunology. Cell Metab 7, 283–285. doi: 10.1016/j.cmet.2008.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyle, W. J., Simonet, W. S., and Lacey, D. L. (2003). Osteoclast differentiation and activation. Nature 423, 337–342. doi: 10.1038/nature01658

PubMed Abstract | CrossRef Full Text | Google Scholar

Breuil, V., Ticchioni, M., Testa, J., Roux, C. H., Ferrari, P., Breittmayer, J. P., et al. (2010). Immune changes in post-menopausal osteoporosis: the Immunos study. Osteoporos. Int. 21, 805–814. doi: 10.1007/s00198-009-1018-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, C. (2017). Osteoporosis: staying strong. Nature 550, S15–S17.

Google Scholar

Chen, E., Liu, G., Zhou, X., Zhang, W., Wang, C., Hu, D., et al. (2018a). Concentration-dependent, dual roles of IL-10 in the osteogenesis of human BMSCs via P38/MAPK and NF-kappaB signaling pathways. FASEB J. 32, 4917–4929. doi: 10.1096/fj.201701256rrr

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Wang, Z., Duan, N., Zhu, G., Schwarz, E. M., and Xie, C. (2018b). Osteoblast-osteoclast interactions. Connect. Tissue Res. 59, 99–107. doi: 10.1080/03008207.2017.1290085

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, K., Qiu, P., Yuan, Y., Zheng, L., He, J., Wang, C., et al. (2019). Pseurotin a inhibits osteoclastogenesis and prevents ovariectomized-induced bone loss by suppressing reactive oxygen species. Theranostics 9, 1634–1650. doi: 10.7150/thno.30206

PubMed Abstract | CrossRef Full Text | Google Scholar

Chin, C. H., Chen, S. H., Wu, H. H., Ho, C. W., Ko, M. T., and Lin, C. Y. (2014). cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst. Biol. 8(Suppl. 4):S11. doi: 10.1186/1752-0509-8-S4-S11

PubMed Abstract | CrossRef Full Text | Google Scholar

Conway, J. R., Lex, A., and Gehlenborg, N. (2017). UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics 33, 2938–2940. doi: 10.1093/bioinformatics/btx364

PubMed Abstract | CrossRef Full Text | Google Scholar

Coughlan, T., and Dockery, F. (2014). Osteoporosis and fracture risk in older people. Clin. Med. (Lond) 14, 187–191. doi: 10.7861/clinmedicine.14-2-187

PubMed Abstract | CrossRef Full Text | Google Scholar

D’Amelio, P., Grimaldi, A., Di Bella, S., Brianza, S. Z. M., Cristofaro, M. A., Tamone, C., et al. (2008). Estrogen deficiency increases osteoclastogenesis up-regulating T cells activity: a key mechanism in osteoporosis. Bone 43, 92–100. doi: 10.1016/j.bone.2008.02.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, R., Yao, R., Du, J., Wang, S., and Fan, Z. (2013). Depletion of histone demethylase KDM2A enhanced the adipogenic and chondrogenic differentiation potentials of stem cells from apical papilla. Exp. Cell Res. 319, 2874–2882. doi: 10.1016/j.yexcr.2013.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Ensrud, K. E., and Crandall, C. J. (2017). Osteoporosis. Ann. Intern. Med. 167, ITC17–ITC32.

Google Scholar

Faienza, M. F., Ventura, A., Marzano, F., and Cavallo, L. (2013). Postmenopausal osteoporosis: the role of immune system cells. Clin. Dev. Immunol. 2013:575936.

Google Scholar

Gao, R., Dong, R., Du, J., Ma, P., Wang, S., and Fan, Z. (2013). Depletion of histone demethylase KDM2A inhibited cell proliferation of stem cells from apical papilla by de-repression of p15INK4B and p27Kip1. Mol. Cell. Biochem. 379, 115–122. doi: 10.1007/s11010-013-1633-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamamura, K., Chen, A., Tanjung, N., Takigawa, S., Sudo, A., and Yokota, H. (2015). In vitro and in silico analysis of an inhibitory mechanism of osteoclastogenesis by salubrinal and guanabenz. Cell Signal 27, 353–362. doi: 10.1016/j.cellsig.2014.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Khosla, S., Farr, J. N., and Kirkland, J. L. (2018). Inhibiting cellular senescence: a new therapeutic paradigm for age-related osteoporosis. J. Clin. Endocrinol. Metab. 103, 1282–1290. doi: 10.1210/jc.2017-02694

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, B. J., Bae, S. J., Lee, S. Y., Lee, Y. S., Baek, J. E., Park, S. Y., et al. (2012). TNF-alpha mediates the stimulation of sclerostin expression in an estrogen-deficient condition. Biochem. Biophys. Res. Commun. 424, 170–175. doi: 10.1016/j.bbrc.2012.06.100

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S. Y., and Long, F. (2018). Notch signaling suppresses glucose metabolism in mesenchymal progenitors to restrict osteoblast differentiation. J. Clin. Invest. 128, 5573–5586. doi: 10.1172/jci96221

PubMed Abstract | CrossRef Full Text | Google Scholar

Lei, S. F., Wu, S., Li, L. M., Deng, F. Y., Xiao, S. M., Jiang, C., et al. (2009). An in vivo genome wide gene expression study of circulating monocytes suggested GBP1, STAT1 and CXCL10 as novel risk genes for the differentiation of peak bone mass. Bone 44, 1010–1014. doi: 10.1016/j.bone.2008.05.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J. J., Wang, B. Q., Fei, Q., Yang, Y., and Li, D. (2016). Identification of candidate genes in osteoporosis by integrated microarray analysis. Bone Joint Res. 5, 594–601. doi: 10.1302/2046-3758.512.bjr-2016-0073.r1

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y. Z., Dvornyk, V., Lu, Y., Shen, H., Lappe, J. M., Recker, R. R., et al. (2005). A novel pathophysiological mechanism for osteoporosis suggested by an in vivo gene expression study of circulating monocytes. J. Biol. Chem. 280, 29011–29016. doi: 10.1074/jbc.m501164200

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirza, F. S., Padhi, I. D., Raisz, L. G., and Lorenzo, J. A. (2010). Serum sclerostin levels negatively correlate with parathyroid hormone levels and free estrogen index in postmenopausal women. J. Clin. Endocrinol. Metab. 95, 1991–1997. doi: 10.1210/jc.2009-2283

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Takayanagi, H., Sato, K., Takaoka, A., and Taniguchi, T. (2005). Interplay between interferon and other cytokine systems in bone metabolism. Immunol. Rev. 208, 181–193. doi: 10.1111/j.0105-2896.2005.00337.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, E. F., and Eferl, R. (2005). Fos/AP-1 proteins in bone and the immune system. Immunol. Rev. 208, 126–140.

Google Scholar

Wang, F., Zhang, C., Ge, W., and Zhang, G. (2019). Up-regulated CST5 inhibits bone resorption and activation of osteoclasts in rat models of osteoporosis via suppression of the NF-kappaB pathway. J. Cell Mol. Med. 23, 6744–6754.

Google Scholar

Wilson, L. M., Rebholz, C. M., Jirru, E., Liu, M. C., Zhang, A., Gayleard, J., et al. (2017). Benefits and harms of osteoporosis medications in patients with chronic kidney disease: a systematic review and meta-analysis. Ann. Intern. Med. 166, 649–658.

Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287.

Google Scholar

Zhou, Y., Gao, Y., Xu, C., Shen, H., Tian, Q., and Deng, H. W. (2018a). A novel approach for correction of crosstalk effects in pathway analysis and its application in osteoporosis research. Sci. Rep. 8:668.

Google Scholar

Zhou, Y., Zhu, W., Zhang, L., Zeng, Y., Xu, C., Tian, Q., et al. (2018b). Transcriptomic data identified key transcription factors for osteoporosis in caucasian women. Calcif. Tissue Int. 103, 581–588.

Google Scholar

Zhou, Y., Xu, C., Zhu, W., He, H., Zhang, L., Tang, B., et al. (2019). Long noncoding RNA analyses for osteoporosis risk in caucasian women. Calcif. Tissue Int. 105, 183–192.

Google Scholar

Keywords: bone mineral density, osteoporosis, microarray, co-expression, enrichment analysis

Citation: Zhang H, Feng J, Lin Z, Wang S, Wang Y, Dai S, Kong W, Wang Y and Zhang Z (2020) Identification and Analysis of Genes Underlying Bone Mineral Density by Integrating Microarray Data of Osteoporosis. Front. Cell Dev. Biol. 8:798. doi: 10.3389/fcell.2020.00798

Received: 02 July 2020; Accepted: 28 July 2020;
Published: 27 August 2020.

Edited by:

Lei Deng, Central South University, China

Reviewed by:

Hao Lin, University of Electronic Science and Technology of China, China
Liang Yu, Xidian University, China

Copyright © 2020 Zhang, Feng, Lin, Wang, Wang, Dai, Kong, Wang and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Zhiyi Zhang, zhangzhiyi2014@163.com

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.