Based on gene expression data, cellular pathways and protein-protein interactions (PPIs) information, the identification of dysregulated pathways using the pathway-interaction network (PIN) included three steps. In the first step, the microarray genes were mapped to the pathways and the principal component analysis (PCA) method (17) was used to compute the pathway activity for each pathway according to the sum of the expression values of all genes of this given pathway, and seed pathway was selected based on the pathway activities (the pathway of activity score with the maximum change between PMOP and normal groups was regarded as the seed pathway). In the second step, a PIN was constructed that relied on gene expression data, PPIs and cellular pathways, in which each node represented a cellular pathway. In the third step, dysregulated pathways were extracted from the PIN based on the seed pathway and increased classification accuracy, which was measured using area under the curve (AUC) indexing based on five-fold cross validation.

Raw gene expression data were downloaded from the ArrayExpress database offered by European Bioinformatics Institute (http://www.ebi.ac.uk/arrayexpress; accession no. E-MEXP-1618) (15), which was based on the A-AFFY-44 Affymetrix GeneChip Human Genome U133 Plus 2.0 (HG-U133_Plus_2) platform. This data set included 84 iliac bone biopsy samples, comprising 45 patients with OP and 39 healthy patients, which were obtained from postmenopausal women. The inclusion criteria were an average spine T-score >-1 for the normal group, and an average spine T-score <-1 for the PMOP group.

Affymetrix CEL files and probe annotation files were downloaded, and the gene expression profile of E-MEXP-1618 was pre-treated using the Affy package of Bioconductor (18). Briefly, primary data in the CEL format was converted into expression measures, followed by background correction using robust multiarray average (19). Next, quartile data normalization was conducted through quantiles (20), and median polish was used to summarize the expression measures. Finally, a gene expression matrix that included 20,545 genes was obtained following the probes were mapped to gene symbols. Subsequently, expression scores of all genes were standardized.

A total of 1,675 predefined biological pathways were obtained from the REACTOME database (http://www.reactome.org), which is a manually curated open-source database of cellular pathways (21). Pathways comprising too few genes may not have sufficient biological information, whereas pathways having too many genes may be too generic (22). Therefore, the present study extracted a set of pathways by excluding the pathways with <5 genes or >100 genes. Following removal of these pathways, 1,189 background pathways were identified for further analysis. PPIs were downloaded from the STRING database (version 9.1; http://string-db.org) (23). The original PPI data set included 787,896 interactions among 16,730 genes. The STRING database uses confidence scores to measure how likely an association will appear. As a result, with the goal of minimizing the ambiguity, only the interactions with confidence score >0.2 were selected to construct the original PPI network. Taking the intersection of the original PPI network and the microarray data, a background PPI data set containing 14,917 proteins and 449,833 PPIs was identified and used for further analysis.

Using gene expression profiles, PPI information and pathway data, a PIN was constructed, in which each node represented a pathway and an edge was laid between two pathways when they met one of two conditions; if not, then the edges were discarded. The first condition was that these two pathways had to share at least one gene, and at least one of the common genes between two pathways had to be a DEG between the PMOP and normal groups. In the present study, DEGs were identified using Student's t-test with a cutoff of P<0.01. The second condition was that the two genes that coded interacting proteins used to lay an edge between the two pathways had to be highly co-expressed, with an absolute value of Pearson's correlation coefficient (|PCC|) >0.8; if not, then the edge between two pathways would be rejected. Based on these requirements, a PIN was constructed. It is considered that if a network is too big that a considerable number of key genes and interactions are probably ignored (24). As a result, to reduce the complicated network, the scores of a pair of pathways in the PIN were calculated; defined as the summation of the absolute values of PCC for the PPIs in every two pathways. The top 5% pathway interactions were chosen to construct an informative PIN for the detection of dysregulated pathways.

PCA is a dimension reduction approach that has been widely used in gene expression analysis (17); it is able to efficiently depict the internal structure of high-dimension data by preserving the variance within the data while converting the data into low-dimension space (25). In the present study, PCA was used to calculate pathway activity according to gene expression data from each pathway. Notably, all the genes were mapped to background pathways, and only those genes that aligned to the background pathways were kept for PCA analysis. The activity of each pathway (that is, the summary of the expression values of all genes in the given pathway) was calculated using the PCA method. The activity score for the corresponding pathway was determined as the first principal component from PCA. The activity score for a pathway between the experimental and the control groups was different, and the difference may demonstrate its correlation to the disease; that is, the greater the difference was, the more relevant this pathway was to the disease. Therefore, in the present study, the pathway with the maximum change of activity score between PMOP and control groups was selected as the seed pathway.

Support vector machines (SVMs) were used to extract the dysregulated pathways. Briefly, an individual pathway that best distinguished disease from normal state was identified as the first pathway biosignature (hereafter called the seed pathway), and the second pathway that was able to integrate with the seed pathway to obtain better classification accuracy was extracted from the pathways that interacted with the seed pathway in the PIN. The process was repeated, and new pathways were added to the obtained pathway set until no more pathways could be assembled into the pathway set to enhance classification results. Classification accuracy was measured using the AUC index, based on five-fold cross validation. In an attempt to obtain robust biomarkers, five-fold cross-validation was repeated 100 times. The mean value was collected as the final result.

Based on the cutoff criteria of P<0.01, a total of 544 DEGs between the PMOP group and the normal group were identified. The top 20 DEGs, in ascending order of P-values, are provided in Table I. The identified DEGs were used to select the interactions for establishing the PIN, as only the interactions in the background PPI data were able to meet at least one of the two conditions and were kept to build the PIN.

Based on the pre-defined threshold, a PIN that was involved in 54,517 pathway-pathway interactions was constructed. As very large networks may easily neglect a few significant interactions (23) they should be reduced in size. In the present study, the interactions with low |PCC| scores were removed, and the top 5% of all interactions were adopted for further analysis. A network comprising the top 5% of interactions was defined as the informative PIN (Fig. 1). The informative PIN included 2,725 interactions between 1,189 pathways and was used to identify dysregulated pathways; specific information about the 1,189 pathways is not provided here. Based on the informative PIN, the strengths among the pathways were distinguishable. The weight of a pathway-pathway interaction was defined as the summation of |PCC| scores of all genes, and interactions having greater weight scores may be more important for PMOP compared with other interactions. The weight scores of the 2,725 interactions ranged between 95 and 323. Notably, in the informative PIN, there were only 11 pairs of pathway interactions with weight scores >300.

As there were differences between pathways in the informative PIN, how to assess the significance of each pathway and extract the significant pathway in the PIN was a challenge. The activity score for the 1,189 pathways was calculated using the PCA method to evaluate the significance. It was revealed that the pathways had different changes to their activity score between PMOP and normal groups (Fig. 2). The pathway with the maximum change in activity score between PMOP and normal groups was regarded as the seed pathway, which, in the present study, was identified as ‘mitotic prometaphase’.

Starting with the seed pathway, a pathway set was identified by adding pathways to increase the classification accuracy, and this procedure continued until the AUC index decreased. Using ‘mitotic prometaphase’ as the seed, a pathway set was identified that was involved in 8 dysregulated pathways with and AUC score of 0.85. This pathway set included ‘mitotic prometaphase’, ‘resolution of sister chromatid cohesion’, ‘transport of mature mRNA derived from an intron-containing transcript’, ‘interactions of Vpr with host cellular proteins’, ‘metabolism of non-coding RNA’, ‘regulation of HSF1-mediated heat shock response’, ‘mRNA splicing’ and ‘mRNA splicing-major pathway’. The effective classification ability of this method indicated that the detected dysregulated pathways may be utilized as robust signatures. A network of these eight dysregulated pathways is provided in Fig. 3.

In addition, the genes enriched in the detected dysregulated pathways were compared with the identified DEGs (Table II). It was observed that only a small fraction, between 8.4 and 21.2%, of the genes in the extracted dysregulated pathways overlapped with the DEGs. This phenomenon further demonstrated that a pathway as an entity may be better than a single gene at diagnosing complicated diseases, even if the genes enriched in the pathway were not differentially expressed.