This article is mentioned in:

Introduction

Osteoporosis refers to a systemic skeletal disease characterized by a low bone mass and deterioration of the microstructure, followed by an increase in bone fragility and fracture (1). As one main type of osteoporosis, postmenopausal osteoporosis (PO) results from a scarcity of endogenous oestrogen in postmenopausal women (2). Clinically, the two major manifestations for PO patients are low bone mass density and fracture, and even the occurrence of vertebral compression fractures during routine activities without a specific fall or injury (3). Consequently, PO imposes a great burden on life and work both for the individual and society. Thus, an effective prevention, diagnosis and treatment before the occurrence of fractures for this disease is urgently needed (4).

Microarray has become a powerful technology that simultaneously monitors the expression levels of thousands of genes, and provides many basic applications, such as functional genomics, molecular pathway modeling and tumor classification (5). The most commonly in use is the comparison of gene expression discrepancy between two different statuses (normal controls versus disease and treated versus untreated samples). Furthermore, under such an experimental setup, how to investigate a gene whose expressed value is remarkably different between the two status is a dominate challenge. A variety of methods have been proposed to detect significant genes across controls and cases in the past few years (6,7). However, it has been revealed that the most significant genes obtained from distinct research are typically inconsistent for a specific disease because of batch effects and method difference (8).

Therefore, to avoid inconsistent results and enhance the feasibility and confidence of the results, an integrated method was proposed in the present study. The integrated method combined three commonly used methods, including the convergent evidence (CE) method, rank product (RP) algorithm and combing P-values, to identify informative genes for PO. Subsequently, hub genes were identified by accessing topological analysis of mutual information network (MIN) which were constructed based on the informative genes. Moreover, hub gene functions were determined by Gene Ontology (GO) enrichment analysis for informative genes.

Materials and methods

Preparing gene data

A gene expression dataset [E-MEXP-1618 (9,10)] was downloaded for PO from the online ArrayExpress database (http://www.ebi.ac.uk/arrayexpress/). In brief, E-MEXP-1618 consisted of 39 normal controls and 45 PO samples, and deposited on the A-AFFY-44-Affymetrix GeneChip Human Genome U133 Plus 2.0 [HG-U133_Plus_2] Platform. Subsequently, duplicated or invalid probes in the data were removed, and the retained probes were converted into gene symbols. Ultimately, 20,544 genes were obtained for subsequent exploitation.

Ranking genes by multiple algorithms

Currently, an effective method to investigate and prioritize candidate genes in complex diseases is integration of gene data from multiple evidence layers (11). Thus, to recognize informative genes generated from multiple independent lines of investigation for PO patients, we integrated 3 commonly used approaches, CE method, RP algorithm and combining P-values.

CE method

CE method, a variant of the famous PageRank algorithm, is implemented in the GenRank Bioconductor package (http://bioconductor.org/packages/GenRank). This method aggregates ranks of genes using a weighted vote counting method that the genes were ranked according to the self-importance and the number of each evidence layer was detected (12). In addition, the algorithm permits us to integrate evidence from studies where P-values and effect sizes are unavailable or the gene sets varies across research. During this process, the weighted arithmetic mean of each gene, a CE score, was calculated as the given formula (13):

of which CE(Li) referred to the self-importance of gene g in evidence layer i, while n(Li) represented the number of genes within evidence layer i.

RP algorithm

To the best of our knowledge, RP method has earlier been widely used to carry out differential expression and meta-analysis in microarray-based gene expression datasets (5,14). Furthermore, it is also a simple, but powerful to rank genes that are consistently ranked high in replicated microarray tests (5). Thus, we used the RP algorithm to identify consistently highly ranked genes across various evidence layers. Supposing that RPi, g was the position of gene g in the list of genes in the ith replicate or layer among k replicates sorted by decreasing fold-change (FC), the RP value for a gene g was computed as following:

Consequently, all genes in the gene lists were sorted according to their RP values. Note that genes possessing the smallest RP values are the most fascinating candidates and researchers could choose some of them for further investigation.

Combining P-values

Estimation of P-values from the same null nypothesis was the precondition of this method. When a null assumption was examined on distinct sets of data, the individual P-values may effectively be integrated to supply a more powerful test than any individual one (15). Each P-value abides by a uniform distribution, in case the null hypothesis is true, and the statistical model is correct. The Fisher's exact test was utilized to combine the P-values, which is a statistical test to compare binary outcomes between two groups (16). In consequence, we obtained P-value for each gene and corrected them by false discovery rate (FDR) implemented in Benjamini-Hochberg (BH) method (17).

Extracting informative genes

As mentioned above, based on 20,544 genes in the gene data, we obtained three different gene lists ranked by CE, RP and combining P-values, respectively. Subsequently, we selected the top 100 genes from the three gene lists separately. Taking intersections between any two of them was a good way to overcome the inconsistent results from the top 100 genes. These common genes might be more significant than the others, and thus we defined them as the informative genes for PO patients. Analyses were conducted based on the informative genes to investigate hub genes and gene functions in the progression of PO.

Constructing network and topological analysis

Constructing network

To investigate the co-operations and interactions between any two of informative genes, a MIN was constructed based on the context likelihood of relatedness (CLR) algorithm. CLR algorithm is an information-theoretic approach that applies a generalization of the pairwise correlation coefficient called mutual information for each pair of genes (18). For a pair of genes × and y, an edge score, Sxy, was derived (19):

of which

where μMx represents the sample mean of the empirical distribution of the values Mxy, and σMx refers to the standard deviation Mxy. Similarly, Sy was calculated. Finally, Cytoscape software (http://www.cytoscape.org/) was used to visualize MIN for PO (20).

Topological analysis

For purpose of deeply investigating biological significance of nodes in the MIN, topological centrality analysis of degree index was conducted. Degree quantifies the local topology of each gene by summing up the number of its adjacent genes (j) and gives a simple count of the number of interplays of a given node (21). Besides, the nodes on the top of degree distribution (≥10% quantile) in the MIN were described as hub genes for PO.

Gene functional enrichment analysis

The Biological Networks Gene Ontology tool (BiNGO), is a tool to identify which GO terms are significantly overrepresented in a set of gene lists, and utilizes Cytoscape's versatile visualization environment to generate an intuitional and customizable visual representation of the results (22). Thus, we applied the BiNGO to perform the GO enrichment analysis for informative genes in this study. The GO terms with P<0.001 were considered to be predominant functional themes, and screened for mapping to the GO hierarchy. Most importantly, the term at the top of P distribution (≥10% quantile) among the significant gene sets were defined as hub gene functions of PO.

Results

Informative genes

In the present study, three commonly used methods of ranking genes were implemented to prioritize significant genes for PO patients. As a result, we obtained three gene lists dependent upon their CE scores, RP scores and P-values, respectively. The top 100 genes which might be more important than others were selected from each gene list for further analysis. Although there were differences among the three kinds of top 100 genes, we found that a certain number of genes were common. As shown in Fig. 1, 82 common genes from two or even all three methods were identified, of which 40 genes were common to the three methods. All 82 common genes were defined as informative genes of PO patients, where 40 genes were attributed to the numbers from 1 to 40, such as LYZ, LCN2 and EEF2 (Table I).

Figure 1. Venn diagram for common genes between any of the two methods.

Table I. Informative genes.

Table I.

Informative genes.

No.	Gene	No.	Gene	No.	Gene	No.	Gene
1	LYZ	22	RPL23A	43	CAMP	64	FTH1P5
2	LCN2	23	RPL37	44	UBC	65	MMP9
3	EEF2	24	PRG2	45	RPL11	66	HLA-A
4	COX1	25	RPL32	46	AHSP	67	RPS5
5	HBB	26	RNASE3	47	NCOA4	68	TALDO1
6	UBB	27	ELANE	48	AZU1	69	RPL3
7	S100A9	28	RNASE2	49	TUBA1C	70	CA2
8	TMSB4X	29	RPS2	50	RPL8	71	CEACAM8
9	LTF	30	GAPDH	51	RPL12	72	PFN1
10	RPL39	31	RPS17	52	RPL10A	73	HLA-B
11	S100A12	32	RPL27	53	RPS25	74	PF4
12	RPS4X	33	HMGN2	54	HBD	75	RPL35
13	ND2	34	RPLP0	55	RPL34	76	COX2
14	DEFA4	35	RPS3	56	RPS7	77	IGLL3P
15	TPT1	36	CLC	57	MPO	78	BLVRB
16	RPS18	37	CTSG	58	UBA52	79	ANP32B
17	ACTB	38	RPL19	59	TMSB10	80	HLA-C
18	RPL9	39	OAZ1	60	RPL6	81	RPL41
19	FTL	40	ND4	61	RPL29	82	PRG3
20	HBM	41	RPL24	62	MNDA
21	RPL30	42	PPBP	63	FAU

Hub genes

Since genes often work together in complex disease rather than individually, a MIN based on the informative genes of PO was constructed using the CLR algorithm, as illustrated in Fig. 2. There were 82 nodes and 1,741 edges in the MIN, which suggested that all of the informative genes were mapped to the PO-associated network. Furthermore, topological centrality of degree analysis was conducted to screen crucial genes in the MIN, and we denoted the nodes at the top 10% of degree distribution as hub genes for PO. In consequence, a total of 8 hub genes were identified, including PFN1 (Degree = 63), EEF2 (Degree = 54), S100A9 (Degree = 53), LCN2 (Degree = 52), ACTB (Degree = 51), FAU (Degree = 50), RPS25 (Degree = 50) and CTSG (Degree = 50).

Figure 2. MIN of PO. There were 82 nodes and 1,741 edges in the MIN. Nodes referred to the informative genes, especially the yellow stood for hub genes. Besides, edges were the interactions between any two nodes. MIN, mutual information network; PO, postmenopausal osteoporosis.

Hub gene functions

GO enrichment analysis was conducted on the informative genes by the BiNGO plugin in order to investigate significant gene functions in the progression of PO. The result displayed that 49 GO BP terms were determined when setting the threshold as P<0.001 (Fig. 3). The darker the node, the more significant the GO term for PO was. In addition, we ordered the 49 terms in descending order of their P-values, and defined the top 10% as hub terms enriched by informative genes. Specifically, the hub gene functions were comprised of translational elongation (P=7.317E-49), translation (P=9.21E-33), cellular macromolecule biosynthetic process (P=3.37E-16), macromolecule biosynthetic process (P=5.92E-16), and gene expression (P=8.44E-14).

Figure 3. Significant GO BP terms with P<0.001 enriched by informative genes of PO. The size of a node indicates the number of genes enriched in this term. The color represented its P-value, the smaller the P, the darker the node is. The arrows stand for progression of BP terms. GO, Gene Ontology; BP, biological process; PO, postmenopausal osteoporosis.

Discussion

A better understanding of pathological mechanism underlying PO may give crucial insights for prevention and treatment in the progression of this disease. Hence, in this study, we predicted hub genes and gene functions for PO by integrating CE method, RP algorithm and combing P-values approach. In brief, the CE method integrates single-case assessment data by converging data gathered across multiple settings, or measures, thereby supplying an overall criterion-referenced outcome (23). Importantly, a conceptually similar gene-level integration of CE method had been smoothly applied to prioritize candidate genes in neuropsychiatric diseases (24). In addition, the RP algorithm provides a powerful statistics for defining significant genes in microarray tests, and needs only a few well-justified hypothesis on the data because of its non-parametric nature (5). Besides, the combing P-values approach is a typical and efficient method to rank genes by their P-values. Above all, we merged the three methods to enhance the feasibility and precision of our research. Consequently, we obtained 82 informative genes for PO patients.

In general, a gene is inclined to interact with other genes, and the co-operated genes often actively participate in the same biological processes. Moreover, network-based approach is able to extract informative and significant genes dependent upon bio-molecular networks and provides a quantifiable characterization of the molecular networks that characterize the complicated interplays and the complex interwoven associations that administer cellular functions (25,26). Thus, a network was constructed on the informative genes, termed with MIN. The hub genes were extracted from the MIN by accessing the degree centrality analysis. As a result, a total of 8 hub genes were identified, PFN1, EEF2, S100A9, LCN2, ACTB, FAU, RPS25 and CTSG. Taking S100A9 as an example, S100A9 (S100 calcium binding protein A9) is a member of the S100 family of proteins containing 2 EF-hand calcium-binding motifs. S100 proteins involve in many biological processes, including proliferation, apoptosis, differentiation, inflammation, Ca2+ homeostasis, as well as migration/invasion by interplays with various target proteins such as receptors, enzymes, transcription factors, cytoskeletal subunits and nucleic acids (27). Specifically, calcium supplementation in PO markedly changed bone mineral and organic matrix quality (28). In addition, S100A9 expressed earlier than the appearance of neutrophilics in the airway of neutrophilic inflammation, and directly activated inflammasome in the airway in asthma (29,30). Due to its expression in inflammation, S100A9 had been indicated to play a significant role in cell proliferation (31). Similarly, the significant role of S100A9 in PO may also attribute to its activities in cell proliferation. Most importantly, this is the first clue of the correlation between S100A9 and PO patients.

Moreover, functional gene sets enriched by informative genes were identified based on the GO enrichment analysis. Consequently, a total of 5 hub gene functions were gained, including cellular macromolecule biosynthetic process, translational elongation, translation, gene expression and macromolecule biosynthetic process. Interestingly, we found that 2 of 5 hub gene functions belonged to the biosynthetic process. Macromolecule biosynthetic process refers to the chemical reactions and pathways resulting in formation of a macromolecule, any molecule of high relative molecular mass, the structure of which essentially comprises the multiple repetitions of units derived, actually or conceptually, from molecules of low relative molecular mass. Mackiewicz et al revealed that macromolecule biosynthesis was a key function of sleep (32). Furthermore, it had been demonstrated that the most obviously involved processes for breast cancer were translation, cellular biosynthesis, and macromolecular biosynthesis (33). The change of hormonal readiness was the common point across breast cancer and PO. Thus we might infer it was the medium for the gene function macromolecular biosynthesis process and PO patients.

In conclusion, we identified 8 hub genes and 5 hub gene functions for PO by an integrated method. Our results provide potential biomarkers for prevention and therapy of PO patients, and further help uncover the pathological mechanism underlying this progression. Whereas the validations of these genes have not finished, and future study should focus on this aspect.

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

HoC and PT designed the study. LZ collected the data. HuC performed the statistical analysis. WZ preparated the figures. HoC, QZ and XL drafted the manuscript and interpreted the data. YG is accountable for the accuracy or integrity of the work. YG and PT revised the manuscrript.

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