This article is mentioned in:

Introduction

Peak bone mass (PBM), the highest bone mass obtained during the life of a person, is a critical determinant of osteoporosis (1,2). It has been reported that gaining a high PBM in the first 30 years could effectively prevent osteoporosis (3,4). Circulating monocytes are important cells that serve as early progenitors of osteoclasts and has been reported to be a suitable model for bone-related research (5,6). Deng et al have identified and confirmed five proteins [ras suppressor protein1 (RSU1), gelsolin (GSN), manganese-containing superoxide dismutase (SOD2), glutathione peroxidase 1 (GPX1), and prolyl 4-hydroxylase β subunit (P4HB)] that may involve in differential osteoclastogenesis by a comparative protein expression profiling study and western blot analysis of circulating monocytes in premenopausal females with extremely discordant bone mineral density (7). Lei et al identified three novel risk genes (GBP1, STAT1 and CXCL10) that may be involved in the differentiation of PBM at the monocyte stage (1). Although there have been some studies on the role of circulating monocytes in osteoporosis development, the relationship between circulating monocytes and the pathogenesis of osteoporosis remains poorly understood.

Gibbs sampling, also called Gibbs sampler, is a Markov chain Monte Carlo (MCMC) algorithm for getting a sequence of observations that are approximated from a specified multivariate probability distribution (8). It can be used to identify perturbed pathways and hub genes which may be closely related to occurrence and development of the disease. In this research, we aimed to identify the differential pathways and hub genes in circulating monocytes of healthy Chinese women with high PBM and low PBM using Gibbs sampling, attempting to promote our understanding of the function of circulating monocytes in osteoporosis.

Materials and methods

Data collection and preprocessing

Transcription profiling of human circulating monocytes was obtained from the ArrayExpress database (http://www.ebi.ac.uk/arrayexpress/) and the accession number was E-GEOD-7158. The database consists of 26 samples, including 14 samples from healthy Chinese females with extremely high PBM and 12 samples from healthy Chinese females with extremely low PBM. Followed by preprocessing (9,10), the expression profile was converted from probe level to gene symbol level and the duplicated symbols were removed. The remaining genes were used for subsequent analysis.

Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathway enrichment

The KEGG database was used in this study (11). A total of 287 pathways, including 6,894 genes were downloaded from the KEGG database. The genes in the transcription profile were enriched to the KEGG pathways. Pathways with gene counts ≥5 were screened out for the next analysis.

Construction of path initial state

The average expression values of the genes in the pathways under high and low PBM conditions were calculated, respectively, and considered as the expression of this pathway. High PBM condition was regarded as the initial state of the pathway and low PBM condition was regarded as a priori value of the pathway.

According to the Markov chain theory, we deduced the initial transition probability from the expression of state 1 and state 2. Then state n is derived from state n-1 and the state n is only associated with state n-1. The original state of the system has no connection with the other Markov processes.

Gibbs sampling

Gibbs sampling was implemented to form a Markov chain. First, an empty Gibbs sampling set was defined followed by storing the initial state and a priori value of the pathways into the Gibbs sample empty set. Secondly, a k-dimensional vector (k =280) was initialized randomly and k-1 of these vectors were fixed accompanied by extracting the remaining one vector. The posterior value of this vector was calculated. After repeating this process for k times, a new k-dimensional vector was obtained and defined as state 3. After repeating this whole process 10,000 times, a Markov chain was obtained. The vectors generated after 2,000 times are considered to follow the distribution of the real samples, which is a satisfactory convergence.

Identifying key pathways

After gaining the posterior values of all the pathways, the probability (α) of occurrence of each pathway was calculated using the following formula 1.

Pi represents the posterior value of this pathway in the ith sample.

According to the expression level of the pathway in different states, the P-value of the pathway is calculated by Student's t-test. All the pathways were then ranked according to their P-value. Subsequently, the correction coefficient c was calculated using the following formula 2.

In formula 2, ranki represents the rank of pathway i and n represents the number of samples.

Next the adjusted probability of each pathway was calculated as follows: αadj = α × c. After accomplishing all the adjusted calculation probabilities, the pathways were all ranked according to their αadj value and pathways with αadj ≥0.8 were selected as differentially expressed pathways (12).

Identifying hub genes

After identifying the differentially expressed pathways, we proceeded to find the hub genes in these pathways by Gibbs sampling. The process of finding hub genes was the same with the identification of differentially expressed pathways. Firstly, the genes in the identified pathways were transformed into Markov chains. Secondly, Gibbs sampling was performed to obtain a new Markov chain and the posterior values of these genes. Thirdly, the probability (α) of occurrence of each gene was calculated. Then, the P-values were calculated by t-test according to gene expression differences between high group and low group and the genes were ranked basing on the P-value. Next, the correction coefficient c was calculated and the adjusted probability of each gene was calculated. Finally, the genes with αadj ≥0.8 were selected as hub genes.

Results

Data preprocessing and pathway enrichment

The transcription profile was downloaded from the ArrayExpress database, including 12 samples from healthy Chinese female with extremely low PBM and 14 samples from healthy Chinese female with extremely high PBM. After pretreatment, we obtained a total of 20,514 genes used for the next analysis.

All human pathways (287) were downloaded from the KEGG database and 280 pathways were selected by screening pathways with a gene intersection ≥5.

Key pathways screened by Gibbs sampling

The posterior probabilities of the 280 pathways were obtained by Gibbs sampling and then the probabilities of the 280 pathways were calculated and adjusted. Fig. 1 shows the probability distribution of the 280 pathways. The blue line indicates the threshold value of probability of 0.80. All 280 pathways were ranked according to their adjusted probability and the top five pathways with the highest αadj are listed in Table I. Under the criterion of αadj ≥0.8, neuroactive ligand-receptor interaction (hsa04080) which ranked at first with αadj = 0.986 was screened out as the key pathway. The expression levels of hsa04080 in the high and low groups are shown in Fig. 2. We can identify that the expression level of hsa04080 was obviously decreased in the low group compared to that in the high group.

Figure 1. The probabilities and density distribution of pathways. (A) The distribution probabilities of all human pathways. (B) The density distribution of pathways. The red line shows the initial probabilities, and the green line shows the adjusted probabilities.

Figure 2. The expression levels of neuroactive ligand-receptor interaction pathway in the high and low group.

Table I. The top 5 pathways ranked according to αadj value.

Table I.

The top 5 pathways ranked according to αadj value.

Pathway	P-value	Rank p	R-value	α	αadj
Neuroactive ligand-receptor interaction (hsa04080)	0.002	1	0.996	0.990	0.986
Natural killer cell-mediated cytotoxicity (hsa04650)	0.068	11	0.961	0.756	0.726
Hedgehog signaling pathway (hsa04340)	0.038	3	0.989	0.701	0.694
Basal cell carcinoma (hsa05217)	0.053	7	0.975	0.605	0.590
PI3K-Akt signaling pathway (hsa04151)	0.054	8	0.971	0.481	0.467

Hub genes in the key pathway

By analyzing the key pathway identified previously, we found 309 genes in the neuroactive ligand-receptor interaction. Gibbs sampling was performed again on these pathway genes. The probability distribution of genes are shown in Fig. 3. The blue line represents the threshold value of probability of 0.80. All 309 genes were ranked according to their αadj and Table II shows the top five genes with the highest αadj value. Neurotensin (NTS), tachykinin receptor 3 (TACR3), follicle-stimulating hormone receptor (FSHR) with αadj value ≥0.8 were considered as the hub genes. Of the three genes, TACR3 gained the highest αadj value (0.930). The expression levels of the three genes were analyzed (Fig. 4). The expression level of NTS was reduced in the low group compared to that in the high group (Fig. 4A). The expression level of TACR3 was significantly increased in the low group compared with that in the high group (Fig. 4B). By contrast, the expression level of FSHR was significantly decreased in the low group compared to the high group (Fig. 4C). Heatmap of the three hubgenes is shown in Fig. 5, and the results were consistent with those of Fig. 4.

Figure 3. Probabilities and density distribution of disturbed genes. (A) Distribution probabilities of genes in neuroactive ligand-receptor interaction pathway. (B) The density distribution of genes in the disturbed pathways. The blue line represents the posterior value of 0.8, the red line shows the initial probabilities, and the green line shows the adjusted probabilities.

Figure 4. The expression levels of (A) NTS, (B) TACR3 and (C) FSHR in human circulating monocytes from healthy Chinese women with high or low PBM, respectively. NTS, neurotensin; TACR3, tachykinin receptor 3; FSHR, follicle-stimulating hormone receptor; PBM, peak bone mass.

Figure 5. Heatmap of the identified hub genes NTS, TACR3 and FSHR. Red represents the expression of the hub genes in low group. Blue represents the expression of the hub genes in high group. NTS, neurotensin; TACR3, tachykinin receptor 3; FSHR, follicle-stimulating hormone receptor.

Table II. The top 5 genes in neuroactive ligand-receptor interaction ranked according to αadj value.

Table II.

The top 5 genes in neuroactive ligand-receptor interaction ranked according to αadj value.

Gene	P-value	Rank p	R-value	α	αadj
TACR3	0.065	18	0.942	0.988	0.930
FSHR	0.00028	1	0.997	0.864	0.861
NTS	0.109	33	0.893	0.904	0.807
HCRTR2	0.029	10	0.968	0.819	0.793
TAAR8	0.028	9	0.971	0.760	0.738

[i] TACR3, tachykinin receptor 3; FSHR, follicle-stimulating hormone receptor; NTS, neurotensin

Discussion

Circulating monocytes, which possess the potential to differentiate into osteoclasts and to produce potent pro-inflammatory and anti-inflammatory cytokines, are likely to play vital roles in bone metabolism (5). Previously, the analysis of key genes that may be associated with osteoporosis in circulating monocytes has been conducted by a traditional method based on Package linear models for microarray data of R (13). Herein, an MCMC algorithm based on Bayesian theory (Gibbs sampling) was performed for the first time to identify differentially expressed pathways and hub genes which probably associated with osteoclasts in circulating monocytes.

After Gibbs sampling, neuroactive ligand-receptor interaction (hsa04080) with the highest αadj (0.986) was identified as a differentially expressed pathway in circulating monocytes that may relate to osteoporosis development. A total of 309 genes were found in this pathway and Gibbs sampling was performed again to identify the key genes among these 309 genes. Finally, three genes (NTS, TACR3 and FSHR) with αadj >0.8 were identified as hub genes that may be associated with osteoporosis.

NTS is a neuropeptide composed of 13 amino acids (14) found in the brain and secreted mainly by neurons and astrocytes (15). Later, researchers found that NT was also expressed in the intestine (16). NTS in the central nervous system was associated with Parkinson's disease, eating disorders, schizophrenia, substance dependence and neurodegeneration, sexually transmitted diseases and other neuropsychiatric disorders (17,18). In the intestine, it involved in regulating intestinal function. In the present study, NTS was differentially expressed in circulating monocytes of high group and low group and was identified as a hub gene by Gibbs sampling. However, only few reports have demonstrate the relationship between NTS, circulating monocytes and osteoporosis. TACR3, a member of genes which function as receptors for TACs, are mainly detected in the central nervous system and play vital roles in physiological development and human reproductive system (19,20). Recently, it has been reported that the expression of TACR3 was significantly evaluated in Human Oral Squamous Cell Carcinoma compared with normal epithelium (19). However, the roles of TACR3 in circulating monocytes remain unknown. Further studies of the role of NTS and TACR3 in circulating monocytes may uncover new mechanism and provide new target for osteoporosis.

FSH is a heterodimeric glycoprotein which plays a role in controlling follicle development and steroidogenesis by binding to the FSHR located in ovarian granulosa cells (21). Allan et al reported that FSH could increase bone mass in female mice (21). By contrast, Wang et al demonstrated FSH serum concentrations were increased in postmenopausal women with osteoporosis suggesting that FSH may increase the risk of postmenopausal osteoporosis (22). Additionally, FSH receptors were found in osteoclasts cultured in vitro and may promote osteoporosis in postmenopausal women (23). Furthermore, it has been reported that lacking either FSHβ or FSHR in female mice could protect from bone loss despite hypogonadism (24). In the present study, FSHR was identified as a hub gene in circulating monocytes and its expression level in low group is significantly lower than that in high group, indicating that this gene may closely relate to the development of osteoporosis and it has the potential to be a therapeutic target for osteoporosis.

In conclusion, we identified a differentially expressed pathway (neuroactive ligand-receptor interaction) and three hub genes (NTS, TACR3 and FSHR) by analyzing the expression profiling of human circulating monocytes using Gibbs sampling. These findings provide insights into understanding the role of circulating monocytes in osteoporosis development. However, further experiments are needed to confirm these findings.

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Authors' contributions

YX conceived the study, acquired and analyzed the data, and drafted the manuscript. XPZ helped to explain the results and revised the manuscript. Both authors approved the final study.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References