A novel method to identify hub pathways of rheumatoid arthritis based on differential pathway networks.

The aim of the current study was to identify hub pathways of rheumatoid arthritis (RA) using a novel method based on differential pathway network (DPN) analysis. The present study proposed a DPN where protein‑protein interaction (PPI) network was integrated with pathway‑pathway interactions. Pathway data was obtained from background PPI network and the Reactome pathway database. Subsequently, pathway interactions were extracted from the pathway data by building randomized gene‑gene interactions and a weight value was assigned to each pathway interaction using Spearman correlation coefficient (SCC) to identify differential pathway interactions. Differential pathway interactions were visualized using Cytoscape to construct a DPN. Topological analysis was conducted to identify hub pathways that possessed the top 5% degree distribution of DPN. Modules of DPN were mined according to ClusterONE. A total of 855 pathways were selected to build pathway interactions. By filtrating pathway interactions of weight values >0.7, a DPN with 312 nodes and 791 edges was obtained. Topological degree analysis revealed 15 hub pathways, such as heparan sulfate/heparin‑glycosaminoglycan (HS‑GAG) degradation, HS‑GAG metabolism and keratan sulfate degradation for RA based on DPN. Furthermore, hub pathways were also important in modules, which validated the significance of hub pathways. In conclusion, the proposed method is a computationally efficient way to identify hub pathways of RA, which identified 15 hub pathways that may be potential biomarkers and provide insight to future investigation and treatment of RA.

Introduction

Rheumatoid arthritis (RA) is a common, complex and long lasting autoimmune disorder that originates in the joints (1). It typically manifests with signs of inflammation, the affected joints becoming swollen, warm, painful and stiff, particularly, early in the morning upon waking or following prolonged inactivity (2). However, the underlying mechanism behind RA remains to be fully elucidated. Increased understanding of the immune mechanisms has led to the development of a considerable number of novel therapeutic agents that alter the progression of the disease and reduce mortality (3). Previous genetic studies on RA, including recent genome wide association studies, have identified 32 risk loci among individuals of European ancestry, including major histocompatibility complex class II DR β 1 (HLA-DRB1), protein tyrosine phosphatase non-receptor type 22 (PTPN22) and other loci with shared autoimmune associations (4,5). Meanwhile, peptidyl arginine deiminase type IV (PADI4) has been identified as a major risk factor in people of Asian descent; however, not in those of European descent (6). Few studies have investigated other biological markers that are associated with RA aside from genetic loci, such as the signaling pathways involved.

Pathway analysis has become a common method to gain insight into the underlying biology of genes and proteins, as it reduces the complexity and increases explanatory power fo analysis (7). However, previous studies primarily focused on identifying the altered pathways between normal and cancer groups, or the common genes between two pathways, whereas investigation of the differential interactions between two pathways across disease samples and normal samples is infequent (8,9). Additionally, network-based methods provide more stable and effective measures to investigate functions of genes in certain diseases, and and to understand the connections betweeen different genes (10). Scoring pathways by evaluating the coherency of gene expression changes and combining gene expression quantification over multiple datasets may a lead to novel insights in pathway-associated studies.

Therefore, the present study proposed a novel method to identify hub pathways of RA based on differential pathway network (DPN), which investigated the differential interactions between pathways in RA. In order to perform this, pathway data was obtained from background protein-protein interaction (PPI) networks and the Reactome pathway database. Subsequently, the differential pathway interactions were extracted from the pathway data and gene-gene interactions were built randomly. The differential pathway interactions were visualized using Cytoscape to construct a DPN. Modules of DPN were mined according to ClusterONE. Topological analysis of DPN was conducted to identify hub pathways that may be potential targets for treatment of RA.

Materials and methods

Gene expression data

The gene expression profile of RA with access no. E-GEOD-45291 (11) was obtained from a public functional genomics data repository ArrayExpress (http://www.ebi.ac.uk/arrayexpress/). E-GEOD-45291, which was deposited on A-GEOD-13158- [HT_HG-U133_Plus_PM] Affymetrix HT HG-U133+PM Array Plate Platform (Affymetrix, Inc., Santa Clara, CA, USA), consisted of 378 RA samples and 20 normal controls. The probe-level dataset in CEL file format was converted into expression measures and a total of 7,352 genes were identified in the expression profile.

Pathway data

In the present study, global PPI interactions were downloaded from the Search Tool for the Retrieval of Interacting Genes/Proteins database (http://string-db.org/), which has catalogued a total of 787,896 PPI in human. By mapping the gene expression profile E-GEOD-45291 on the global PPI network, the present study identified a background PPI network with 6,666 nodes and 196,304 interactions.

Information from gene sets representing different biological pathways was downloaded from the Reactome pathway database (www.reactome.org) and a total of 1,675 pathways were obtained. The Reactome database is an online-curated resource for human pathway data and provides infrastructure for computation across biological reaction networks (12). The number of interacting genes between each pathway and the background PPI network were calculated, and pathways with intersections <5 were discarded in order to improve pathway stability. Pathways with a small number of genes are more easily understood; therefore pathways with gene sets >100 were also excluded (13). A total of 855 pathways were reserved as pathway data for further analysis. An ID for each pathway was assigned in ascending order.

Pathway network construction

To evaluate the interactions between different pathways, the present study used the genes in each pathway to construct random gene interactions. When the same gene interaction was established between two pathways, it was considered that the two pathways were correlated. In order to determine the strength of pathway interactions, Spearman correlation coefficient (SCC) (14) was used for intersected interactions in two pathways. Their absolute differences in RA and normal controls were calculated, and the mean value of absolute differences for all intersected interactions was denoted as the weight between two pathways. Taking pathway 1 for example, genes enriched in pathway 1 and pathway 2 were used to build gene interaction and then these interactions were integrated with the background PPI network, the intersections were considered to be pathway interactions between pathway 1 and pathway 2. Subsequently, the pathway intersections under RA and normal conditions were weighted by SCC. The SCC of a pair of interactions (x and y), was defined as:

Where s was the number of interactions, g(x, i) or g(y, i) was the expression level of interaction x or y in the pathway i under a specific condition (RA or normal), ¯g (x) or ¯g (y) represent the mean expression level of interaction x or y. Calculating the absolute difference of an interaction in RA and normal conditions, the mean value of absolute differences of all interactions was defined as the weight value between pathway 1 and pathway 2. If there was no intersection between pathway 1 and pathway 2, the two pathways had no interaction.

Subsequently, gene interactions were constructed based on genes in pathway 1 and pathway 3, the aforementioned method used to determine the interaction and weight value between pathway 1 and pathway 3. Therefore, all interactions and their weight values for any two pathways were obtained and interactions of weight value >0.7 were identified as differential interactions and used to construct a DPN which was visualized by Cytoscape (15).

Topological analysis

To further investigate the functions and importance of pathways in the DPN, the biological importance of pathways was characterized using indices of topological analysis. Degree quantified the local topology of each gene, by summing the number of its adjacent genes (16). This provides a simple count of the number of interactions of a given node. The degree of a pathway was the sum of its adjacent pathways. The pathways at the top 5% of degree distribution (the ≥95% quantile) in the significantly perturbed networks were defined as hub pathways. The degree D(v) of a node v was defined as:

Modules mined from DPN

The emergence of high-throughput techniques for inferring protein interactions on a large-scale has fueled the development of computational techniques to systematically mine for potential complexes from the interaction networks (17). Therefore, the present study mined potential complexes from the network of pathway interaction using computational techniques, and the resulting complexes were defined as modules. The present study used ClusterONE (18), a method for detecting overlapping pathway complexes from weighted differential interactions based on seeding and greedy growth, was implemented to identify the modules of DPN. It used a cohesiveness measure to determine how likely is it for a group of pathways to form a complex, based on the weight of the interactions within the group and with the rest of the network (19).

In the first step, ClusterONE identified seed pathways and greedy growth them into groups with high cohesiveness. When the greedy growth for a group could not progress any more, the next seed pathway was selected and the procedure was repeated until no more seed pathways remained. In the second step, ClusterONE detected highly overlapping cohesive groups and merged them into potential complex candidates or modules. Furthermore, ClusterONE allowed the identification of overlapping complexes if each of the merged groups represented individual complexes that shared pathways. The modules that met the following thresholds: Size ≥20, density ≥0.3 and overlap ≥0.5) were selected as significant modules for the RA DPN.

Results

DPN construction

A total of 855 pathways were identified to establish pathway-pathway interactions and 295,509 interactions between pathways were obtained. Each pathway interaction was assigned a weight value. Fig. 1 presents the distribution of weight values. It was determined that weight values of the majority interactions ranged between 0 and 0.4, particularly 0.1–0.2. This low weight value may indicate that the two pathways had low interaction strength and there was little change in the two pathways between the RA and normal controls; therefore, interactions of with weight values <0.7 were filtered out, and the remaining pathway interactions were selected for further analysis. When inputting these differential interactions into the Cytoscape software, a DPN with 312 nodes and 791 interactions was visualized (Fig. 2). Nodes represented pathways, and edges represented pathway interactions.

Figure 1.

Weight distribution of pathway interactions. (A) Complete data and (B) zoomed (0.7–1.1 weight).

Figure 2.

Differential pathway network of rheumatoid arthritis. Nodes represent pathways and edges show the interactions between pathways. The pink vertices denote hub pathways with degree distribution in the top 5%. The color bar represents the association between color and weight, where the darker the color was the greater the weight.

Hub pathways

In order to determine the functions and importance of 312 pathways, topological analysis of degree for the DPN was conducted. The degree of a pathway was the sum of all connected pathways, and pathways with top 5% degree distribution in the pathway network were defined as hub pathways. Fig. 2 presents hub pathways with pink vertices, and the detailed degrees for 15 hub pathways are presented in Table I. Each pathway was represented by its matching ID number. Heparan sulfate/heparin-glycosaminoglycan (HS-GAG) degradation with degree 35 (ID: 336), HS-GAG metabolism (ID: 325) with degree 33, keratan sulfate degradation (ID: 382) with degree 33, DNA-binding transcription factor RAP1 signaling (ID: 583) with degree 28 and interleukin-3, 5 and granulocyte-macrophage colony-stimulating factor signaling (ID: 367) with degree 26 were the top 5 significant hub pathways of RA based on the DPN. Notably, the top two hub pathways, HS-GAG degradation and HS-GAG metabolism, were both associated with HS-GAG, this may indicate that HS-GAG may be important in RA.

Table I.

Hub pathways of rheumatic arthritis based on differential pathway network.

ID	Pathway name	Degree
336	HS-GAG degradation	35
325	HS-GAG metabolism	33
382	Keratan sulfate degradation	33
583	DNA-binding transcription factor RAP1 signaling	28
367	Interleukin-3,5 and granulocyte-macrophage colony-stimulating factor signaling	26
294	Glucagon signaling in metabolic regulation	24
67	Aquaporin-mediated transport	23
836	Vasopressin regulates renal water homeostasis via aquaporins	23
732	Sphingolipid de novo biosynthesis	22
295	Glucagon-like peptide-1 regulates insulin secretion	21
71	Assembly of the human immunodeficiency virus virion	21
383	Keratan sulfate/keratin metabolism	20
376	Iron uptake and transport	18
128	Clathrin derived vesicle budding	17
309	Golgi-associated vesicle biogenesis	17

HS-GAG, heparan sulfate/heparin-glycosaminoglycan.

Modules

ClusterONE was used to mine modules, which had similar biological processes and functions from the DPN. A total of 4 modules were identified with thresholds of size ≥20, density ≥0.3 and overlap threshold ≥0.5 (Fig. 3). There were 7, 6, 5 and 1 hub pathways in Module 1, Module 2, Module 3 and Module 4, respectively. The weight values of pathway interactions and edges in Module 3 were higher compared with the other modules. In addition, pathways 336 (HS-GAG degradation), 383 (keratan sulfate/keratin metabolism), 325 (HS-GAG metabolism) and 376 (iron uptake and transport) participated in more than one module. For Module 4, only one hub pathway (71, assembly of the HIV virion) was involved. Based on the ClusterONE analysis, hub pathways were important for modules, suggesting that they may have potential roles in the progression of RA.

Figure 3.

Modules extracted from differential pathway network. (A) Module 1. (B) Module 2. (C) Module 3. (D) Module 4. Nodes represent the pathways and the edges represent the interactions of pathways. The pink vertices denote hub pathways with degree distribution in the top 5%. The color bar represents the association between color and weight, where the darker the color was the greater the weight was.

Discussion

The present study proposed a novel method to identify hub pathways of RA using DPN, which was composed of differential pathway interactions based on a background PPI network, the Reactome pathway database and gene expression profiles. It was determined that 15 hub pathways, including HS-GAG degradation, HS-GAG metabolism and keratan sulfate degradation, were significant pathways involved in RA. Modules of DPN were mined depending on which ClusterONE and hub pathways they were involved into validate the feasibility of this novel approach for identifying hub pathways involved in RA.

Diagnostic or prognostic markers are usually obtained by identification of the most significant differentially expressed genes (DEGs) in high-throughput case-control studies of a disease. However, a previous study determined that the most significant DEGs obtained from different studies for a particular disease are frequently inconsistent (20). To overcome this problem, significant genes and biological processes for disease-association may be evaluated using a network strategy, such as PPI networks (21). A network may provide significant instructions for mining unknown connections in incomplete networks. Although the data of large-scale protein interaction has accumulated with the development of high throughput testing technology, a certain number of significant interactions are not tested, such as key genes in significant pathways (22). This may be resolved to some extent by using pathway-associated networks, such as DPN (23). Therefore, the present study proposed a novel method to construct DPN, which determined differential pathway interactions across RA patients and normal controls and identified hub pathways based on DPN. These findings were validated using ClusterONE modules and determined that the method was an efficient and feasible approach.

From the 15 hub pathways identified, HS-GAG degradation, HS-GAG metabolism and keratan sulfate degradation were significant. HS-GAG is a member of the glycosaminoglycan family and consists of a variably sulfated repeating disaccharide unit, the most common one (50% of the total) is glucuronic acid associated with N-acetylglucosamine (24). In addition, GAG and HS have a similar molecular structure, and the polymer chains are composed of repeating disaccharide units of glucosamine and hexuronic acids are sulfated at various positions (25). It has been previously reported that HS-GAG may bind to a core protein and regulate various biological processes, including angiogenesis, blood coagulation and tumor metastasis (26). Li et al (27) revealed that heparanase activity was elevated in the synovial fluid of patients with RA, indicating that heparanase may be a reliable prognostic factor for RA progression and a potential target for RA treatment. Taking into account that heparanase is responsible the degradation of heparin (28) and HS (29), it is possible that degradation of macromolecular heparin in cells may modulate the release of granule proteases that are involved in inflammation (30). Increased heparanase activity in activated inflammatory cells may accelerate the turnover of HS production and thereby induce alterations in HS structure, such as over-sulfation, the increased sulfation of HS induced by upregulated heparanase may facilitate its significant role in disease (31). A previous study suggested that degradation of HS-GAG was required to maintain a natural turnover of GAGs, which accumulate, rather than being broken down by degradative enzymes in RA (32). Degradation of the HS side chains represents an important mechanism underlying chronic inflammation of the arthrosis and associated tumorigenesis (33). Furthermore, heparanase is expressed by the vascular endothelium at the site of inflammation, which results in degradation of the subendothelial basement membrane and subsequent vascular leakage (34). Therefore, RA, as a common autoimmune disease, may be closely associated with HS-GAG closely.

Belcher et al (35) demonstrated that GAG concentration in knee synovial fluid was reduced in RA, and that altered concentrations of chondroitin sulfate and keratan sulfate were be detected in RA. Keratan sulfate is one of several sulfated glycosaminoglycans (structural carbohydrates) that have been identified, in the cornea, cartilage and bone, and consists of large, highly hydrated molecules, which acts as a cushion to absorb mechanical shock in joints (36). The use of circulating keratan sulfate as a marker of metabolic changes of cartilage proteoglycan in juvenile idiopathic arthritis indicated that treatment which modified inflammation simultaneously did not contribute to total regeneration of articular matrix components, and signalized the need for further treatment (37). Additionally, keratan sulfate suppressed cartilage damage and ameliorated inflammation in an experimental arthritis mouse model (38). Previous studies have not identified a correlation between keratan sulfate degradation and RA, however, the present study indicated that keratan sulfate was important for RA progression.

In conclusion, the proposed method used in the current study was computationally efficient to identify hub pathways of RA and 15 hub pathways, including HS-GAG degradation, HS-GAG metabolism and keratan sulfate degradation were identified. These pathways may be potential biomarkers for RA and provide insights for the future study and treatment of RA.