Problem formulation

Suppose that we have a pair of PPI networks with the graph representations \({\mathcal {G}}_{X} = \left ({{\mathcal {U}},{\mathcal {D}}}\right)\) and \({\mathcal {G}}_{Y} = \left ({{\mathcal {V}},{\mathcal {E}}} \right)\), in which nodes represent proteins in each PPI network (i.e., \(u_{i} \in \mathcal {U}\) or \(v_{j} \in \mathcal {V}\)), and edges (\(d_{ij} \in \mathcal {D}\) or \(e_{ij} \in \mathcal {E}\)) indicate that the corresponding protein u _i (or v _i ) binds with the protein u _j (or v _j ). The edge weights in the PPI networks can indicate the strength or confidence of the interactions between the proteins. Given a pair of nodes across the PPI networks, we assume that the pairwise node similarity score s(u _i ,v _j ), \(u_{i} \in {\mathcal {U}}\) and \(v_{j} \in {\mathcal {V}}\) can be computed, for example, based on the sequence similarity between the proteins. In this study, we utilized BLAST bit scores between proteins as the pairwise node similarity scores. However, other types of similarity measurements (or their combinations) could be also used as the pairwise node similarity score in case such measurements can be easily obtained.

according to the MEA criterion. This MEA approach has been widely used by many multiple sequence alignment algorithms [15–19] and it has been shown to be useful for network alignment [10, 11] and network querying [20] as well.

Motivation and overview of the proposed method

Based on the above problem setting, to construct a confident network alignment, it is crucial to accurately estimate the pairwise node alignment probabilities. To obtain biologically meaningful alignment results, it is necessary that the pairwise node alignment probability is proportional to both the pairwise node similarity (i.e., sequence similarity) and the topological similarity between the subnetwork regions surrounding the nodes in the respective networks. This is based on the observation that orthologous proteins typically have a high level of compositional similarity and often display similar interaction patterns to their neighboring nodes [1, 2]. To accurately estimate the pairwise node alignment probability by effectively integrating these two different types of similarities, we propose to utilize the concept of steady-state network flow (i.e., the amount of ‘water’ that flows through a given channel in the network). Similar concepts have been previously adopted in various engineering applications to find the solutions to similar assignment problems. For example, in digital communication systems, the water-filling algorithm [21] is utilized to compute the optimal allocation of resources. Conceptually, it pours ‘water’ into an OFDM (orthogonal frequency division multiplexing) channel, and the ‘water level’ in the OFDM channel is utilized to find the optimal solution of the transmit power for each subcarrier. In digital image processing, the so-called watershed method [22] is used to find edges or contours of objects in the given image. The watershed method assumes that ‘water’ flows along the image gradient (e.g., intensity differences) and eventually reaches the local minima so that the ‘water level’ in the image provides the solution for the desired image segmentation.

In the proposed method, we measure the steady-state network flow in an integrated network that is obtained by combining the PPI network pair to be aligned. More specifically, edges are inserted between nodes in different networks that have positive pairwise node similarity, and the pairwise node similarity score is assigned as the edge weight (see Fig. 2). Suppose we pour ‘water’ on the integrated network and that the amount of water flow is proportional to the edge weight. If a given pair of nodes in different PPI networks have higher pairwise node similarity and if their neighboring nodes also have higher pairwise node similarity, there would be a larger water flow between the pair of nodes in the long run. However, if the nodes have a similar topological structure (i.e., in terms of the number of interacting nodes in the respective networks) but if their neighboring nodes are not similar, there will be relatively small water flow between the pair of nodes (see Fig. 3). As a result, the water flow between nodes across different PPI networks provides an intuitive way of measuring the overall similarity of the nodes – or functional correspondence between the proteins. As will be shown later, the resulting node correspondence score obtained based on the concept of water flow in the integrated network can serve as an effective building block for constructing an accurate and biologically meaningful network alignment.

Estimating the node correspondence through a Markov random walk model

In order to effectively estimate the node correspondence by integrating both the pairwise node similarity and topological similarity using a Markov random walk model, we first construct the integrated network G=(V,E) by combining \({\mathcal {G}}_{X}\) and \({\mathcal {G}}_{Y}\). Nodes of the integrated network G are the union of the nodes of \({\mathcal {G}}_{X}\) and \({\mathcal {G}}_{Y}\) (i.e., \(V={\left \{ {{\mathcal {U}},{\mathcal {V}}} \right \}}\)), and edges are the union of the edges of \({\mathcal {G}}_{X}\), \({\mathcal {G}}_{Y}\), and additional weighted edges \({\mathcal {F}}\), where \({\mathcal {F}} = \left \{ {s\left ({u_{i},v_{j}} \right)|u_{i} \in {\mathcal {U}},v_{j} \in {\mathcal {V}}} \right \}\) (i.e., \(E={\left \{ {{\mathcal {D}},{\mathcal {E}},{\mathcal {F}}} \right \}}\)). On this integrated network G, we allow the random walker to randomly move from the current node to any of its neighboring nodes at each time step. We define two different types of random moves based on their starting and ending points. First, if the random walker moves from a node in \({\mathcal {U}}\) to a node in \({\mathcal {U}}\) (or from a node in \({\mathcal {V}}\) to a node in \({\mathcal {V}}\)), we define it as an intra-network random move, as the random walk takes place in the same PPI network. Second, if the random walker moves from a node in \({\mathcal {U}}\) to a node in \({\mathcal {V}}\) (or from a node in \({\mathcal {V}}\) to a node in \({\mathcal {U}}\)), we refer to this as a cross-network random move. The intra-network random move mainly aims to capture the topological similarity between the two PPI networks while the cross-network random move aims to incorporate the pairwise node similarity between nodes that originally belong to different PPI networks.

where D _X is a \(\left | {\mathcal {U}} \right | \times \left | {\mathcal {U}} \right |\) dimensional diagonal matrix such that \(D_{X} \left [ {i,i} \right ] = \sum \limits _{\forall j} {A_{X} \left [ {i,j} \right ]} \), and D _Y is a \(\left | {\mathcal {V}} \right | \times \left | {\mathcal {V}} \right |\) dimensional diagonal matrix such that \(D_{Y} \left [i,i\right ] = \sum \limits _{\forall j} {A_{Y} \left [ {i,j}\right ]}\).

where S ^T is a \(\left |{\mathcal {V}}\right |\times \left |{\mathcal {U}}\right |\) dimensional matrix for the pairwise node similarity score, and \( {\textbf {D}_{S^{T}}}\phantom {\dot {i}\!}\) is a \(\left | {\mathcal {U}} \right | \times \left |{\mathcal {U}}\right |\) dimensional diagonal matrix such that \(D_{S^{T}}\left [{i,i}\right ] = \sum \limits _{\forall j} {s^{T}\left [{i,j}\right ]}\). In fact, the transition probability matrices P _X→Y and P _Y→X are normalized pairwise node similarity score matrices in the row-wise and column-wise manner.

Based on the proposed random walk protocol, the random walker transits more frequently between the pair of nodes (u _i ,v _j ) if the node u _i and the node v _j have a higher pairwise node similarity and also if their neighboring nodes also have higher pairwise node similarity (i.e., higher topological similarity). So, as a result, the random walker will spend more time on an edge that connects a pair of nodes (u _i ,v _j ), \(u_{i} \in {\mathcal {U}}\) and \(v_{j} \in {\mathcal {V}}\) as their overall similarity (or node correspondence) increases. Hence, we can effectively estimate the pairwise node alignment probability – which should be proportional to the desired node correspondence – by measuring the steady-state network flow through each edge (u _i ,v _j ), \(u_{i} \in {\mathcal {U}}\) and \(v_{j} \in {\mathcal {V}}\).

To compute the steady-state network flow, we first compute the steady-state probability π(x) of the random walker for every node \(x \in {\mathcal {U}}\cup {\mathcal {V}} \) in the integrated network. This is equivalent to the long-run proportion of time that the random walker spends at a given node x. The steady-state probability distribution is equivalent to the eigenvector of the transition probability matrix P that corresponds to unit eigenvalue. This eigenvector, hence the steady-state probability, can be easily obtained through the power method, as the transition probability matrix P will be generally sparse for real PPI networks [10, 11]. The steady-state probability π(x) can be viewed as the amount of ‘water’ at the node x in the long-run, and since the amount of the water flow is proportional to the edge weight, we can obtain the steady-state network flow along the edge (u _i ,v _j ) as follows (see Fig. 4):

$$ \begin{aligned} c\left({u_{i},v_{j}} \right) &= \pi \left({u_{i}} \right) \cdot \Pr \left[ {v_{j} |u_{i}} \right] + \pi \left({v_{j}} \right) \cdot \Pr \left[ {u_{i} |v_{j}} \right] \\ &= \pi \left({u_{i}} \right) \cdot \frac{{s\left({u_{i},v_{j}} \right)}}{{\sum\limits_{\forall v_{j}} {s\left({u_{i},v_{j}} \right)} }} + \pi \left({v_{j}} \right) \cdot \frac{{s\left({u_{i},v_{j}} \right)}}{{\sum\limits_{\forall u_{i}} {s\left({u_{i},v_{j}} \right)} }}. \\ \end{aligned} $$

(10)

where C is a \(\left |{\mathcal {U}}\right |\times \left |{\mathcal {V}}\right |\) dimensional matrix for the steady-state network flow (i.e., pairwise node correspondence scores), π _X is a \(\left |{\mathcal {U}}\right |\times \left |{\mathcal {U}}\right |\) dimensional diagonal matrix such that π _X [i,i]=π(u _i ), \(u_{i} \in {\mathcal {U}}\), and π _Y is a \(\left | \mathcal {V}\right |\times \left |{\mathcal {V}} \right |\) dimensional diagonal matrix such that π _Y [j,j]=π(v _j ), \(v_{j} \in {\mathcal {V}}\).

After transforming and removing node correspondence scores lower than a specific threshold using Eq. (13), we obtain the final node correspondence scores \(\bar {\mathbf {C}}\), which will be used to construct the network alignment.

Constructing the pairwise network alignment

Finally, to find the optimal solution of Eq. (1) based on the derived pairwise node alignment probability, we construct the maximum weighted bipartite matching (MWBM) between \(\mathcal {G}_{X}\) and \(\mathcal {G}_{Y}\), using an efficient implementation of the MWBM algorithm included in the GAIMC library [23].

Datasets and experimental set-up

We assessed the performance of CUFID-align based on the IsoBase dataset [24], which includes PPI networks of five different species: H. sapiens (human), M. musculus (mouse), D. melanogaster (fly), C. elegans (worm), and S. cerevisiae (yeast). PPI networks in the IsoBase dataset were constructed by integrating five different databases: BioGRID [25], DIP [26], HPRD [27], MINT [28], and IntAct [29]. In IsoBase, the H. sapiens network has 22,369 proteins and 43,757 interactions; the M. musculus network has 24,855 proteins and 452 interactions; the D. melanogaster network has 14,098 proteins and 26,726 interactions; the C. elegans network has 19,756 proteins and 5,853 interactions; and the S. cerevisiae network has 6,659 proteins and 38,109 interactions.

We assessed the quality of the predicted network alignment based on the following metrics: correct nodes (CN), specificity (SPE), gene ontology consistency (GOC) scores, conserved interactions (CI), conserved orthologous interactions (COI), and computation time. Note that CN, SPE, and GOC scores assess the biological significance of the alignment, and CI and COI assess the topological quality of the alignment. If the aligned nodes have the same functional annotation based on the KEGG Orthology (KO) group annotations [30], we considered the node alignment to be correct. CN counts the total number of correctly aligned nodes in a given network alignment. SPE is the relative ratio of the total number of correctly aligned node pairs to the total number of aligned node pairs.

where |c| is the number of proteins having the particular GO term c, and |r o o t(c)| is the total number of proteins under the root GO term of the particular GO term c, where three root GO terms are molecular function (MF, GO:0003674), biological process (BP, GO:0008150), and cellular component (CC, GO:0005575). Note that if at least one protein in the aligned protein pair does not have a functional annotation such as KO group annotations or GO terms, the aligned protein pair was removed before computing the performance metrics CN, SPE, and GOC scores.

We compared the performance of CUFID-align against a number of state-of-the-art alignment methods: IsoRank [6], SMETANA [10], SMETANA-CSRW [11], PINALOG [13], and HubAlign [14]. Additionally, to verify the effectiveness of the network-based approach over the conventional approach that uses sequence similarity alone, we compared the various network-based methods and with a method that finds the best mapping between networks solely based on the sequence similarity between the proteins. More specifically, given a network pair, we tried to predict the network alignment by using maximum weighted bipartite matching based on BLAST bit scores. Since both SMETANA and SMETANA-CSRW yield many-to-many mappings by default, we used the parameter n _max =1 to obtain one-to-one mappings. Other than this, the default parameters were used in our experiments (i.e., α=0.9 and β=0.8). For HubAlign, we used the default parameters (i.e., λ=0.1, d=10, and α=0.7). For IsoRank, we set the parameter α=0.6 as recommend in the original paper. For CUFID-align, we set the parameter α=0.9 and β=90 percentile of the transformed correspondence score. We performed all experiments on a desktop computer equipped with a 3.2 GHz Intel i5 quad-core processor and 8 GB memory.

Performance assessment based on the IsoBase dataset

When it comes to the specificity of the alignment results, random walk based methods (CUFID-align, SMETANA-CSRW, and SMETANA) achieve relatively higher SPE compared to PINALOG and HubAlign. SPE of HubAlign appears to be more sensitive than the other methods with respect to the average degrees of the PPI networks. CUFID-align, SMETANA-CSRW, and SMETANA achieve similar SPE, often higher than those of PINALOG and HubAlign. This means that CUFID-align can in general more accurately align protein pairs that have the same functional annotations compared to PINALOG and HubAlign.

Since proteins can have multiple functions, we further evaluated the functional consistency of the alignment results based on the GOC scores, where higher GOC scores indicate that the obtained alignments are functionally more coherent. As we can see in Fig. 5, CUFID-align achieves higher GOC scores than the other compared algorithms in all test cases. Again, if the network pairs have higher average degrees, PINALOG and HubAlign show comparable GOC scores. However, probably due to the restricted search space of the seed-and-extend approach, GOC scores of PINALOG and HubAlign tend to be smaller than the other methods when the average degree of one of the PPI networks is relatively smaller than that of the other. In comparison, CUFID-align is more robust to the change of topological properties such as the average degrees of the PPI networks to be aligned.

The above results show that CUFID-align can accurately predict matching proteins in different species that have similar functionalities, according to the functional annotations of proteins that are currently available. The results also imply that the proposed algorithm may provide a useful tool for predicting the functions of unknown proteins in less studied species through network alignment with species that have been better studied.

We also compared the network-based approaches with the sequence-similarity-based approach. As we can see in Table 1 and Fig. 5, a simple sequence-similarity-based approach can construct network alignments with high functional coherence, and that the node similarity score may provide useful guidelines for identifying orthologous proteins. However, this results should be taken with a grain of salt, since they are likely due to the fact that the current functional annotations of proteins are often based on sequence similarity between proteins. As shown in Tables 2 and 3, BLAST-MWBM – which uses BLAST bit score and MWBM without using any network information – can identify a much smaller number of CIs and COIs compared to the network-based methods. These results imply that strong dependence on sequence similarity for constructing a network alignment has the potential risk of getting biased results that may fail to capture important protein interactions that are conserved across different species, which may be critical in deciphering the underlying cellular mechanisms that involve those interactions. In contrast, network-based methods, including CUFID-align, that incorporate topological information for constructing network alignments can make accurate and balanced predictions that identify both orthologous proteins as well as conserved interactions. Our results clearly show the importance of effective integration of node similarity and topological similarity for effective comparative analysis of PPI networks.

Extension of CUFID-align for the alignment of multiple networks

In this work, we have focused on the steady-state network flow approach and its application to the pairwise network alignment problem. However, the problem of multiple network alignment has been gaining wide interest in the research community and its practical importance has been increasing as the number of available PPI networks for different species continue to increase. Although it is beyond the scope of the current paper, we expect the extension of CUFID-align for multiple network alignment will be relatively straightforward. First of all, to this aim, we can modify the transition probability matrix in Eq. (9) by concatenating the normalized adjacency matrices and node similarity score matrices for the multiple PPI networks to be aligned. Following the construction of this extended transition probability matrix, the steps for computing the node correspondence scores – shown in Eqs. (11) and (13) – can be modified by constructing diagonal matrices and inserting corresponding the matrices into the diagonal terms. The extended version of CUFID-align for multiple PPI network alignment is expected to have distinctive advantages over other existing multiple PPI network alignment algorithms. First, it may be able to estimate the ‘global’ node correspondence scores more accurately. Currently, most multiple PPI network alignment algorithms estimate the node correspondence scores for every PPI network pair in the interest of computational complexity. The estimated pairwise node correspondence scores are later updated based on additional transformations to make them more suitable for multiple network alignment. However, considering that the ultimate goal is in constructing the alignment of multiple networks, it would be preferable to estimate the node correspondence scores (or equivalently, node alignment probabilities) Pr[u _i ∼v _j |G] considering all networks, rather than just estimating \(\Pr \left [ {u_{i} \sim v_{j} |{\mathcal {G}}_{X},{\mathcal {G}}_{Y}} \right ]\) based on the given network pair, where \(u_{i} \in {\mathcal {G}}_{X} \), \(v_{j} \in {\mathcal {G}}_{Y} \), and G is the set of all PPI networks including \({\mathcal {G}}_{X}\) and \({\mathcal {G}}_{Y}\). Since the aforementioned extension of CUFID-align estimates the node correspondence scores based on an integrated network that combines all networks in G, it has the potential to accurately compute the posterior node-to-node alignment probability given all the networks. Computation of such ‘global’ node correspondence score may lead to improved multiple network alignment results. Second, the extended version of CUFID-align will still be computationally very efficient, as most steps in CUFID-align only require simple matrix operations even if extended to multiple networks. Finally, the extended approach will require relatively low computational resources (especially, in terms of memory). For example, suppose that there are N PPI networks, where the number of nodes in the i-th network G _i is V _i . To align the N PPI networks, IsoRankN will need the pairwise node correspondence scores for each of the \(\binom {N}{2}\) network pairs, where for each pair, the algorithm will need to construct a |V _i ·V _j |×|V _i ·V _j | dimensional matrix. However, CUFID-align can compute the global node correspondence scores by constructing a single \(\left | {\sum \limits _{i = 1}^{N} {V_{i}} } \right | \times \left | {\sum \limits _{i = 1}^{N} {V_{i}}}\right |\) dimensional matrix. We are currently working on extending CUFID-align for multiple network alignment.