Inferring pleiotropy by network analysis: linked diseases in the human PPI network

Background Earlier, we identified proteins connecting different disease proteins in the human protein-protein interaction network and quantified their mediator role. An analysis of the networks of these mediators shows that proteins connecting heart disease and diabetes largely overlap with the ones connecting heart disease and obesity. Results We quantified their overlap, and based on the identified topological patterns, we inferred the structural disease-relatedness of several proteins. Literature data provide a functional look of them, well supporting our findings. For example, the inferred structurally important role of the PDZ domain-containing protein GIPC1 in diabetes is supported despite the lack of this information in the Online Mendelian Inheritance in Man database. Several key mediator proteins identified here clearly has pleiotropic effects, supported by ample evidence for their general but always of only secondary importance. Conclusions We suggest that studying central nodes in mediator networks may contribute to better understanding and quantifying pleiotropy. Network analysis provides potentially useful tools here, as well as helps in improving databases.


Background
The systems perspective on complex biological systems emphasizes that individual genes act in genetic networks and individual proteins play their roles in protein-protein interaction (PPI) networks [1]. There is increasing interest in these networks, as their analysis helps to understand the relationship between the components (i. e. genes, proteins) and how these are positioned in the whole system. Well-connected hubs seem to be of high functional importance [2,3]. Consequently, studies on diseases based on PPI networks had the starting point by analysing the centrality of disease proteins. Genes associated with a particular phenotype or function are not randomly positioned in the PPI network, but tend to exhibit high connectivity; they may cluster together and can occur in central network locations [4,5].
Beyond focusing on the number of neighbours of graph nodes (their degree), wider neighbourhoods, indirect effects and larger subsets of nodes can also be analyzed by the rich arsenal of network analytical tools. This non-local information may help, for example, to quantify the structural relationships between different sets of proteins. In an earlier paper [6], we have determined proteins that mediate indirect effects between sets of proteins causing five diseases in the human PPI network. Their mediator role was quantified and they were ranked according to structural importance. Their functional role may be of high interest, as proteins involved in certain pairs of diseases have no direct interactions among them [6]. These findings motivated an appealing problem: "which proteins connect diseases in the human PPI network?".
To be connected to diverse regions of the PPI network may lend a functionally pleiotropic character to a protein in a classical, genetic sense: it has been demonstrated that high connectivity correlates well with pleiotropic effects [7,8]. The most central mediators are especially important in connecting apparently distant nodes in the human PPI network. Specific network positions may render strange but characteristic behaviour (expression pattern) to different proteins [9,10]. Instead of being exceptional, these epistatic effects may be of primary importance in physiology [11] and in better understanding animal development and adaptation.
In this paper, (1) we compare two interaction networks of mediators (mediating indirect effects between heart disease and obesity, and between heart disease and diabetes), (2) we analyse the structure of these two networks and their aggregated total network, (3) we study the overlap between the two mediator networks, and (4) we infer biological functions for some proteins and provide supporting literature data. All in all, we illustrate that network analysis is an excellent tool for identifying pleiotropy and epistasis from complex networks extracted from multiple databases.

Network analysis
We obtained 9 proteins involved in heart diseases (H), as well as 44 and 20 involved in diabetes (D) and obesity (O), respectively. The HD network contains N = 2142 nodes and L = 3537 links, while the HO network contains N = 1746 nodes and L = 2567 links and the total network contains N = 2221 nodes and L = 3686 links. Figure 1 provides a schematic illustration for how the networks had been constructed (see Methods). Figure 2 shows the relationships between mediator proteins in the HD (Figure 2a The distributions of individual structural indices are very similar for all of the three analyzed networks. Additional file 1 shows all values of the six network indices for all nodes in the three networks. Figure 5 shows these distributions only for the total network. We can observe that almost all indices follow a strongly leftskewed distribution where only a few nodes are extremely important. While degree (D), topological importance (TI) and betweenness centrality (BC) have really only one or a few hubs, topological overlap (TO) indicates several key nodes. Closeness centrality (CC) has a unimodal, normal-like distribution. For each network, there seem to be strong and positive rank correlation between all centrality indices but not for the overlap indices (TO 3 0.01 and TO 3 0.005 ). TO indices correlate positively and weakly with other centrality indices whereas they correlate negatively and weakly with CC (see Table 1). D best correlates with TI 3 . The TO measure offers different, complementary information than the centrality indices. Table 2 summarizes the results of the randomization  test (note that only the means are shown in the table, for simplicity). The observed rank correlation coefficients are all significantly lower than those for the random networks (with 95% confidence interval). This suggests that there are stronger rank correlations between different centrality indices in the random networks, in comparison to the results obtained from the HD, HO and total networks. One possible explanation for this discrepancy is that, beyond the mathematical properties, real networks are structured also by biological constraints. Thus, different centrality indices can capture different aspects of network topology, therefore correlation between different indices are weaker for real networks. This provides more support on using various network indices to capture different topological properties embedded in real networks.

Biological results
We now examine more closely the rank order of the top nodes in each network. The degree ranks for the three networks are almost identical (see Tables 3, 4 and 5). The most central nodes are P62993 (Growth factor receptor-bound protein 2), P63104 (14-3-3 protein zeta/ delta) and P06241 (Tyrosine-protein kinase Fyn). The 9 shared proteins rank in the same order in HD and HO and there is no change in rank order also in the total network. In the HO network, the 12 mediators lead the ranking, and then come their neighbours. However, in the HD rank (and also in the total network), there is one non-mediator protein in the top 26 of the rank (among the 25 HD mediators); this is P00533 (Epidermal growth factor receptor) in the 23 rd position.
The betweenness ranks correspond quite well to the degree ranks with some exceptions. For example in the HD network, P06241 (Tyrosine-protein kinase Fyn) is three positions lower in betweenness ranking when compared to its degree rank position. In the HD network, instead of one, now five non-mediators are mixed with HD mediators in the top of the list, while some HD mediators such as Q99616 (C-C motif chemokine 13) lose their high degree-based rank completely. In contrast, the degree rank order seems to be consistent with its betweennes counterpart for HO and total networks.
Despite the large overlap between the HD and HO networks, the rank positions of HD and HO mediator proteins are quite different in the two networks. For example, both P17302 (Gap junction alpha-1 protein) and P43405 (Tyrosine-protein kinase SYK) rank high in the HD network but not in the HO network. As it is shown on Figure 2b, O14908 (PDZ domain-containing protein GIPC1) is the only protein among the three exclusive HO mediators that is part of the interaction network of HO mediators.
Additional File 2 shows the extracted GO terms of proteins ranked by different structural indices for the HD network, the HO network and the total network. For example, by considering the top 30 proteins ranked by degree in the HD network, we found that half of them are related to the processes 'intracellular signaling cascade' (GO:0007242) and 'protein amino acid phosphorylation' (GO:0006468), meanwhile in the HO Figure 1 The HD and HO mediator networks and their subnetworks. Red, blue and yellow proteins are involved in three diseases (H: heart diseases, D: diabetes, O: obesity). Pink proteins mediate indirect effects between the red and the blue ones, while orange proteins mediate between the red and the yellow ones. Black proteins mediate between both pairs. White proteins are the non-mediator neighbours of the mediator proteins. We analyzed five networks: the HD mediator network (pink and black nodes with their white neighbours), the HO mediator network (orange and black nodes with their white neighbours), the total mediator network (pink, orange and black nodes with their white neighbours), the subnetwork of interactions among HD mediators (pink and black nodes) and the subnetwork of interactions among HO mediators (orange and black nodes). network 16 of them are located in 'plasma membrane' (GO:0005886) and 13 of them are related to process 'cell surface receptor linked signal transduction' (GO:0007166).
The p-values of proteins quantify their average fit to the studied GO-terms (i.e. to what extent they can be characterized by certain functionality). By comparing those p-values to centrality and overlap indices used in this study, we can conclude that the performance of different indices vary strongly. In the total network, only the TO 3 0.01 index correlates significantly with biological function ( Table 6). Note that the performance of  depends on the t threshold used. Proteins in unique positions are, thus, typically involved in the above-mentioned key functions. The other relatively well-performing index is CC, whereas D and BC correlate with function only once each. TI 3 and TO 3 0.005 do not correlate with functions defined by GO terms. Furthermore, functional roles are best predictable by these structural indices in the HD network and less so for the HO network.

Discussion
Based on its centrality ranks, P63104 (14-3-3 protein zeta/delta) corresponding to gene YWHAZ seems to be the second most important protein in these mediator processes. This is in concert with the literature, stating that P63104 is a chaperon [12] and is richly connected to several kinds of other molecules with mostly weak links [13]. Specifically, it is involved in cell growth and carcinogenesis [14], breast cancer reoccurrence after chemotherapy resistance [15], luteal sensitivity to PGF [16] and, finally, it is part of antiapoptotic (P13K/AKT) and cell proliferation (ERK/MAPK) pathways [17]. Its connecting position has been demonstrated by network analysis, showing its involvement in several HSNs (highscoring subnetworks [18]). Ogihara et al. [19] suggested that the association with 14-3-3 protein may play a role in the regulation of insulin sensitivity by interrupting the association between the insulin receptor and IRS1. It means that P63104 probably mediate HD and HO through the regulation process of insulin (as insulin is a crucial hormone in human metabolic system). Typically it is not directly responsible for diseases (not assigned to any disease in the OMIM database) but very frequently mentioned as a candidate protein in the background, requiring further investigation [14].
The most important protein, P62993 (Growth factor receptor-bound protein 2) corresponding to gene GRB2 leads in all of the six structural importance ranks. It appears in the mammalian Grb2-Ras signaling pathway with SH2/SH3 domain interactions and several functions in embryogenesis and cancer [20]. Zhang et al. [21] also found that GRB2 is essential for cardiac hypertrophy and fibrosis in response to pressure overload and that different signaling pathways downstream of GRB2 regulate fibrosis, fetal gene induction, and cardiomyocyte growth. Yet, in the subgraph of the HO mediators, P62993 does not seem to occupy a central position but its phenotypic traits are likely to be affected through the links to non-mediators instead of other HO mediators. This kind of structural arrangement is advantageous for information integration, while a strongly connected mediator subnetwork implies functional redundancy.
Among the three exclusive HO (non-HD) mediators, O14908 (PDZ domain-containing protein GIPC1) corresponding to the gene GIPC1 appears in the HO mediator subgraph, while the other two are isolated (Q14232 -Translation initiation factor eIF-2B subunit alpha corresponding to gene EIF2B1; Q5JY77 -G-protein coupled receptor-associated sorting protein 1 corresponding to gene GPRASP1). This may suggest also that O14908 is an HD mediator. Its connection to heart disease is clear but its interaction with diabetes-related proteins is not a b documented in the OMIM databases (also not for the other two proteins). However, this inferred function is well supported by Klammt et al. [22] reporting on the role of O14908 in diabetes. A possible outcome of network analysis is to suggest potential updates in the databases. The only protein that ranks higher than HD mediator proteins in the degree-based centrality rank of the HD network is P00533 (Epidermal growth factor receptor), corresponding to the gene EGFR. We could speculate that this protein might also mediate between H and D proteins. In the total PPI network, it is linked to two D proteins (Q9UQF2 -JNK-interacting protein 1; Q9UQQ2 -Signal transduction protein Lnk) but not to H protein. EGFR and its ligands are cell signaling molecules involved in a wide range of cellular functions, including cell proliferation, differentiation, motility, and tissue development [23]. Research on EGFR's pathogenesis have been focused on lung cancer [24] and have not discovered its link to heart diseases. However, Iwamoto and his colleagues observed the role of ErbB signaling in heart functions [25]. Also, it has been shown to be a central protein according to other sophisticated network analysis techniques [26], dominating the clique composition of certain pathways. Based on our static, structural inference, it is not easy to decide whether a protein is "strongly linked to a disease" or it is a "disease protein". The definitions are very poor here. Is P00533 a H protein (causing heart diseases) or HD protein (mediating between H and D proteins)? The solution is to use inference for generating new hypotheses, improving databases and designing experiments, instead of regarding the inferred findings as results.

Conclusions
Our study focused on only a few diseases but the approach and the methods used can be generalized. It may be interesting to extend this research to other diseases and to study the pleiotropic effects of mediators linking other disease pairs. The mediator proteins analyzed in this study typically have pleiotropic effects. They connect several pathways and influence several phenotypic traits. The reason why their inferred structural roles miss from the OMIM database is exactly that they act in a non-Mendelian way. They are typically not the singular elements of important pathways but weak connectors among several pathways of high importance. This way, their effects can be fundamental. Their understanding needs a multi-locus, systems-based, network view. As individual pathways are linked to networks, our non-Mendelian knowledge on linkage, epistasis and pleiotropy becomes larger. If network analysis makes these epistatic and pleiotropic effects quantifiable and predictable, we are getting closer to better understand delegated complexity [27]. From an application perspective, it would be interesting to see whether a healthy (intact and well-connected) network of mediators could contribute to healthy phenotypes or, in contrary, disconnecting the mediator network could be used to isolate diseases and reduce side-effects of drugs.

Data
We have analyzed human protein-protein interaction network (PPI) data extracted from the I2D database. I2D (Interologous Interaction Database) is an on-line database of known and predicted mammalian and eukaryotic protein-protein interactions [28]. It is one of the most comprehensive sources of known and predicted eukaryotic PPIs. We carefully considered the completeness of the PPI network by investigating various human PPI databases. In their database, the Authors have collected data from almost all of the well-known human protein interaction databases including HRPD http://www.hprd.org/, BIND http://bind.ca/, MINT http://mint.bio.uniroma2.it/mint/ and Intact http://www.ebi.ac.uk/intact/, among others. Those databases are built by arrange of methods, some are experimental ones, some are predicted ones, and some are curated from the literature. By using the I2D database, we could thus construct the network integrated from multiple data sources. We investigated other databases not included in the I2D database, particularly the STRING database http://string.embl.de/ and we found that almost all high-scoring interactions in STRING were covered in our data set. Combining data from various sources is supposed to be more comprehensive for analyzing the PPI network than studying each data source separately. To obtain a more reliable set of protein interactions, we excluded all the interactions obtained by homology methods: only experimentally verified ones were included in our analysis. For the disease phenotypes, the clinical Online Mendelian Inheritance in Man database (OMIM, [29]) was investigated. We have checked whether we need to update our database used in Nguyen and Jordán [6] and found that we can use the same data set as the number of updates is negligible.

Analysis
From the human PPI network data, we constructed: (1) a network of proteins mediating indirect effect between heart disease (H) and diabetes (D) proteins (i.e. HD mediators) and their direct neighbours (i.e. HD network); (2) a network of proteins mediating indirect effect between heart disease (H) and obesity (O) proteins (i.e. HO mediators) and their direct neighbours (i. e. HO network); and (3) an aggregated network of the two previous networks (i.e. total network). We considered only two-step mediator proteins, directly connected to two proteins related to different diseases and being otherwise unconnected (so, we do not consider chains of mediators). We have also studied the subnetworks of (1) and (2) without non-mediator neighbours. See Figure  1 for schematically illustrating the relationships between these five networks. Figure 2 shows the subnetworks without non-mediator neighbours (Figure 2a for HD and Figure 2b for HO). Figure 3a shows the HD and Figure 3b shows the HO network. The total network is shown in Figure 4.
Earlier we have determined the identity of these HD and HO mediators and quantified the strength of their mediator effect [6]. Here, we focus on the networks of

0.775
The Spearman rank correlation coefficients between each pair of centrality indices for the HD and HO networks as well as the total network.

0.945
The mean Spearman rank correlation coefficients between each pair of centrality indices obtained from 1000 random networks of the same size as the HD, HD and the total network. mediators. Links in these networks are undirected (if protein i is linked to protein j, then j is also linked to i) and unweighted (we have no data for the intensity or strength of the interactions). We have characterized each network by some simple network statistics. (i) The simplest index that provides the most local information about node i is its degree (D i ). This is the number of other nodes connected directly to node i. We have calculated the normalized degree: where N is the number of nodes in the network.
(ii) A measure of positional importance quantifies how frequently a node i is on the shortest path between every pair of nodes j and k. This index is called "betweenness centrality" (BC i ) and it is used routinely in network analysis [30]. The normalized betweenness centrality index for a node i (nBC i ) is: where i ≠ j and k; g jk is the number of equally shortest paths between nodes j and k, and g jk (i) is the number of these shortest paths to which node i is incident (g jk may equal one). The denominator is twice the number of pairs of nodes without node i. This index thus measures how central a node is in the sense of being incident to many shortest paths in the network.
(iii) "Closeness centrality" (CC i ) is a measure quantifying how short are the minimal paths from a given node i to all others [30]. The normalized index for a node i (nCC i ) is: The rank of the most central 30 nodes in the HD network, based on the six importance indices analyzed.
where i≠j, and d ij is the length of the shortest path between nodes i and j in the network. This index thus measures how close a node is to others. The larger nCC i is for node i, the more directly its deletion will affect the majority of other nodes.
(iv) Topological importance can also be quantified by general matrix algebra. In an undirected network, we define a n,ij as the effect of j on i when i can be reached from j in n steps. The simplest way of calculating a n,ij is when n = 1 (i.e. the effect of j on i in 1 step): where D i is the degree of node i (i.e. the number of its direct neighbours). We assume that indirect effects are multiplicative and additive. For instance, we wish to determine the effect of j on i in 2 steps, and there are two such 2-step pathways from j to i: one is through k and the other is through h. The effects of j on i through k is defined as the product of two direct effects (i.e. a 1,kj ×a 1,ik ), therefore the term multiplicative. Similarly, the effect of j on i through h equals to a 1,hj,1 ×a 1,ih . To determine the 2-step effect of j on i (a 2,ij ), we simply sum up those two individual 2-step effects: and therefore the term additive. When the effect of step n is considered, we define the effect received by node i from all other nodes in the same network as: which is equal to 1 (i.e. each node is affected by the same unit effect.). Furthermore, we define the n-step effect originated from node i as: The rank of the most central 30 nodes in the HO network, based on the six importance indices analyzed. The rank of the most central 30 nodes in the total network, based on the six importance indices analyzed.
which may vary among different nodes (i.e. effects originated from different nodes may be different). Here, we define the topological importance of node i when effects "up to" n step are considered as: which is simply the sum of effects originated from node i up to n steps (one plus two plus three...up to n) averaged over by the maximum number of steps considered (i.e. n). This TI n index measures the positional importance of a node by considering how effects originated from such a given node can spread through the whole network to reach all nodes after a pre-defined n step length [31]. Calculations were performed by the CosBiLAB Graph software [32].
(v) Basically every node in a network is connected to each other, but it still matters how strongly they are connected (whether two nodes are neighbors in the network, second neighbors or more distant ones). Thus, it is of interest to study the indirect neighborhood of particular nodes, considering more than only the neighbors but less than the whole network. For a given step length n and a given network, there is an interaction matrix presenting the relative strengths of interactions between each pair of nodes i and j. We note that interaction strength is used here in a totally structural sense, with no dynamical component. If n exceeds 2 or 3, and the network is not very large, then there is non-zero interaction strength between each pair of nodes (everything is connected to everything else). Thus, an effect threshold (t) can be set, determining the "effective range" of the interaction structure of a given graph node i, and nodes within this effective range are defined as strong interactors of i (i.e. effects received from i being greater than t) whereas nodes outside this range are defined as i's weak interactors (effects received from i is less than t). Since the sets of strong interactors of two or more nodes may overlap, it is possible to quantify this overlap (the number of shared strong interactors) in order to measure the positional uniqueness of individual graph nodes. The topological overlap between nodes i and j up to n steps (TO n t, ij ) is the number of strong interactors appearing in both i's and j's effective ranges determined by the threshold t. The sum of all TO-values between node i and others provides the summed topological overlap of node i: TO n t,ij (i = j).
For simplicity of representation, we drop the subscript i for all indices. A more detailed description of this index can be found in [33]. Calculations were performed by the CosBiLAB Graph software [32]. Two thresholds have been used, t 1 = 0.01 and t 2 = 0.005.
Each of the six above mentioned structural indices were determined for every node in the networks. The 30 most central ones are presented for the HD network (Table 3), the HO network (Table 4) and the total network (Table 5). Additional File 2 presents all index values for all nodes in these networks.
Since different network indices provide different rankings, it is a question of how similar these rankings are. Similarity refers to robust importance ranks (irrespective to the index), while dissimilarity refers to the complementary information content of the different indices. For statistical analysis, we calculated the Spearman rank correlation coefficient for each pair of the indices in the three major networks (Table 1).
In order to better understand the ranking of nodal indices, we determined the distribution of each structural index for each network. We present these distributions for the total network in Figure 5. To test the significance of the observed rank correlation coefficients, we have constructed random networks. For each of our observed networks (i.e. HD, HO, total), we calculated the probability of two nodes being linked together: We have constructed 1000 random networks with fixed N and a p link probability. For each random network, we calculated the same centrality indices and determined the Spearman rank correlation coefficient for each pair of centrality indices. Since we have 1000 random networks, for each pair of centrality indices we thus have 1000 Spearman rank correlation coefficients. From their distribution, we determined the mean and the 95% confidence intervals. Results are summarized in Table 2.
For the top 30 nodes ranked by a particular index in a particular network, we quantified their biological function by calculating the p-values of GO terms [34]. Specifically, we determined the ratio of the top 30 nodes that can be characterized by a certain GO term and computed the associated p-values (Table 6). Bold numbers mean p < 0.05.