Fast calculation of the quartet distance between trees of arbitrary degrees

Background A number of algorithms have been developed for calculating the quartet distance between two evolutionary trees on the same set of species. The quartet distance is the number of quartets – sub-trees induced by four leaves – that differs between the trees. Mostly, these algorithms are restricted to work on binary trees, but recently we have developed algorithms that work on trees of arbitrary degree. Results We present a fast algorithm for computing the quartet distance between trees of arbitrary degree. Given input trees T and T', the algorithm runs in time O(n + |V|·|V'| min{id, id'}) and space O(n + |V|·|V'|), where n is the number of leaves in the two trees, V and V are the non-leaf nodes in T and T', respectively, and id and id' are the maximal number of non-leaf nodes adjacent to a non-leaf node in T and T', respectively. The fastest algorithms previously published for arbitrary degree trees run in O(n3) (independent of the degree of the tree) and O(|V|·|V'|·id·id'), respectively. We experimentally compare the algorithm with existing algorithms for computing the quartet distance for general trees. Conclusion We present a new algorithm for computing the quartet distance between two trees of arbitrary degree. The new algorithm improves the asymptotic running time for computing the quartet distance, compared to previous methods, and experimental results indicate that the new method also performs significantly better in practice.


Background
The evolutionary relationship for a set of species is conveniently described by a tree in which the leaves correspond to the species, and the internal nodes correspond to speciation events. The true evolutionary tree for a set of species is rarely known, so inferring it from obtainable information is of great interest. Many different methods have been developed for this, see e.g. [1] for an overview.
Different methods often yield different inferred trees for the same set of species, and even the same method can give rise to different evolutionary trees for the same set of species when applied to different information about the species. To study such differences in a systematic manner, one must be able to quantify differences between evolutionary trees using well-defined and efficient methods.
One approach for comparing evolutionary trees is to define a distance measure between trees and compare two trees by computing this distance. Several distance measures have been proposed, e.g. the symmetric difference metric [2], the nearest-neighbour interchange metric [3], the subtree transfer distance [4], the Robinson and Foulds distance [5], and the quartet distance [6]. Each distance measure has different properties and reflects different aspects of biology.
This paper is concerned with calculating the quartet distance. A quartet is a set of four species, the quartet topology induced by an evolutionary tree is determined by the minimal topological subtree containing the four species. The four possible quartet topologies of four species are shown in Fig. 1. Given two evolutionary trees on the same set of n species, the quartet distance between them is the number of sets of four species for which the quartet topologies differ in the two trees.
Steel and Penny [7] pointed at Doucettes unpublished work [8] which presented an algorithm for computing the quartet distance in time O(n 3 ), where n is the number of species. Bryant et al. in [9] presented an improved algorithm which computes the quartet distance in time O(n 2 ) for binary trees. Brodal et al. in [10] showed how to compute the quartet distance in time O(n log n) considering binary trees. For arbitrary degree trees, the quartet distance can be calculated in time O(n 3 ) or O(n 2 d 2 ), where d is the maximum degree of any node in any of the two trees, as shown by Christiansen et al. [11].

Results and discussion
In [11], we presented an algorithm for computing the quartet distance between trees of arbitrary degree. It runs in time O(n 2 d 2 ) and space O(n 2 ), where n is the number of leaves in each tree and d is the maximal degree found in either of the trees. In this paper, we present an improved algorithm running in time O(n + |V||V'| min{id, id'}) and space O(n + |V||V'|), where |V| and |V'| are the number of internal (non-leaf) nodes in the two input trees, and id and id' are the maximal degree of an internal node, when disregarding edges to leaves, in the two trees.

Time analysis for different types of trees
The terms |V|, id, |V'| and id' are all clearly O(n), but on the other hand neither |V| and id nor |V'| and id' are independent. Intuitively, if there are a lot of internal nodes in a tree, they will not have a very large internal degree. We address in this section, how this dependency will affect the running time for different types on trees.
The worst theoretical running time of the algorithm for calculating the quartet distance presented above is O(n 3 ).
Consider a tree with an internal node of degree , connected to internal nodes of degree three each con-  calculating the quartet distance will be O (n 2 ). In the following section we will do practical verification of the theoretical results in this section.

Experimental running times
The graphs in Fig. 3 show the running time for comparing worst case trees (see Fig. 2), (d-ary trees and random trees. There are six types of (d-ary trees; binary, 6-ary, 15-ary, and 30-ary and two types of random trees; r8s-based (see [12]) and trees with random topologies. The trees generated by r8s are binary, but by contracting edges, we can get trees of arbitrary degree (contracting an edge e connecting nodes u and v means removing u and e and attaching the rest of u's edges to v). Each edge is contracted with a probability that is inversely proportional with its length, i.e. a short edge has a higher probability of being contracted than a long edge. The trees with random topology are gen-erated by adding leaves one by one, starting with a tree of size 2. A leaf can be added by attaching it to a random inner node or by spliting a random edge with a new node, to which the leaf is attached.
The running time for worst case input trees (as described in the previous section) is O(n 3 ), because such trees have O(n) internal nodes and min{id, id'} is O(n). This is supported by the first graph in Fig. 3, which shows that the plot of the polynomial n 3 (representing the best sum-ofsquares fit of the polynomial c·n 3 to the data-points) is closest to the plot of the running times with regard to slope.
The running time on the algorithm on d-ary trees is

O(
). The plots of the running times in the second graph are parallel, and one of them is plotted directly on top of a plot of the polynomial n 2 (here c·n 2 is fitted to the data-points for each d separately; the different colors match the d colors). This supports that they all have a running time of O(n 2 ) for fixed d's. The graph also shows that higher degrees give lower running times, which is also expected. The reason why the algorithm is more than twice as fast on 6-ary trees than it is on binary trees, is that the number of internal nodes in 6-ary trees is less than in binary trees, and even though |V| is O(n) in both cases, that difference has an impact on the running time. The last graph shows the running time of the algorithm on trees created as either random trees (each topology is equally likely) or trees simulated using r8s (with edge contraction as described above). We have no theoretical running time for this data, but the graphs show that the running time is O(n 2 ). Even though the plotted data is only a small random sample, this indicates that many pairs of trees actu- Therefore it is not unreasonable to expect that our algorithm runs in time O(n 2 ) on trees used in practice. All experiments were performed on a standard PC (Pentium 4, 3 GHz, 1 Gb Ram) running Linux Fedora Core 3.

Comparison with existing algorithms
In Fig. 4 we compare the running time of the new algorithm with the O(n 2 d 2 ) and O(n 3 ) time algorithms from [11] on random and r8s simulated trees. In Fig. 5 we compare the running time of the new algorithm with the other two algorithms on Buneman and refined Buneman trees built for a range of Pfam [13] derived distance matrices using the tool in [14]. Buneman and refined Buneman trees are not binary unless this is well supported by the input distance matrix, and thus represent the kind of trees   The observed running times on random trees thus indicate that over selection of random trees consists of trees of low degree, whereas the r8s simulated, Buneman, and refined Buneman trees are trees with a few nodes of high degree which more than compensate for the additional overhead of dealing with nodes of low degree.
In conclusion, we find that the experimental comparison of the new algorithm with the previously developed algorithms indicate that the new algorithm not only improves on the theoretical asymptotic running time, but also improves the running time in practice if the input trees contain a few nodes of high degree.

Conclusion
We have constructed an algorithm for finding the quartet distance between two trees of arbitrary degree. It

Comparison of new and existing algorithms
Number of leaves These theoretical running times have been validated by running a series of tests using a Java implementation of the algorithm, available at [15]. We also done a series of tests on random trees, trees generated by the program r8s, Buneman trees, and refined Buneman trees. Running the algorithm on these trees gives an impression on how it performs on trees used in practice. On both types of trees the running time appears to be O(n 2 ). It is however still an open problem to develop an algorithm running in time O(n 2 )for all types of trees.

Methods
Consider two input trees, and assume that a quartet has butterfly topology in both trees, i.e. that one pair of the four leaves is separated from the other pair by an edge in the tree in both trees. We say that the butterfly quartet is shared, if it has the same butterfly topology in both trees. Otherwise, we say that the butterfly quartet is nonshared. We let shared(T, T') denote the number of butterflies shared between tree T and tree T', i.e. the number of quartets that are butterflies with the same topology in tree T and tree T', and let nonshared (T, T') denote the number of quartets that are butterflies in both T and T' but with different topology. By our definition of shared, the number of butterfly quartets in a single tree can be stated as the number of butterfly quartets shared between the d min Comparisons with earlier algorithms on Buneman and refined Buneman trees 10 20 50 100 200 500 1000 5e−01 5e+00 5e+01 5e+02

Comparison on Buneman and refined Buneman trees
tree and itself, i.e. shared(T, T) or shared(T', T') for the number of butterfly quartets in T and T' respectively. (This notation also emphasizes that computing the number of butterfly quartets in a single tree by our algorithm is performed as a comparison of the tree against itself.) In [11] we argue that the quartet distance between T and T', qdist(T, T'), can be found by focusing only on the computation of the number of shared and nonshared butterfly quartets between two trees, i.e. it is unnecessary to consider non-butterfly quartets explicitely. More specifically, we show that: The proof of this formula is a follows. Let Q denote the number of quartets which have butterfly topology in T and non-butterfly topology in T'. Symmetrically, let Q' denote the number of quartets which have butterfly topology in T' and non-butterfly topology in T. A butterfly quartet in T is either a butterfly quartet in T' or a non-butterfly quartet in T'. The number of butterfly quartets in T, shared(T, T), can thus be expressed as the sum shared(T, T') + nonshared(T, T') + Q. Similarly, the sum shared (T', T') = Q' + shared (T, T') + nonshared (T, T'). It is now straightforward to verify that the righthand side of (1) adds up Q + Q' + nonshared(T, T') which is the number of quartets where the quartet topologies differ in T and T', i.e. qdist(T,T').
In the section below, we describe how to use (1)

Terminology
Let T and T' be two unrooted trees. In this paper we will explicitly refer to the leaves of a tree as leaves and the nonleaf nodes as internal node. We will assume that T and T' each has n labelled leaves numbered 1,..., n such that the leaf numbered x in T has the same label as the leaf numbered x in T'. The leaf sets are denoted L and L' for T and T' respectively, note that L = L'. We will use V and V' to denote the internal nodes in T and T' respectively. The degree of an internal node v is the number of subtrees connected to it, and is denoted d v . The internal degree of an internal node v, id v , is the number of non-leaf subtrees connected to it. We will assume that no internal node in T and T' has degree two, and we will denote the maximal internal degree of all internal nodes in T and T' by id and id' respectively. Let v be an internal node in T, and let F 1 , ..., be the subtrees connected to it, as shown in Fig. 6 We call these the subtrees of v. We say that v claims all butterfly quartets ab|cd where a,b ∈ F i , c ∈ F k and d ∈ F m for i ≠ k ≠ m (see Fig. 7). With this definition, each butterfly quartet is claimed by exactly two internal nodes.
Adding the subscript yz|w|x to an internal node claiming the butterfly quartet wx|yz, indicates that the leaves y and z are found in a single subtree of the internal node, while the leaves w and x are found in different subtrees. For example, considering the quartet ab|cd, v and v' in Fig. 7 are written as v ab|c|d and .
Given a subtree F of T, and a subtree G of T', we call the intersection F ∩ G a shared leaf set, i.e. the set of leaves present in both F and G. The size of the shared leaf set, |F ∩ G|, then denotes the number of leaves present in both F and G. The size of a single subtree F is similarly denoted |F|. We will use to represent the subtree of T containing all leaves not in F and similarly for and G in T', see

T T T T T T T T
An internal node v ∈ T with subtrees F 1 ,...,F d , here d v = 6.
entation of the algorithms below, we will assume that we have access to |F|, |G| and |F ∩ G| for all subtrees F of T and G of T' in time O (1). At the end of the section we will describe how this can be achieved by an O(n) time preprocessing step, which does not affect the asymptotic worst case running time of the presented algorithms.

Counting shared butterfly quartets
For each pair of internal nodes v, v' from T, T' we want to count the number of shared butterfly quartets claimed by both internal nodes, shared(v, v'). Assume that F 1 , ..., are subtrees of v and G 1 , ..., are subtrees of v'. We wish to count all quartets on the form ab|cd where a, b ∈ F i ∩ G j , c ∈ F k ∩ G l and d ∈ F m ∩ G n , i ≠ k ≠ m, j ≠ l ≠ n (see Fig.   7). Counting the possible combinations of a and b is expressed by the following double sum, which sums over all pairs of subtrees of v and v': Given that a and b are in F i ∩ G j , we need to find c and d Any pair in F k ∩ G l is counted twice, once in |F k ∩ j | and once in | i ∩ G l |, therefore these pairs are subtracted once using the double sum above. (2) expresses the number of ways c and d can be found in different subtrees, given that a and b are found in F i ∩ G j : A rooted subtree F, and its complement rooted subtree Figure 8 A rooted subtree F, and its complement rooted subtree .

F F
Internal nodes v ∈ T and v' ∈ T', each claiming the quartet ab|cd Figure 7 Internal nodes v ∈ T and v' ∈ T', each claiming the quartet ab|cd.
We can now compute the number of shared butterfly quartets between two internal nodes, ie. the number of butterfly quartets claimed by both internal nodes with the same topology: If the trees, T and T', have a shared quartet ab|cd, then there are two internal nodes in each tree that claims this quartet: v ab|c|d and v cd|a|b in T and and in T'.
Since both shared(v ab|c|d , ) and shared(v cd|a|b , ) will count the quartet, the total number of shared quartets between the two trees is: It is straightforward to observe that calculating shared(v, Provided that the sums j , , S j , and S have been calculated, (4) can be calculated in time O(l). Since calculation of the sums is independent on the calculation of shared(v, v'), these calculations can be done serially as shown in the algorithm below, thereby reducing the time usage of shared(T, T') to:

ALGORITHM -CALCULATING THE NUMBER OF SHARED BUTTERFLY QUARTETS BETWEEN T AND T'
Requires: T, T' two input trees with the same leaf set.

Ensures: Res = shared(T, T')
Res ← 0 for v internal node in T do

end for end for
Res ←

Counting nonshared butterfly quartets
For each pair of internal nodes v, v' we want to count the number of nonshared butterfly quartets claimed by both internal nodes, nonshared(v, v'). Such quartets have the property that a pair of leaves found in the same subtree of v will be found in different subtrees of v' and vice versa, i.e. a nonshared quartet with leaves a, b, c and d, has a ∈ F i ∩ G j , b ∈ F i ∩ G l , c ∈ F k ∩ G n and d ∈ F m ∩ G j (see Fig. 9). The following expression counts all nonshared quartets related to a pair of nodes v and v', obeying that if two leaves of the quartet are in one subtree of v they are in different subtrees of v' and vice versa: Even though (6) satisfies the property of nonshared quartets, it possibly counts more than the number of nonshared quartets claimed by an internal node in each tree. The problem is that given two internal nodes, they do not nescessarily claim the quartets counted by (6). If we denote the leaves of an nonshared quartet a, b, c and d, the first, second, third and fourth factors in (6)  Similarly there are cases where b and c are chosen from the same subtree G l , l ≠ j of v', which we must also subtract. These can be counted as: The cases where both c and d are chosen from the same subtree F k , k ≠ i of v and b and c are chosen from the same subtree G l , l ≠ j of v' are included in both the expressions above and therefore they must be added again. The following expression counts the number of these cases: Combining equations (6), (7), (8) and (9), gives a way of calculating the number of nonshared quartets between two internal nodes v and v': Assuming that the trees have a nonshared quartet with topology form ab|cd in T, and ad|bc in T', there are two internal nodes in each tree claiming the quartet: v ab|c|d and v cd|a|b in T and and in T'. All of the four combinations of these will identify the quartet as nonshared. Therefore the total number of nonshared quartets between the two trees is: Each of the subtrees has to be intersected with two disjoint sets in each of the sums. This means that if a subtree is only a leaf, at least one of these intersections will be zero, and the term will be zero. Therefore we can ignore subtrees that consist of a single leaf, just like when computing shared (T, T By substituting (9) with (12) in (10), we can calculate Since calculation of the sums is independent of the calculation of nonshared(v, v'), these calculations can be done serially as shown in the algorithm below.
we can now express (13) as S can be calculated in time O(id v )and using the precomputed S, (14) can be also calculated in time O(id v ). Summing the results of (14) for all nodes in T gives the number of quartets in the tree, shared(T, T). The total time usage is

Calculating the shared leaf set sizes
The algorithms presented above all rely on O(1) time access to the size of the shared leaf set |F ∩ G| for any pair of subtrees F of T and G of T', where F and G each has size at least two, i.e. contains more than a single leaf. We will refer to |F ∩ G| as the intersection size of subtree F and G. In [9] and O(n 2 ) time and space algorithm is presented for computing the intersection sizes of all pairs of subtrees F and G of two binary input trees. A straightforward generalization of this algorithm to two input trees T and T' of arbitrary degrees results in an O(n 2 dd') time and O(n 2 ) space algorithm, which gives a worst case running time of O(n 4 ).
In this section we will present an improved algorithm for computing the intersection sizes of all pairs of subtrees of T and T' which runs in time O(n + |V||V'|) and space O(|V||V'|). We will assume that the size of each subtree F of T and G of T', i.e. |F| and |G|, is available in time O (1). This can be achieved as presented in the next section. Our algorithm for computing all intersection sizes is as follows. Choose an arbitrary node r in T and an arbitrary node r' in T'. Rooting the trees in r and r' respectively gives rise to two rooted trees T r and , Fig. 10 shows an example of rooting a tree. Calculating the shared leaf set sizes of T r and and all subtrees in both trees can be done using: where F i are all subtrees of T r and G j are all subtrees of .
This can be calculated using dynamic programming in time O(n 2 ): Except |T r ∩ |, the shared leafset sizes calculated by (15) are the shared leaf set sizes of all rooted subtrees of T and T' that do not contain the nodes r and r'. Assuming that the subtree F of T does not contain r, then the subtree does contain r and similarly for r' and subtrees G and of T'. The shared leaf set sizes of these trees that do contain r and r' can be calculated from the intersection sizes that we have available using (16): An arbitrary node, r, is chosen as the root in the tree, leading to a rooted tree Figure 10 An arbitrary node, r, is chosen as the root in the tree, leading to a rooted tree.
In other words all shared leaf set sizes can be calculated in time O(n 2 ). First the shared leaf set sizes for subtrees that to not contain an arbitrary node r and r' from each tree are can be defined recursively as: The first term counts all leaves directly connected to both r and r'. The second term counts all leaves connected directly to r', that are also in T r , but not directly connected to r. The third term counts all leaves connected directly to r, that are also in , but not directly connected to r'. Summing these three terms counts all leaves present in both subtrees, but leaves not connected directly to the roots are counted twice, and are subtracted by the last term.
Since (17) is only summing over the non-leaf children of a given internal node, calculating the shared leaf set sizes of all these pairs of subtrees can be done using dynamic programming in time: By the same arguments as above, the rest of the shared leaf set sizes can be computed in time O(|V||V'|) and space O(|V||V'|). Therefore the total running time of the algorithm is O(n + |V||V'|) and space usage is O(|V||V'|).
Calculating the sizes of all subtree leaf sets All of the above algorithms make use of the sizes of the leaf sets of the rooted subtrees of the input trees, either directly or indirectly. Rooting T in an arbitrary node r gives rise to the rooted tree T r . Every subtree F x of T r is a rooted subtree of T, and x is also a rooted subtree of T.
Note that the set of subtrees F x ∪ x , since one tree is the complement of the other, contains all subtrees of T and that | x | = n -|F x |. By using dynamic programming the sizes of all subtrees, F x , can be computed by a single traversal of T r . For each F x the size of x can be computed in constant time, since n is known. This means that all leaf set sizes of a tree of arbitrary degree can be calculated in time O(n).