A fast tool for minimum hybridization networks

Simulated data

To generate simulated datasets, we use a program due to Beiko and Hamilton [12]. To obtain a pair (T, T’) of trees, their program first generates T randomly and then obtains T’ from T by performing a specified number r (say, 20) of random rSPR operations on T. Recall that an rSPR operation on a tree T first removes an edge (p,c) from T, then contracts p (the vertex of out-degree 1 resulting from the removal of edge (p,c)), and further re-attaches the subtree rooted at c to an edge (p’,c’) of T (by introducing a new vertex m’, splitting edge (p’,c’) into two edges (p’,m’) and (m’,c’), and adding a new edge (m’,c)). So, the actual rSPR distance of T and T’ is at most r. Moreover, the hybridization number of T and T’ can be r, smaller than r, or larger than r.

We first use Beiko and Hamilton’s program to generate 120 pairs of trees each of which has 100 leaves. The first (respectively, second) 60 pairs are generated by setting r = 14 (respectively, r = 17). It turns out that among the 120 generated tree-pairs, 6 (respectively, 22, 33, 11, 21, or 27) tree-pairs have hybridization number 12 (respectively, 13, 14, 15, 16, or 17). Figure 1 summarizes the average running time of the programs for the generated tree-pairs, where each average is taken over those tree-pairs with the same hybridization number. As can be seen from the figure, FastHN with Heuristic 1 and/or Heuristic 2 is much faster than HybridNet and Dendroscope 3. This difference in speed becomes more significant as the hybridization number becomes larger. Moreover, Heuristic 1 contributes the most to the saving of running time. Indeed, when Heuristic 1 is used, both Heuristic 2 and initial cluster reductions do not help much. It is worth noting that Beiko and Hamilton’s program tends to create a pair of trees without a relatively large common clusters. This is why initial cluster reductions do not help much for tree-pairs randomly generated by their program.

The comparison is done on 120 randomly generated tree-pairs with relatively small hybridization numbers, where each tree has 100 leaves. If the running time of a program for a tree-pair is more than 600 seconds, then it has been rounded down to 600 seconds. Among the 120 pairs, Dendroscope 3 takes more than 600 seconds for 28 pairs, FastHN without Heuristic 1, 2, or initial cluster reductions takes more than 600 seconds for 20 pairs, FastHN with only initial cluster reductions takes more than 600 seconds for 18 pairs, and FastHN with Heuristic 1 or 2 takes less than 60 seconds for every pair.

To test how the number of leaves in an input tree influences the running time of the algorithms, we next use Beiko and Hamilton’s program to generate 120 pairs of trees each of which has 50 leaves. The first (respectively, second) 60 pairs are generated by setting r = 14 (respectively, r = 17). It turns out that among the 120 generated tree-pairs, 3 (respectively, 17, 26, 26, 20, 20, or 6) tree-pairs have hybridization number 11 (respectively, 12, 13, 14, 15, 16, or 17). Moreover, for each h∈{9,10}, there is exactly one generated tree-pair with hybridization number h. Figure 2 summarizes the average running time of the programs for those generated tree-pairs with hybridization number in the range [12 .. 17], where each average is taken over those tree-pairs with the same hybridization number. As can be seen from the figure, the superiority of FastHN over HybridNet and Dendroscope 3 remains the same (as in Figure 1) if Heuristic 1 or 2 is used. Moreover, Heuristic 1 contributes the most to the saving of running time.

The comparison is done on 120 randomly generated tree-pairs with relatively small hybridization numbers, where each tree has 50 leaves.

To compare the performance of the algorithms for tree-pairs with relatively large hybridization numbers, we further use Beiko and Hamilton’s program to generate 60 pairs of trees by setting r = 25, where each tree has 100 leaves. It turns out that among the 60 generated tree-pairs, 4 (respectively, 10, 15, 18, or 12) tree-pairs have hybridization number 21 (respectively, 22, 23, 24, or 25). Moreover, there is exactly one generated tree-pair with hybridization number 20. Figure 3 summarizes the average running time of the two best versions of FastHN for the generated tree-pairs, where each average is taken over those tree-pairs with the same hybridization number. As can be seen from the figure, both versions take less than 16 minutes (on average) even when the hybridization number is as large as 25, while FastHN with both Heuristics 1 and 2 is the faster version. In contrast, Dendroscope 3 fails to complete for each of the 60 tree-pairs.

The time is measured on 60 randomly generated tree-pairs with relatively large hybridization numbers, where each tree has 100 leaves. (Note: Dendroscope 3 fails to complete for each of the 60 pairs.)

Biological data

As can be seen from Table 1, for most of the tree-pairs, there is not much difference in speed between Dendroscope 3 and FastHN with initial cluster reductions. This is because most of the tree-pairs have small hybridization numbers. For the tree-pair (ndhf, ITS), FastHN with cluster reductions runs substantially faster than FastHN without cluster reductions. This is because the problem for this pair can be reduced to two tree-pairs of roughly equal sizes by initial cluster reductions in the preprocessing stage of the algorithm.

Discussion

Roughly speaking, FastHN consists of two stages, namely, the preprocessing stage and the exhaustive-search stage. In the preprocessing stage, FastHN performs only subtree reductions and cluster reductions. Indeed, other kinds of reductions are also known. One of them is chain reduction [15]. Performing chain reductions on the input trees results in trees whose nodes are weighted. Unfortunately, it seems that Whidden et al.’s O(2.4^2dn)-time algorithm for computing the rSPR distance d of two given phylognetic trees with n leaves does not work when the trees are weighted. This is why FastHN does not perform chain reductions.

In the exhaustive-search stage, FastHN also performs subtree reductions whenever possible, but does not perform cluster reductions. The main reason of not performing cluster reductions in the 2nd stage is that performing a cluster reduction is too time-consuming (namely, takes O(^{n 2}) time, where n is the number of leaves in the trees).

When running FastHN, one can decide whether to adopt initial cluster reductions, Heuristic 1, or Heuristic 2. If two input trees have relatively large common clusters, performing initial cluster reductions on them lead to solving independent and significantly smaller subproblems. So, we should always choose to adopt initial cluster reductions. Moreover, as can be seen from our simulated results, we should always choose to adopt Heuristic 1 because it enables the algorithm to save a lot of time by cutting more branches of the search tree in the exhaustive-search stage. Our simulated results also show that adopting both Heuristics 1 and 2 makes FastHN run faster (on average) than adopting only Heuristic 1. Thus, in general, we should choose to adopt Heuristic 2 as well. However, in our experiments, we have found some tree-pairs for which FastHN with Heuristic 1 but without Heuristic 2 runs significantly faster than FastHN with both Heuristics 1 and 2. Hence, as long as Heuristic 1 is adopted, there is still room to decide whether to adopt Heuristic 2 as well.

Hybridization networks and phylogenetic trees

Let X be a set of existing species. A hybridization network on X is a directed acyclic graph N in which the set of nodes of out-degree 0 (still called the leaves) is X, each non-leaf node has out-degree 2, there is exactly one node of in-degree 0 (called the root), and each non-root node has in-degree larger than 0. Note that the in-degree of a non-root node in N may be larger than 1. A node of in-degree larger than 1 in N is called a reticulation node of N. Intuitively speaking, a reticulation node corresponds to a reticulation event. The hybridization number of a reticulation node in N is its in-degree in N minus one. The hybridization number of N is the total hybridization number of reticulation nodes in N.

A phylogenetic tree on X is a rooted binary tree whose leaf set is X. A hybridization network N on X displays a phylogenetic tree T on X if N has a subgraph M such that M is a rooted tree, the root of M has exactly two children in M, and modifying M by contracting its unifurcate nodes yields T. A hybridization network of two phylogenetic trees T₁ and T₂ on X is a hybridization network N on X such that N displays both T₁ and T₂. A hybridization network of T₁ and T₂ is minimum if its hybridization number is minimized among all hybridization networks of T₁ and T₂. Obviously, if N is a minimum hybridization network of T₁ and T₂, then the in-degree of every reticulation node in N is exactly 2 and hence the hybridization number of N is equal to the number of reticulation nodes in N. For convenience, we define the hybridization number of T₁ and T₂ to be the minimum hybridization number of a hybridization network of T₁ and T₂.

We are now ready to define one problem studied in this paper:

Hybridization Network Construction (HNC):

· Input: Two phylogenetic trees T₁ and T₂ on the same set X of species.

· Goal: To construct a minimum hybridization network of T₁ and T₂.

Agreement forests

Throughout this subsection, let T₁ and T₂ be two phylogenetic trees on the same set X of species. If we can apply a sequence of detaching operations on each of T₁ and T₂ so that they become the same forest F, then we refer to F as an agreement forest (AF) of T₁ and T₂. A maximum agreement forest (MAF) of T₁ and T₂ is an agreement forest of T₁ and T₂ whose size is minimized over all agreement forests of T₁ and T₂. The size of an MAF of T₁ and T₂ is called the rSPR distance between T₁ and T₂. The following lemma is shown in [16].

Extending Whidden et al.’s Algorithm

Throughout this subsection, let T₁ and T₂ be two phylogenetic trees on the same set X of species. We sketch the fastest known algorithm (due to Whidden et al.[16]) for computing an MAF of T₁ and T₂, and then state a slight extension of the algorithm that will be used in our algoirthm for HNE.

The basic idea behind Whidden et al.’s algorithm is as follows. For k = 0, 1, 2, …(in this order), we try to find an AF of T₁ and T₂ of size k and stop immediately once such an AF is found. To find an AF of T₁ and T₂ of size k, we start by setting _{F 1}=_{T 1} and _{F 2}=_{T 2} and associating a label set{x} to each leaf x of F₁ and F₂. We then repeatedly modify F₁ and F₂ (until either |_{F 1}|>k or F₁ becomes a forest without edges) as follows. We find two arbitrary sibling leaves u and v in F₂. If u and v are also siblings in F₁, then we modify F₁ and F₂ separately by merging the identical subtrees of F₁ and F₂ rooted at the parent of u and v each into a single leaf whose label set is the union of the label sets of u and v. On the other hand, if u and v are not siblings in F₁, then we distinguish three cases as follows. Case 1: u and v are in different component trees of F₁. In this case, in order to transform F₁ and F₂ into an AF of T₁ and T₂, we have two choices to modify them, namely, by either detaching the subtree rooted at u or detaching the subtree rooted at v. Case 2: u and v are in the same component tree of F₁ and either (1) u and the parent of v are siblings in F₁ or (2) v and the parent of u are siblings in F₁. In this case, if (1) (respectively, (2)) holds, then we modify F₁ by detaching the subtree rooted at the sibling of v (respectively, u). Case 3: u and v are in the same component tree of F₁ and neither (1) nor (2) in Case 2 holds. In this case, in order to transform F₁ and F₂ into an AF of T₁ and T₂, we have three choices to modify them. The first two choices are the same as those in Case 1. In the third choice, we modify F₁ by detaching the subtrees rooted at those non-root nodes w such that the parent of w appears on the (not necessarily directed) path between u and v in F₁ but w does not.

By the above three cases, we always have the following:

· |_{F 1}|≥|_{F 2}|.

· All component trees of F₂ except at most one have no edges.

· For each component tree _{Γ 2}of F₂ without edges, F₁ has a component tree _{Γ 1}without edges such that the label sets associated with the unique leaves of _{Γ 1}and _{Γ 2}are identical.

Once |_{F 1}| becomes larger than k, we know that F₁ and F₂ have no AF of size k. On the other hand, once F₁ becomes a forest without edges, we can use the label sets L(v) of the leaves v of F₁ to obtain an AF of T₁ and T₂ of size |_{F 1}| by modifying T₁ as follows. For each leaf v of F₁ such that L(v) does not contain the dummy leaf, detach the subtree of T₁ rooted at the LCA of the leaves in L(v).

Now, we are now ready to make a key observation in this paper. By (b) and (c) in the above, Whidden et al.’s MAF algorithm can actually be used to solve the following slightly more general problem in O(2.4^2kn) time:

rSPR Distance Checking (rSPRDC):

· Input:(_{T 1},_{T 2},k,_{F 1},_{F 2}), where T₁ and T₂ are two phylogenetic trees on the same set X of species, k is an integer, F₁ (respectively, F₂) is a rooted forest obtained from T₁ (respectively, T₂) by performing zero or more detaching operations, and every component tree of F₂ except at most one is identical to a component tree of F₁.

· Goal: To decide if performing k more detaching operations on F₁ leads to an AF of T₁ and T₂.

Finally, if we want to enumerate all MAFs of T₁ and T₂, then we need to modify Whidden et al.’s algorithm as follows. First, we do not distinguish Cases 2 and 3 because modifying F₁ as in Case 2 may lose some MAF of T₁ and T₂. Moreover, whenever an AF of T₁ and T₂ of size k is found, we do not stop immediately and instead continue to find other AFs of T₁ and T₂ of size k. The resulting algorithm runs more slowly, namely, in O(^3kn) time.

Speeding up HybridNet

Throughout this subsection, let T₁ and T₂ be two phylogenetic trees on the same set X of species. We first sketch how HybridNet enumerates all MAAFs of T₁ and T₂, and then explain how to speed it up.

First, we need several definitions. For a rooted forest F, we use (F) to denote the family of the leaf sets of the component trees of F. Let F and F’ be two forests each obtained by performing zero or more detaching operations on T₁. If F≠ ^{F ″} and for every set $Y\in \mathcal{ℒ}\left(F\right)$, there is a set $Y^{″}\in \mathcal{ℒ}\left(F^{″}\right)$ with $Y\subseteq Y^{″}$, then we say that F is finer than F’ and F’ is coarser than F.

To enumerate all MAAFs of T₁ and T₂, the idea behind HybridNet is to design an algorithm for the following problem:

Generalized Agreement Forest (GAF)

· Input:(_{T 1},_{T 2},k,_{F 1}), where T₁ and T₂ are two phylogenetic trees on the same set X of species, k is an integer, and F₁ is a rooted forest obtained from T₁ by performing zero or more detaching operations.

· Goal: To find a sequence of AFs of T₁ and T₂ including all AFs F of T₁ and T₂ such that (1) F can be obtained by performing at most k detaching operations on F₁ (or equivalently, at most |_{F 1}| + k detaching operations on T₂) and (2) no AF of T₁ and T₂ is finer than F₁ and coarser than F.

In the supplementary material of [19], an O(^3kn)-time algorithm for solving GAF is detailed. The algorithm differs from Whidden et al.’s algorithm for enumerating all MAFs of T₁ and T₂ only in that we start with F₁ (as it is given) and _{F 2}=_{T 2}(instead of starting with _{F 1}=_{T 1}and _{F 2}=_{T 2}) and then repeatedly modify F₁ and F₂ until either |_{F 1}|>k + _{k 0}or F₁ becomes a forest without edges, where k₀ is the original size of F₁. Now, we are now ready to make two other key observations in this paper. To speed up Chen and Wang’s algorithm for solving GAF, we modify it as follows:

· Heuristic 1: Every time before we start to make multiple choices of modifying F₁ and F₂, we call the algorithm for rSPRDC in Lemma 1 on input (_{T 1},_{T 2},k−|_{F 1}| + _{k 0},_{F 1},_{F 2}) to check if performing k−|_{F 1}| + _{k 0}more detaching operations on F₁ leads to an AF of T₁ and T₂.

As the result, if we know that performing k−|_{F 1}| + _{k 0}more detaching operations on F₁ does not lead to an AF of T₁ and T₂, then no more choice of modifying F₁ and F₂ is necessary; otherwise, we proceed to make multiple choices of modifying F₁ and F₂ the same as before but with the following difference:

· Heuristic 2: Instead of selecting two arbitrary sibling leaves u and v in F₂ (cf. the Extending Whidden et al.’s Algorithm section), we select two sibling leaves u and v in F₂ such that they are as far apart as possible in F₁.

The intuition behind Heuristic 2 is that if u and v are far apart in F₁, then either u and v fall into two different connected components of F₁ so that we do not have to try Case 3 in the Extending Whidden et al.’s Algorithm section, or u and v fall into the same connected component of F₁ and we can detach a lot of subtrees from F₁ in Case 3.

Note that Step 11b is nontrivial. As described in the supplementary material of [19], Lemma 2 is very helpful for this purpose. More specifically, we first find a cycle C in G_F″ in O(|F″|²) time. By Lemma 2, in order to make F’ acyclic, we have to choose one node r of C and modify F’ by detaching the subtree of F’ rooted at an (arbitrary) child of r. Note that since r is a root of F’, detaching the subtree of F’ rooted at a child of r is achieved by simply deleting r from F’ and is hence independent of the choice of the child. Moreover, if the parent r’ of the dummy leaf in F’ is a node of C, then by Lemma 3, we can exclude r’ from consideration when choosing r. So, we have at most | ^{F ″} |≤k−1 ways to break C. After modifying F’ in this way, we again construct G _F’ and test if it is acyclic. If it is acyclic, then we can output F’; otherwise, we again find a cycle C in G _F’ and use it to modify F’ as before. We repeat modifying F’ in this way, until either F’ becomes acyclic, or | ^{F ″} |=k and G _F’ is still cyclic. Once F’ becomes acyclic, we output it. The total time taken by Step 11b is O(²(k−1)^k−|F″|), because we make a total number of at most O((k−1)^k−|F″|) choices for breaking cycles.

Experiments show that Heuristics 1 and 2 help us speed up the algorithm substantially. However, the two heuristics may not help in the worst case. That is, we are unable to prove that the two heuristics improve the worst-case time complexity of the algorithm which is O(^3d|X| + ^{3d(k−1)k−d + 2}) (as shown in [19]), where d is the size of an MAF of T₁ and T₂. We note that k and d are usually quite close.

The new algorithm for HNE

In this subsection, we only design an algorithm for HNE. Note that it is trivial to obtain a faster algorithm for HNC by modifying the algorithm for HNE so that it stops immediately once an MAAF is found.

Throughout this subsection, let T₁ and T₂ be two phylogenetic trees on the same set X of species. As mentioned before, we can easily use an MAAF of T₁ and T₂ to construct a hybridization network displaying T₁ and T₂[6]. So, we only explain how to enumerate all MAAFs of T₁ and T₂.

In the last subsection, we have explained how to speed up HybridNet so that it can enumerate all MAAFs of T₁ and T₂ within shorter time. Indeed, we can make HybridNet even faster. The idea is to preprocess T₁ and T₂ so that the given trees become smaller or the problem becomes to solve two or more smaller independent subproblems. More specifically, we perform the following two reductions on T₁ and T₂ until neither of them is available.