Seed selection strategy in global network alignment without destroying the entire structures of functional modules

Global network alignment

EC is the percentage of edges from a smaller network that are correctly aligned to edges in a bigger network. Naturally, when aligning two networks, we want to achieve as high EC as possible. GNA is NP-complete meaning as the underlying subnetwork isomorphism problem and heuristic approaches must be devised to get thus approximate solutions, especially for the large-sized PPI networks.

On the other hand maximizing the EC should not be the unique goal for biological networks alignment and we must strike a balance between topological and biological significances of our result. So our MHA improve the EC while maintaining the structure of functional modules. Although being the same as a seed-and-extend approach, MHA can solve the significant seeds chosen problem of such heuristics to a certain extent.

The seed selection strategy and our algorithm

Figure 2 shows an example of applying our method to the GNA problem. Firstly, we determine the centralities of nodes in each network and show the two networks in Figure 2(a) with the sizes of nodes proportional to their centrality values. Then the multiple hub seed set {(1, a), (6, k)} can be constructed for their local maximum similarities. It can be seen obviously that nodes {1, 6, a, k} are local hubs of networks respectively. Secondly, in Figure 2(b), beginning with the seed (1, a), the fractional alignment result {(1, a), (2, b), (3, c), (4, e), (5, d), (6, g)} is obtained as shown in Figure 2(b); Simultaneity, beginning with the seed (6, k), another fractional alignment result {(6, k), (7, l), (8, m), (9, n) } is obtained as shown in Figure 2(c). During this alignment process the dynamic membership similarity is changeable with the current aligned seed. MHA enables us to get a more overall consistent alignment shown as Alignment2 of Figure 1(c). MHA aligned two dense regions in network1 (marked with colours red and purple) according to distinct dense regions in network2 (marked with colours green and blue).

What follows is the detailed description of the major steps of the MHA method. In order to construct hub seeds for the seed-and-extend process, we should compute centralities and membership values of nodes in each network. An integrative network module identification and key nodes determination method family, called ModuLand [19], can be used to gain the values of the steps 1-3 straightway as follows:

Respectively, ${\tilde{C}}_G^i=\mathsf{\text{lg}}\left(C_G^i\right)$, ${\tilde{C}}_G^i=\mathsf{\text{lg}}\left(C_G^i\right)$ and θ is the average centrality value of all nodes in networks G and H. There are always plentiful pairs of nodes and the differences between their centrality values are not an order of magnitude. The logarithmic value $C_G^i$ and ${\tilde{C}}_H^j$ should be employed in practice. Note that SC(i, j) will be large if node i and j are both significant nodes ($C_G^i$ and ${\tilde{C}}_H^j$ are all higher) in their own network and has little gap $\left(\left|{\tilde{C}}_G^i-{\tilde{C}}_H^j\right|\mathsf{\text{ is small}}\right)$ of their centrality values. θ 's presence is to address the issue that a pair of similar nodes (i, j) with little gap should have a certain magnitude SC even if they are of very small roles in topology.

Step 2: Determination of multiple hub seeds. Here we present one concise approach suitable for the determination of hub seeds of networks. Local maxima based hub seeds is defined as: A hub seed is the pair of nodes having the locally maximal SC value, while all of their neighbouring nodes have lower SC values. The result of this step is the set SEEDS containing all hub seeds.

Functional modules in PPI are associated with a single central hub protein. Beginning with a hub (h _g , h _h ) ∈ SEEDS to extend and align the neighbourhood of the seed, this seed-and-extend process is to construct functional modules having similar biological function and topology around the hubs respectively. For every hub seed does the same process, MHA has avoided the phenomenon that multiple functional modules in G are aligned to a few dense regions in H and can align the networks by modules.

SM (i, j) is the harmonic mean of $L_{h_g}\left(i\right)$ and $L_{h_h}\left(j\right)$. It is not necessary to calculate SM (i, j) between all pairs of nodes (i, j) ∈ V _G × V _H relative to all hub seeds in H _G and H _H . For some current seed (h _g , h _h ), MHA only needs to get the SM values between neighbouring nodes around h _g and h _h during alignment process (in step 5) dynamically.

S is N × M -sized matrix and S(i, j) is an kind of dynamic similarity. Different seed considered among the alignment process leads to different SM (i, j) and S(i, j). The weight α can be adjusted to assign relative importance to biological and topological data, depending upon the confidence level attributed to them and the type of results sought. In our implementation we assign the weight (α = 0.6). This parameter has been discussed in the following section.

Then we present the detailed description of MHA based on the matrix S, and in the following part we define the specific concepts used.

Algorithm MHA (G, H)

Construct the matrices SC, SE and the set SEEDS.

Initialize alignment A to an empty set and alignment score vector B equals to 0.

for a hub seed (h _g , h _h ) ∈ SEEDS do

   Add (h _g , h _h ) to alignment A, B(h _g , h _h ) = S(h _g , h _h ).

   for all k ∈ { 1, ..., D} do

      Construct a bipartite graph $BP\left(N_G^k\left(h_g\right),N_H^k\left(h_h\right),E\right)$

      Compute SM (i, j) relative to the seed and assign the weight ω(i, j) = S (i, j).

      Solve the Maximum Weight Bipartite Matching Problem by the Hungarian algorithm [16].

      To each optimal matching (u, v) found above, if and only if S (u, v) >B(u), then add (u, v) to the current alignment A, B(u) = S(u, v).

   end for

end for

return alignment A.

The k ^th neighbourhood of node h _g in network G, $N_G^k\left(h_g\right)$, is defined as the set of nodes of G that are at distance ≥ k from h _g . Hence, $N_G^k\left(h_g\right)$ is the k ^th still unaligned neighbourhood and can be thought of as the "ball" of nodes around h _g up to and including nodes at distance k. This allows us to model insertions and deletions of nodes in the paths conserved between two networks. D is the longest distance restricted.

A bipartite graph, BP (V ₁, V ₂, E), is a graph with a node set V consisting of two partitions, V = V ₁ ∪ V ₂, so that every edge e ∈ E connects a node from V ₁ with a node from V ₂ ; that is, there are no edges between nodes of V ₁ and there are no edges between nodes of V ₂ -- all the edges "go across" the node partition.

The set A contains pairs of nodes which are the optimal matching results during the process of MHA. The matrix B records the alignment score when a matching is added to A. While using the multiple hub seeds, the object node in network H to which MHA has matched a node in G can be changeable. Membership values are related to the current aligned seeds. Beginning with different seeds, we will get different similarity score influenced by membership values, then a node in G should be aligned to the node that they can get the highest similarity score together. The significance of dynamical and changeable membership value is that a node can select the best associated hub seed to form a functional module, and also several smaller modules in one network will never be covered constrainedly by a larger one in another network along with the alignment process.

Computational complexity of MHA

NodeLand [19] algorithm for determination of centralities of nodes is structurally similar to a breadth-first search, therefore the worst-case runtime complexity is O(N*(N + |E|)), respectively, where N is the number of nodes and |E| is the number of edges in the network. The presented ProportionalHill [19] method has a runtime complexity of O(d*|E|*h), with d being the average node degree and h being the number of the identified hubs. Solving the alignment for bipartite graph BP (V ₁, V ₂, E) takes O((|V ₁|+|V ₂|)*(|E|+(|V ₁|+|V ₂|)*log((|V ₁|+|V ₂|)))). Therefore, the total time complexity of MHA algorithm for aligning networks G = (V _G , E _G ) and H = (V _H , E _H ) is smaller than O(h*|V _G |*(|E _G |+|V _G |*log(|V _G |))) and the space complexity is O(|V _G |*|V _H |+|E _G |+|E _H |) clearly.

Data sets and summarized statistics for our seed selection strategy

For each pair of networks table 2 shows our alignment results. When we align from the S. cerevisiae (S) to D. melanogaster (D) PPI networks, the set SEEDS consisting of 255 hub seeds is constructed. And there are 106 modules with obvious biological significance in yeast and 45 functional modules in fly formed around the seeds after the alignment process eventually. Performing our method for aligning M. musculus to H. sapiens PPI networks, the set SEEDS consisting of 91 hub seeds is obtained. There are 55 functional modules in mouse and 45 functional modules in human formed around the seeds after the alignment process eventually. All these functional modules have obvious biological process annotation of GO. The conserved modules shown in table 2 are those who execute the same functions in their respective networks and have the same annotation of GO terms in our alignment results. Conserved edge is the percentage of edges from a smaller network that are correctly aligned to edges in a bigger network. These results show that we can maintain the functional modules undamaged during the heuristic alignment process and align the networks by modules. Simultaneously our method is efficient in achieving a high percentage of conserved edge. In Figure 3, we give some functional modules of yeast that formed during our alignment process from yeast to fly PPI network. From a large-size seed some nodes painted same colour will be aligned. And the same colour nodes will form a functional module around the seed node.

	Hub seeds	Conserved edge	Functional modules of first species	Functional modules of second species	Conserved modules
S to D	255	2649(15.91%)	106(42%)	45(18%)	8(17.78%)
H to S	175	515(26.40%)	76(43%)	51(29%)	14(27.45%)
M to S	57	134(26.12%)	40(70%)	18(32%)	9(50%)
M to H	91	226(44.15%)	55(60%)	45(49%)	17(37.78%)
H to D	264	476(24.42%)	116(44%)	27(10%)	7(25.93%)

The best currently published network alignment algorithm, MI-GRAAL, uses a cost function relying on a highly constraining measure of topological similarity and achieved higher percentage of conserved edge (usually at about 50%) than our results. But the only one object for MI-GRAAL is to optimize the EC no fear of destroying the functional modules in original networks. For proteins always associate with each other to perform biological functions, it is more important to maintain the functional modules unchanged rather than maximize the EC value. As shown in Figure 4 we can improve our EC values shown in table 2 simply by decreasing the parameter α. And also our strategy of seed selection can be used as the reference and complement of other heuristic methods to seek more meaningful alignment results.

Enriched functional modules of alignment network

As a typical case study we use MHA to identify the common alignment network between PPI networks the S. cerevisiae (yeast) and D. melanogaster (fly). We construct the set SEEDS consisting of 255 hub seeds. The common graph corresponding to the global alignment between the yeast and fly PPI networks is comprised of 2649 edges with edge correctness of 15.08% (α = 0.6). The alignment graph consists of many disconnected components, with the largest component having 48 edges amongst 49 proteins (in Figure 5(a)). Table 3 presents the results of our algorithm and IsoRank [14], in terms of conserved interactions, size of the largest component and P-value. We know that MHA produces a more optimal solution, which can be found at no additional computational cost. Retrieving in particular 86% more conserved interactions than IsoRank for a given level of sequence similarity between two networks denotes that our alignment result is more topologically similar. Also, its largest component with P-value of 5.83E-13 has a significant function of DNA binding achieved congruously in both yeast and fly (in table 4). In contrast the function of IsoRank algorithm's largest component is unnoticed with P-value of 0.00098. As shown in Figure 5, the components discovered simultaneously span various topologies, from linear pathways (Figure 5(d), Figure 5(e)) to components corresponding to protein complexes (Figure 5(a), Figure 5(b) and Figure 5(c)). Each one has coherent function (identified by same GO term) in the two species. The detailed GO terms and P-value of these components are presented in table 4. Finally, we give some components having functionality only in yeast or fly detected by MHA in table 5. The most significant module with P-value of 4.30E-21, comprised of 49 aligned proteins, is obtained by our alignment in yeast PPI.

We emphasize that our components are rich in biological functions notwithstanding they are just subnetworks of the sparse alignment graph. Enrichment functions in our alignment graph profit from the selection of initialization multiple hubs. Around these hubs neighbourhood proteins form consistent and conserved functional modules in two PPI networks simultaneously. Our alignment result is more biologically similar.

Identification of functional orthologs

Effect of the parameter α

To measure the biological quality of the alignment result, we count the fraction of aligned pairs that have at least 1 GO term [22] in common, estimated by blue curve in Figure 4. Simultaneously we compute the edge correctness (EC) to evaluate the topological quality of the alignment result, marked by red curve. When choosing the most appropriate value of the free parameter α, we rejected the choice corresponding to the largest common alignment graph size (the largest EC); thus, the α leading to the largest-size alignment graph may not be a biologically appropriate choice. Instead, for the choice of α, we compared the biological quality with the topological quality and chose the α = 0.6 that is used to adjust the assignment of the relative importance to the topological and biological data.

The error tolerance of MHA

Our simulations indicate that the algorithm is tolerant to errors in the input networks (Figure 6); this is valuable since PPI networks have high false positive and false negative rates. To evaluate the algorithm's error-tolerance we randomized a fraction of its edges using the Maslov-Sneppen trick that preserves node degrees [24]: we randomly choose two edges (a, b) and (c, d), remove them, and introduce new edges (a, d) and (c, b). For each choice of the percentage of randomized edges, we compute the fraction of nodes that are mapped to exactly same objects in the original graph and randomized graph after GNA. Clearly, the algorithm's performance degrades smoothly and very slowly. These simulations suggest our results are quite robust to errors in PPI data.