# The Hub Genes Related to Osteoporosis Were Identified by Bioinformatics Analysis

Authors: Mengdie Hu; Hong Ding; Rui Chao; Zhidong Cao
Journal: BioMed Research International
 (2023)

---

## Abstract

Osteoporosis (OP) is commonly encountered, which is a kind of systemic injury of bone mass and microstructure, leading to brittle fractures. With the aging of the population, this disease will pose a more serious impact on medical, social, and economic aspects, especially postmenopausal osteoporosis (PMOP). This study is aimed at figuring out potential therapeutic targets and new biomarkers in OP via bioinformatics tools. After differentially expressed gene (DEG) analysis, we successfully identified 97 upregulated and 172 downregulated DEGs. They are mainly concentrated in actin binding, regulation of cytokine production, muscle cell promotion, chemokine signaling pathway, and cytokine-cytokine receiver interaction. According to the diagram of protein-protein interaction (PPI), we obtained 10 hub genes: CCL5, CXCL10, EGFR, HMOX1, IL12B, CCL7, TBX21, XCL1, PGR, and ITGA1. Expression analysis showed that Egfr, Hmox1, and Pgr were significantly upregulated in estrogen-treated osteoporotic patients, while Ccl5, Cxcl10, Il12b, Ccl7, Tbx21, Xcl1, and Itga1 were significantly downregulated. In addition, the analysis results of Pearson’s correlation revealed that CCL7 has a strong positive association with IL12b, TBX21, and CCL5 and so was CCL5 with IL12b. Conversely, EGFR has a strong negative association with XCL1 and CXCL10. In conclusion, this study screened 10 hub genes related to OP based on the GEO database, laying a biological foundation for further research on new biomarkers and potential therapeutic targets in OP.


---

## Body

## 1. Introduction

Osteoporosis (OP) is a common systemic disease infiltrating bones, and it is age-related with bone microstructure destruction, bone mass reduction, and worsening bone fragility, as a result of the substantial increase in fracture occurrence risk [1]. OP were also accompanied with a constellation of complications that can greatly influence the quality of life of the victims, namely, fracture, chronic pain, or disability [2, 3]. Among all OP complications, a fracture affects the most with an estimated over 8.9 million fractures associated with OP every year, resulting in considerable health, social, and economic burden [4]. With the aggravation of population aging, OP has emerged as a heavy medical burden around the world [5].The past decade sees a big increase in OP prevalence in China, and more than 1/3 of people over 50 years have been affected [6]. OP incidence is higher in females than in males, and PMOP accounts for the majority of female patients [7]. PMOP is a kind of primary OP with a high incidence rate. Postmenopausal women are accompanied by a sharp decline in hormone levels in addition to age-related factors [8]. Hence, the improvement of diagnosis and treatment of postmenopausal OP is of great significance in clinical practice.OP is a highly insidious disease. Due to the lack of overt symptoms and sensitive biomarkers, many patients are only diagnosed after a fracture has occurred. Studies have now shown the presence of biomarkers in the OP, such as some potential biomarkers of OP that have been found in metabolomics studies. These metabolites are mainly concentrated in fats and amino acids [9]. Studies have shown that circulating miR-181c-5p and miR-497-5p may serve as potential biomarkers to monitor the effect of anti-OP therapy or diagnostic methods [10]. By integrating feature selection of SVM-RFE and RF, six genes (HNRNPC, PFDN2, PSMC5, RPS16, TCEB2, and UBE2V2) were selected as genetic biomarkers for smoking-related postmenopausal OP (SRPO) [7]. However, clinically effective markers are rarely reported; so, we should find more effective markers to improve the diagnosis of OP.In addition, the first-line drugs for the treatment of OP link to multiple complications produce an unsatisfactory effect on the overall management of this disease. With the development of several pharmaceutical technologies, novel medications, calcitonin, bisphosphonates, denosumab, and teriparatide, have been widely introduced for the prevention and treatment of OP by targeting osteoblast bone formation and resorption activity [11, 12]. Unfortunately, such kinds of therapies can cause severe adverse events, hindering the application of the drugs and compliance of patients [13]. The etiology of OP and its therapeutic targets required has not been fully understood yet. It is therefore that the development of alternative methods for OP management is in urgent need to improve the efficacy and reduce the toxicity.


## 2. Method

### 2.1. Data Download

Gene expression profile data GSE163849 was screened and downloaded from the GEO database using the following link (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE163849). It is used to screen DEGs in OP. Eight specimens (four experimental group samples and four for control) were included.


### 2.2. Data Analysis

Gene differential expression analysis usesR package DESeq2 to analyze the differential genes among the samples. After P values were calculated, multiple hypothesis test correction was carried out. A false discovery rate (FDR) was employed to determine the threshold of the P values. Q value represented the corrected P values. Subsequently, we calculated differential expression multiple according to FPKM values, that is, fold change. The screening index of this analysis is P value <0.05, log2FC>1, or < -1.


### 2.3. DEG Function Enrichment Analysis

DEGs with statistical significance was analyzed subsequently via the R3.6.1 software with “biological conductor” and “goplot” software packages, Gene Ontology (GO) functions, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways which were enriched as well [14, 15]. The adjusted P value <0.05 is considered the critical standard.


### 2.4. PPI Analysis

The search tool for retrieving interacting genes (STRING) (http://string.embl.de/) was used to construct the (protein-protein interaction) PPI network of the identified DEG and to present this network using Cytoscape visualization software version 3.6.1. The protein-protein interaction (PPI) network utilized to construct identified DEGs was presented using Cytoscape visualization software version 3.6.1. Then, modules with the most significance in the PPI diagram are determined by the molecular recombination detection tool (MCODE). The standard was set as score >2, cut-off degree = 2, and cut-off degree of node score = 0.2.


### 2.5. Cytoscape Hub Gene Analysis

Genes play an important role in biological processes. We constructed an interaction network through coexpression or protein interaction and screened out hub genes according to the network topology. CytoHubba, in Cytoscape, is primarily used to rank nodes in a network through its network function. CytoHubba offers 11 topological analysis methods, including degree, edge penetration component, maximum neighborhood component, maximum neighborhood component density, maximum group centrality, and six centroids (Bottleneck, EcCentricity, Affinity, Radiance, Intermediacy, and Stress). Of the 11 methods, MCC had better performance in terms of accuracy in predicting essential proteins from PPI networks.


## 2.1. Data Download

Gene expression profile data GSE163849 was screened and downloaded from the GEO database using the following link (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE163849). It is used to screen DEGs in OP. Eight specimens (four experimental group samples and four for control) were included.


## 2.2. Data Analysis

Gene differential expression analysis usesR package DESeq2 to analyze the differential genes among the samples. After P values were calculated, multiple hypothesis test correction was carried out. A false discovery rate (FDR) was employed to determine the threshold of the P values. Q value represented the corrected P values. Subsequently, we calculated differential expression multiple according to FPKM values, that is, fold change. The screening index of this analysis is P value <0.05, log2FC>1, or < -1.


## 2.3. DEG Function Enrichment Analysis

DEGs with statistical significance was analyzed subsequently via the R3.6.1 software with “biological conductor” and “goplot” software packages, Gene Ontology (GO) functions, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways which were enriched as well [14, 15]. The adjusted P value <0.05 is considered the critical standard.


## 2.4. PPI Analysis

The search tool for retrieving interacting genes (STRING) (http://string.embl.de/) was used to construct the (protein-protein interaction) PPI network of the identified DEG and to present this network using Cytoscape visualization software version 3.6.1. The protein-protein interaction (PPI) network utilized to construct identified DEGs was presented using Cytoscape visualization software version 3.6.1. Then, modules with the most significance in the PPI diagram are determined by the molecular recombination detection tool (MCODE). The standard was set as score >2, cut-off degree = 2, and cut-off degree of node score = 0.2.


## 2.5. Cytoscape Hub Gene Analysis

Genes play an important role in biological processes. We constructed an interaction network through coexpression or protein interaction and screened out hub genes according to the network topology. CytoHubba, in Cytoscape, is primarily used to rank nodes in a network through its network function. CytoHubba offers 11 topological analysis methods, including degree, edge penetration component, maximum neighborhood component, maximum neighborhood component density, maximum group centrality, and six centroids (Bottleneck, EcCentricity, Affinity, Radiance, Intermediacy, and Stress). Of the 11 methods, MCC had better performance in terms of accuracy in predicting essential proteins from PPI networks.


## 3. Results

### 3.1. Screening of DEGs

GSE163849 was selected, and the “limma” package inR software 3.6.1 was used for DEG analysis. There were 269 DEGs identified with 97 genes upregulated and 172 downregulated (log2fc>1, adjusted P values <0.05). Figure 1(a) presents all plotted 269 DEGs. Blue indicates downregulated genes, red indicates upregulated genes, and gray indicates other DEGs. In addition, the expression levels of all DEGs are shown in the heat map (Figure 1(b)), and these genes are well clustered between the estrogen treatment group and the control.Figure 1
Analysis of DEGs. (a) Volcanic map of DEGs. Blue indicates downregulation, red indicates upregulation, and gray indicates the rest of the DEGs. (b) DEG heat map. Columns represent data sets and rows represent genes. Green indicates genes downregulated and red indicates genes upregulated. EB: estrogen treatment group; CB: control group. DEGs: differentially expressed genes.
(a)(b)


### 3.2. Functional Enrichment Analysis of DEGs

Subsequently, we obtained the top 10 GO enriched terms consisting of 3 aspects: cellular component (CC), biological process (BP), and molecular function (MF) (Figures2 and 3). CC analysis showed that DEGs significantly enriched myofibril, sarcomere, I band, and Z disc. The first three significantly enriched are actin-binding, serine-type endopeptidase activity, and growth factor binding during MF analysis. In biological process (BP), DEGs significantly enriched regulation of smooth muscle cell promotion, response to a toxic substance, and negative regulation of cytokine production (Figure 2). The chord diagram shows that ligp1, ll12b, Xcl1, Cxcl10, Ccl7, Ccl5, Mmp12, and F3 are involved in cytokine-mediated signaling pathways. The genes involved in lymphocyte migration are Tbx21, Xcl1, Cxcl10, Ccl7, and Ccl5 (Figure 3). In addition, enriched KEGG pathways include cytokine-cytokine receiver interaction, Toll-like receiver signaling pathway, chemokine signaling pathway, and HIF-1 signaling pathway (Figure 4).Figure 2
GO analysis of DEGs (P value <0.01, Q value <0.05). DEGs: differentially expressed genes; GO: Gene Oncology.Figure 3
The first five enriched gene ontology terms.Figure 4
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (P values <0.05 and Q value <0.05).


### 3.3. PPI Network Diagram Analysis

There were 49 nodes (DEGs) in the PPI diagram. Ten key genes CCL5, CXCL10, EGFR, HMOX1, IL12B, CCL7, TBX21, XCL1, PGR, and ITGA1 were identified (Figure5(a)). Two important modules are presented: one contains three gene nodes MMP12, CXCL10, and CCL5, with 23 edges (interactions) (Figure 5(b)), and the other contains four upregulated genes, including IL12B, CCL7, TBX21, and XCL1, with 16 edges (Figure 5(c)).Figure 5
PPI diagram of differentially expressed genes and module analysis. (a) PPI diagram. (b) Gene nodes in module 1 with 23 edges. (c) Gene nodes in module 2 with 16 edges. PPI: protein-protein interaction.
(a)(b)(c)


### 3.4. Hub Gene Expression

The 10 hubs are presented using a plug-in “CytoHubba” of Cytoscape version 3.7.1 according to the degree (Figure6(a)). Compared with the control group, EGFR, HMOX1, and PGR were significantly upregulated in estrogen-treated osteoporotic patients, while CCL5, CXCL10, IL12B, CCL7, TBX21, XCL1, and ITGA1 were significantly downregulated (Figure 6(b)). The expression levels of 10 hub genes (CCL5, CXCL10, EGFR, HMOX1, IL12B, CCL7, TBX21, XCL1, PGR, and ITGA1) were analyzed using Pearson’s correlation analysis. A strong positive association was revealed between CCL7 and IL12B, TBX21, and CCL5. CCL5 had a positive association with IL12B. EGFR was significantly negatively correlated with XCL1 and CXCL10 (Figure 6(c)).Figure 6
Analysis of the hub gene expression. (a) Network of 10 hub genes. The values from high to low of the ten nodes are expressed from red to yellow. (b) Expression of 10 hub genes. (c) Map of correlation analysis among 10 hub genes. The correlation in red indicates positive, and blue indicates negative, respectively. EB: estrogen treatment group; CB: control group.
(a)(b)(c)


## 3.1. Screening of DEGs

GSE163849 was selected, and the “limma” package inR software 3.6.1 was used for DEG analysis. There were 269 DEGs identified with 97 genes upregulated and 172 downregulated (log2fc>1, adjusted P values <0.05). Figure 1(a) presents all plotted 269 DEGs. Blue indicates downregulated genes, red indicates upregulated genes, and gray indicates other DEGs. In addition, the expression levels of all DEGs are shown in the heat map (Figure 1(b)), and these genes are well clustered between the estrogen treatment group and the control.Figure 1
Analysis of DEGs. (a) Volcanic map of DEGs. Blue indicates downregulation, red indicates upregulation, and gray indicates the rest of the DEGs. (b) DEG heat map. Columns represent data sets and rows represent genes. Green indicates genes downregulated and red indicates genes upregulated. EB: estrogen treatment group; CB: control group. DEGs: differentially expressed genes.
(a)(b)


## 3.2. Functional Enrichment Analysis of DEGs

Subsequently, we obtained the top 10 GO enriched terms consisting of 3 aspects: cellular component (CC), biological process (BP), and molecular function (MF) (Figures2 and 3). CC analysis showed that DEGs significantly enriched myofibril, sarcomere, I band, and Z disc. The first three significantly enriched are actin-binding, serine-type endopeptidase activity, and growth factor binding during MF analysis. In biological process (BP), DEGs significantly enriched regulation of smooth muscle cell promotion, response to a toxic substance, and negative regulation of cytokine production (Figure 2). The chord diagram shows that ligp1, ll12b, Xcl1, Cxcl10, Ccl7, Ccl5, Mmp12, and F3 are involved in cytokine-mediated signaling pathways. The genes involved in lymphocyte migration are Tbx21, Xcl1, Cxcl10, Ccl7, and Ccl5 (Figure 3). In addition, enriched KEGG pathways include cytokine-cytokine receiver interaction, Toll-like receiver signaling pathway, chemokine signaling pathway, and HIF-1 signaling pathway (Figure 4).Figure 2
GO analysis of DEGs (P value <0.01, Q value <0.05). DEGs: differentially expressed genes; GO: Gene Oncology.Figure 3
The first five enriched gene ontology terms.Figure 4
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (P values <0.05 and Q value <0.05).


## 3.3. PPI Network Diagram Analysis

There were 49 nodes (DEGs) in the PPI diagram. Ten key genes CCL5, CXCL10, EGFR, HMOX1, IL12B, CCL7, TBX21, XCL1, PGR, and ITGA1 were identified (Figure5(a)). Two important modules are presented: one contains three gene nodes MMP12, CXCL10, and CCL5, with 23 edges (interactions) (Figure 5(b)), and the other contains four upregulated genes, including IL12B, CCL7, TBX21, and XCL1, with 16 edges (Figure 5(c)).Figure 5
PPI diagram of differentially expressed genes and module analysis. (a) PPI diagram. (b) Gene nodes in module 1 with 23 edges. (c) Gene nodes in module 2 with 16 edges. PPI: protein-protein interaction.
(a)(b)(c)


## 3.4. Hub Gene Expression

The 10 hubs are presented using a plug-in “CytoHubba” of Cytoscape version 3.7.1 according to the degree (Figure6(a)). Compared with the control group, EGFR, HMOX1, and PGR were significantly upregulated in estrogen-treated osteoporotic patients, while CCL5, CXCL10, IL12B, CCL7, TBX21, XCL1, and ITGA1 were significantly downregulated (Figure 6(b)). The expression levels of 10 hub genes (CCL5, CXCL10, EGFR, HMOX1, IL12B, CCL7, TBX21, XCL1, PGR, and ITGA1) were analyzed using Pearson’s correlation analysis. A strong positive association was revealed between CCL7 and IL12B, TBX21, and CCL5. CCL5 had a positive association with IL12B. EGFR was significantly negatively correlated with XCL1 and CXCL10 (Figure 6(c)).Figure 6
Analysis of the hub gene expression. (a) Network of 10 hub genes. The values from high to low of the ten nodes are expressed from red to yellow. (b) Expression of 10 hub genes. (c) Map of correlation analysis among 10 hub genes. The correlation in red indicates positive, and blue indicates negative, respectively. EB: estrogen treatment group; CB: control group.
(a)(b)(c)


## 4. Discussion

Bones with OP are fragile and tend to break easily. It is common in the elderly, especially in postmenopausal women [16]. Primary OP is frequently seen in females who are in menopause due to this disease linked to estrogen deficiency. Estrogen level decreases in this period, causing extensive resorption of bones instead of bone formation and consequently, OP occurs [17]. Animal studies in mice showed that under estrogen deficiency, callus formation and angiogenesis were impaired during fracture healing. Estrogen deficiency is associated with reduced fractal neoangiogenesis and delayed endochondral ossification [18].At present, bisphosphonates are introduced as the first-line drugs against OP, targeting OP management due to elevated levels of bone resorption [19]. However, many studies have reported severe complications caused by first-line medications including the application of chloromethyl bisphosphonates [20]. Hence, the improvement of diagnosis and treatment of postmenopausal OP is worthwhile. In this study, bioinformatics analysis was carried out based on the GEO database to identify the hub genes of OP associated with estrogen deficiency. Through DEG analysis, 97 significantly upregulated and 172 significantly downregulated genes were successfully identified. GO analysis shows that DEGs are mainly enriched in actin binding, regulation of cytokine production, muscle cell promotion, etc. KEGG results revealed that DEGs were mostly related to the chemokine signaling pathway, Toll-like receiver signaling pathway, cytokine-cytokine receiver interaction, and HIF-1 signaling pathway.Cytokine (CK) acts as molecular polypeptides or glycoproteins, mediating intercellular interactions. It is known to have participated in cell growth, maturation, differentiation, tumor growth, and decline [21]. Among many cytokines, IL-1, IL-4, IL-6, IL-8, and M-CSF can directly act on osteoclasts and cause corresponding biological effects. Other cytokines act on osteoblasts, regulate the contact between osteoblasts and osteoclasts, and indirectly produce effects [22]. Chemokines are a member of the CK group, with light molecular weight, which plays a role in mediating leukocyte migration thereby affecting inflammation [23]. Toll-like receptors (TLRs) can identify specific molecular patterns of pathogens in multiple microorganisms, and they are important for innate immunity [24]. Tumor necrosis factor- (TNF-) α can trigger the differentiation of osteoclast precursors into mature osteoclasts secreted by LPS/TLR4 signaling [25]. CX3CL1 chemokine is crucial for OP pathogenesis [26, 27]. HIF-1α signal transduction in B cells is essential for controlling osteoclast production [28].It was found that p53 may play a central role in the development of OP [29]. RBBP4 and DHTKD1 may also be potential regulators of PMOP by interacting with ESR1 and regulating mitochondrial dysfunction, respectively. OSTF1, ADRB1, and NEO1 are involved in PMOP by regulating the balance between bone formation and bone resorption [30]. Six central genes, including EGF, TP53, ICAM1, PRKACA, PXN, and NGF, may be potential targets for the treatment of OP [31]. In this study, 10 key genes Ccl5, Cxcl10, Egfr, Hmox1, Il12b, Ccl7, Tbx21, Xcl1, Pgr, and Itga1 were identified by PPI network analysis.CCL5 is recognized as a downstream gene of NF-κB and locates on mouse chromosome 11 [32]. This chemokine can increase the activity of myoblast migration [33]. CXCL10 is reported to promote osteoclast differentiation [34]. PRMT5 is an activator of osteoclast genesis. PRMT5 inhibition inhibits osteoclast differentiation by downregulating CXCL10 and RSAD2 [35]. As a member of the ErbB family, the epidermal growth factor receptor (EGFR) acts as a transmembrane glycoprotein. When activated, it triggers many downstream signaling pathways involved in regulating the proliferation and differentiation of cells [36]. Transgenic mouse models have shown that EGR2 is expressed in osteoblasts and chondrocytes and is important for skeletal development [37]. The EGFR-induced EGR2 expression is critical for osteoprogenitor maintenance and new bone formation [38].Withdrawal estrogen increased CCR2, CCL7, and CCL2 expression but reduced PMOP in mice after inhibition of CCL7 and CCL2 synthesis by administration of band [39]. Hmox1 (HO-1) has been identified as a new potential therapeutic target against OP [40, 41]. TBX21 gene is a T-box expressed in T cells. This gene is a core feedback loop regulating Th1 and Th2, which can promote the differentiation of Th1 cells, induce Th1 cells to secrete IFN-γ, and participate in cell-mediated immune responses [42]. Studies have shown that progesterone can be used to prevent and treat OP in women [43]. The study showed that EGFR, HMOX1, and PGR were significantly upregulated in estrogen-treated osteoporotic patients, while CCL5, CXCL10, IL12B, CCL7, TBX21, XCL1, and ITGA1 were significantly downregulated. Although most of these genes have been proved to be markers of OP, none of Tbx21, PGR, and Itga1 have fully demonstrated their roles in OP. Therefore, we can further confirm its role in OP in the future, which will provide a new direction for the selection of biomarkers and the exploration of therapeutic targets for OP.This study has some limitations. Firstly, the microarray dataset sample size was not large enough. Secondly, despite several enriched pathways and key DEGs being identified, they have not been verified in vivo and in vitro experiments. Therefore, potential OP-related biomarkers remain to be further determined using a larger sample size.


## 5. Conclusion

In conclusion, this study identified 10 genes significantly associated with OP, providing an essential foundation for investigating the molecular mechanism of OP biomarkers and therapeutic target selection. However, biological studies with larger sample sizes and validation are lacking; so, further experiments are needed to verify the diagnostic ability of these genes for OP before clinical application.



---

*Source: 6726038-2023-01-30.xml*