Maximum expected accuracy (MEA) alignment of biological networks

Let us assume that we have a set of N PPI networks $G=\left\{{\mathcal{G}}_1,{\mathcal{G}}_2,\dots ,{\mathcal{G}}_N\right\}$. Each network ${\mathcal{G}}_n=\left({\mathcal{V}}_n,{\mathcal{E}}_n\right)$ has a set of nodes ${\mathcal{V}}_n=\left\{v_1,v_2,\dots \right\}$ and edges ${\mathcal{E}}_n=\left\{e_{i,j}\right\}$, where e _i,j represents the interaction between nodes v _i and v _j in the network ${\mathcal{G}}_n$. For each pair of PPI networks ${\mathcal{G}}_{\mathcal{U}}=\left(\mathcal{U},\mathcal{D}\right)$ and ${\mathcal{G}}_{\mathcal{V}}=\left(\mathcal{V},\mathcal{E}\right)$, we denote the pairwise node similarity score for a node pair (u _i , v _j ), where $u_i\in \mathcal{U}$ and $v_j\in \mathcal{V}$, as s(u _i , v _j ). In this study, we use the BLAST bit score between proteins as their node similarity score, but other types of similarity scores based on structural or functional similarity can be also utilized if available.

A similar MEA approach [18] has been formerly adopted by a number of multiple sequence alignment algorithms, including ProbCons [17], ProbAlign [19], and PicXAA [20–22]. The MEA framework has been shown to be very effective in constructing accurate alignment of multiple biological sequences, making it one of the most popular approaches for sequence alignment. Recently, the MEA approach has been also applied to comparative network analysis, where RESQUE [4] performs MEA-based network querying and SMETANA [12] performs MEA-based multiple network alignment.

Comparing and aligning networks based on context-sensitive random walk

In order to find the alignment that maximizes the expected accuracy defined in (2), we first need an accurate method for estimating the posterior node alignment probability P (u _i ~ v _j | G). For this purpose, we adopt a context-sensitive random walk (CSRW) model, motivated by the pair hidden Markov model (pair-HMM) that has been widely used in sequence alignment [23]. The pair-HMM provides a simple, yet very effective, mathematical framework for estimating the alignment probability between symbols in different biological sequences. Unlike the traditional HMM, which generates a single symbol sequence, the pair-HMM generates a pair of aligned symbol sequences. Pair-HMM makes transitions between three different internal states M, I _X , and I _Y , where the M state emits an aligned pair of symbols, one symbol in sequence X and the other in sequence Y, while I _X and I _Y emit an unaligned symbol in sequence X and sequence Y, respectively. Given two biological sequences, the pair-HMM can be used to estimate the probability whether a given symbol pair was jointly emitted at state M, hence should be aligned to each other. This probability can be computed using the forward and backward algorithms and the resulting alignment probability provides us with a measure of confidence about the (biological) relevance between the given symbols (i.e., nucleotides, amino acids).

One of the most important features of pair-HMM is that it properly recognizes that conserved sequence patterns and motifs in different species may contain inserted and/or deleted symbols (often referred to as "indels") and therefore it specifically tries to model these indels. In a similar manner, a mathematical model that can recognize node insertions and deletions in different biological networks that contain conserved subnetwork regions and network motifs may be useful for obtaining a reliable posterior node-to-node alignment probability. Recently, random walk models have been shown to be effective for estimating the node correspondence in different networks [7, 12, 15] in a way that seam-lessly integrates both node similarity and topological similarity. However, the random walk models that were used in previous network alignment algorithms did not explicitly consider indels.

In this work, we adopt a novel context-sensitive random walk model that has been recently proposed to improve on existing models by taking such indels into account [24]. In a way that is conceptually similar to the pair-HMM, the CSRW has three different internal states M, $I_{\mathcal{U}}$, and $I_{\mathcal{V}}$, each of which corresponds to a different mode of random walk. At the M state, the random walker simultaneously moves on both networks to enter a pair of "matching" nodes. On the other hand, at the $I_{\mathcal{U}}$ state, the random walker only moves on network ${\mathcal{G}}_{\mathcal{U}}$ to enter a potentially "inserted" node in ${\mathcal{G}}_{\mathcal{U}}$ that may not have a corresponding node in the network ${\mathcal{G}}_{\mathcal{V}}$. Similarly, at the $I_{\mathcal{V}}$ state, the random walker only moves on ${\mathcal{G}}_{\mathcal{V}}$ to enter a potentially inserted node in ${\mathcal{G}}_{\mathcal{V}}$. Transitions between states take place in a context-sensitive manner, where the random walker examines the neighboring nodes to determine the mode of random walk. For example, if there are node pairs with significant node similarity (i.e., potential orthologous nodes) in the immediate neighborhood, the CSRW switches to the M state to make a simultaneous move on both networks and randomly enter one of these node pairs. Otherwise, the CSRW switches to either $I_{\mathcal{U}}$ or $I_{\mathcal{V}}$ and performs an individual random walk only on one of the networks. Based on this random walk model, we compute the long-run proportion of time that a given pair of nodes will be simultaneously visited (i.e., at the M state), which can be used to compute a probabilistic correspondence score between these two nodes, as we will describe in the following section.

Estimation of node correspondence scores

Let the current position of the random walker in the product graph be (u _c , v _c ), where $u_c\in \mathcal{U}$ and $v_c\in \mathcal{V}$. In each time step, the random walker examines the set of similar neighboring nodes $\mathcal{N}\left(u_c,v_c\right)=\left\{\left(u_i,v_j\right)|u_i\in \mathcal{N}\left(u_c\right),v_j\in \mathcal{N}\left(v_c\right),\left(u_i,v_j\right)\in \mathcal{M}\right\}$ to determine its mode of random walk (corresponding to one of the three possible internal states), where $\mathcal{N}\left(u_c\right)$ is the set of neighbors of the node u _c in the network ${\mathcal{G}}_{\mathcal{U}}$ and $\mathcal{N}\left(v_c\right)$ is the set of neighbors of the node v _c in the network ${\mathcal{G}}_{\mathcal{V}}$. If there are similar node pairs among the neighboring node pairs, hence $\mathcal{N}\left(u_c,v_c\right)$ is not empty, the random walker switches its internal state to the M state and performs a simultaneous walk on both networks, moving from (u _c , v _c ) to one of the nodes

for $v_j\in \mathcal{N}\left(v_c\right)$.

Based on the transition probabilities given by (6), (7a), and (7b), we can construct the transition probability matrix P for the random walk on the two networks ${\mathcal{G}}_{\mathcal{U}}$ and ${\mathcal{G}}_{\mathcal{V}}$. Given P, we can estimate the longrun proportion of time that the random walker spends in each pair of nodes (u _i , v _j ) by computing the steady state distribution π. In practice, since real PPI networks typically have a relatively small number of interactions (therefore only few edges for most nodes), the resulting transition probability matrix for the CSRW is sparse, which makes it relatively straightforward to compute the steady state distribution using the power method.

In order to increase the computational efficiency of the proposed network alignment method, instead of using the original transition probability matrix P, we use a reduced matrix $\tilde{P}$. The reduced matrix $\tilde{P}$ is obtained by removing the rows and columns in P that correspond to node pairs in $\mathcal{I}$ while keeping only the rows and columns that correspond to node pairs in $\mathcal{M}$. After the reduction, $\tilde{P}$ is re-normalized to make it a legitimate stochastic matrix. In practice, since the CSRW is designed to spend more time at node pairs with higher similarity, the random walker spends a relatively small amount of time at node-pairs that belong to the set $\mathcal{I}$, and using the reduced matrix $\tilde{P}$ instead of P only minimally affects the estimated long-run proportion of time spent at $\left(u_i,v_j\right)\in \mathcal{M}$. As a result, the difference in terms of network alignment performance that results from replacing the original matrix P by this reduced matrix $\tilde{P}$ appears to be small as shown in the supplementary material (see Section S1).

The above formulation, obtained by allowing the CSRW to restart the random walk at a new position, is especially useful when comparing real PPI networks, which are often incomplete and contain many isolated nodes. Simulation results show that the incorporation of the restart scheme can make our CSRW-based alignment method more robust, especially when the available topological data are either unreliable or insufficient for detecting the similarities between networks (see Section S2).

In order to determine the restart probability λ, we first analyze the structure of the reduced product graph of ${\mathcal{G}}_{\mathcal{U}}$ and ${\mathcal{G}}_{\mathcal{V}}$ that contains only similar node pairs included in $\mathcal{M}$. Intuitively, it is desirable to increase the restart probability λ if the networks are disconnected and decrease the probability if the networks are well connected. For example, if all the nodes in the reduced product graph are completely disconnected, it is desirable to restart the random walker at every step. Additionally, when we consider the following two cases - (i) most nodes in the product graph are connected and there are only a few disconnected nodes; (ii) the product graph is equally divided into N connected subnetworks of identical size - it would be desirable to assign a higher λ to the latter case. Based on these intuitions, we set the restart probability λ as the ratio of the total number of nodes in the top K% smallest subnetworks to the total number of nodes in the reduced product graph. In this work, we used K = 99% to determine the restart probability λ.

Constructing the multiple network alignment

Once we have computed the node correspondence scores in (8) for every pair of networks in G, we take a greedy approach as in [12] to construct the multiple network alignment. The overall alignment process is as follows. First, in order to improve the reliability of the node correspondence scores, we selectively apply the probabilistic consistent transformation (PCT) defined in [12]. If λ is larger than a predefined threshold λ _t , we do not apply PCT to the node correspondence scores. A large λ implies that the product graph is ill connected (e.g., containing a large number of isolated nodes), in which case applying the PCT would not be helpful and may in fact make the scores less reliable. This is because the PCT in [12] was developed based on the assumption that the product graphs for all network pairs are relatively well connected. After the potential score refinement step through PCT, we begin with an empty alignment and greedily add aligned node pairs (u _i , v _j ) to the network alignment, starting from the pairs with the highest node correspondence scores, until there is no other node pair left that can be added without creating inconsistencies in the network alignment. Assuming that the node correspondence scores in (8) obtained by the context-sensitive random walk model with restart accurately reflect the true correspondence between nodes - such that the score is proportional to the posterior node alignment probability - the proposed network alignment scheme can be viewed as a heuristic way to find the MEA alignment of the networks in G.

Datasets and experimental set-up

To assess the performance of the proposed method, we tested the proposed network alignment method based on PPI networks in NAPAbench [25] and IsoBase [26]. NAPAbench is a network alignment benchmark that consists of 3 different datasets, referred to as the pairwise alignment dataset, 5-way alignment dataset, and 8-way alignment dataset. Each dataset contains three different subsets of 10 network families, each subset created using a different network growth model - CG (crystal growth), DMC (duplication-mutation-complementation), and DMR (duplication with random mutation). Each network family consists of 2, 5, or 8 PPI networks depending on the alignment dataset. For network families in the pairwise alignment dataset, each family contains one network with 3,000 nodes and the other with 4,000 nodes. In the 5-way network alignment dataset, a network family consists of 5 networks with 1,000, 1,500, 2,000, 2,500, and 2,500 nodes. Finally, in the 8-way alignment dataset, every network family consists of 8 networks, where each network contains 1,000 nodes. To evaluate the performance of the proposed method on real PPI networks, we utilized IsoBase datasets [26], which was constructed by integrating the following databases: BioGRID [27], DIP [28], HPRD [29], MINT [30], and IntAct [31]. IsoBase contains the PPI networks of five species: H. sapiens, M. musculus, D. melanogaster, C. elegans, and S. cerevisiae. Currently, the PPI network of H. sapiens in [26] has 22,369 proteins and 43,757 interactions, the PPI network of M. musculus has 24,855 proteins and 452 interactions, the PPI network of D. melanogaster has 14,098 proteins and 26,726 interactions, the PPI network of C. elegans has 19,756 proteins and 5,853 interactions, and the PPI network of S. cerevisiae has 6,659 proteins and 38,109 interactions. In our analysis, we excluded the M. musculus network as it currently contains only a small number of interactions.

Based on our simulations, we report the following performance metrics: correct nodes (CN), specificity (SPE), mean normalized entropy (MNE), conserved interaction (CI), coverage, and computation time. CN is the total number of nodes in the correct equivalence classes. Given a network alignment, an equivalence class is defined as the set of aligned nodes, and if all nodes in the equivalence class have the same functionality the given equivalence class is said to be correct. SPE is the relative number of correct equivalence classes to the total number of equivalence classes in a network alignment. For each equivalence class C, the normalized entropy can be computed by $H\left(\text{C}\right)=-\frac{1}{\mathsf{\text{log}}d}{\sum }_{i=1}^d{p}_i\mathsf{\text{log}}{p}_i$, where p _i is the relative proportion of nodes in C with functionality i and d is the total number of different functionalities in the given equivalence class. As a result, a network alignment that accurately maps functionally similar nodes, hence being functionally consistent, will have lower mean normalized entropy. CI is defined as the total number of edges between equivalence classes. We also count the total number of edges between correct equivalence classes, which we refer to as the conserved orthologous interactions (COI), to assess the biological relevance of the conserved interactions that have been identified by the network alignment method. Finally, for 5-way and 8-way alignment datasets, we measure the equivalence class coverage and the node coverage, where the former is the number of equivalence classes that include nodes from k different networks, and the latter is the number of nodes in an equivalence class whose equivalence class coverage is k. For the performance evaluation based on real PPI networks in IsoBase, we determined the functionality of each protein using the KEGG protein annotation [32, 33]. Note that nodes without any functional annotation in each equivalence class and equivalence classes that consist of a single node or nodes from a single network were removed before computing the performance metrics.

We compared the performance of the proposed multiple network alignment method against a number of state-of-the-art algorithms: SMETANA [12], PINALOG [11], BEAMS [14], NetCoffee [13], and IsoRankN [7]. NetCoffee was not included in pairwise network alignment experiments, since it requires at least 3 networks. For multiple network alignment experiments, PINALOG was excluded as the algorithm can only handle pairwise alignments. For IsoRankN, we set the parameter α to 0.6 as in the original paper [7]. For BEAMS, we set the filtering threshold to 0.4 for IsoBase and 0.2 for NAPAbench as in the original paper [14], and set the parameter α to 0.5. The parameter α for NetCoffee was set to 0.5. We used default parameters for SMETANA (i.e., n_max = 10, α = 0.9, and β = 0.8), and the same parameters were used in the proposed network alignment method as well. Finally, in the proposed method, we used λ _t = 0.7 to determine whether or not to apply PCT to the estimated node correspondence scores.

All experiments were performed on a personal computer with a 2.4 GHz Intel i7 processor and 8 GB memory.

Performance assessment based on NAPAbench network alignment benchmark

As we can see in Table 1 in most cases, the proposed algorithm yields a significantly higher CN and SPE compared to other algorithms, which shows that the algorithm is capable of finding conserved nodes with both high sensitivity and specificity. Furthermore, the mean normalized entropy (MNE) is also much lower, indicating that the proposed algorithm yields network alignment results that are more functionally coherent. This table shows that BEAMS yields higher CN for the CG dataset, although its SPE is lower and its MNE is higher than the proposed method. Both SMETANA and the proposed algorithm shows similar performance on the CG dataset, but we can also see that the proposed algorithm consistently outperforms SMETANA on the DMC/DMR datasets. Multiple network alignment results obtained using the 5-way alignment dataset and the 8-way alignment dataset show similar trends. Tables 2 and 3 show that, in most cases, our proposed algorithm outperforms other algorithms with higher CN, higher SPE, and lower MNE. For multiple network alignment, we further compared different network alignment algorithms based on their capability of predicting equivalence classes that span all networks, since one of the main goals of multiple network alignment is to find functionally homologous proteins that are conserved in the networks of all target species. Simulation results show that, in most cases, our proposed method also yields much higher CN and SPE as well as lower MNE for equivalence classes that span all networks.

Next, we compare the number of conserved (orthologous) interactions identified by different network alignment algorithms. As Figure 1 shows, the proposed method was able to identify the largest number of conserved interactions as well as conserved orthologous interactions in most cases, resulting in higher CI and COI. The performance of SMETANA was comparable to the proposed method, while other algorithms typically resulted in lower CI and COI. It is worth noting that more than 95% of the conserved interactions that were detected by our proposed network alignment algorithm were between correct equivalence classes (i.e., conserved orthologous interactions). This certainly shows that our method can effectively detect biologically meaningful conserved interactions through network alignment.

We also analyzed the overall coverage of the predicted alignment results for the 5-way and 8-way network alignments. The results are shown in Figure 2 for the 5-way alignment and in Figure 3 for the 8-way alignment. For the 5-way network alignment, we can see that around 40% of the equivalence classes predicted by the proposed method contained nodes from all 5 networks. SMETANA shows a similar level of coverage, while for the remaining algorithms, only about 30% of the predicted equivalence classes included nodes from all 5 networks. The overall node coverage also shows similar trends. The 8-way alignment results summarized in Figure 3 show that the proposed algorithm can effectively find equivalence classes with good coverage, which include nodes from a large number of networks. For example, we can see that around 40% of the equivalence classes predicted by the proposed method contained nodes from all 8 networks.

Performance assessment based on protein-protein interaction networks in IsoBase

Figure 4 shows the computation time for aligning the PPI networks in IsoBase. SMETANA required the least computation time for pairwise network alignment and NetCoffee was the fastest among all for aligning the PPI networks of 3 species. Although IsoRankN yielded accurate alignment results for real PPI networks in IsoBase, it also required the largest amount of computation time in most cases. Figure 4 shows that our proposed network alignment algorithm requires relatively longer running time compared to other algorithms, in exchange for the improved alignment accuracy. Currently, the main bottleneck is the time required to construct the transition probability matrix $\tilde{P}$ of the context-sensitive random walker, and we are currently optimizing the code for our algorithm to make it computationally more efficient.