Abstract
Background
Gene tree  species tree reconciliation problems infer the patterns and processes of gene evolution within a species tree. Gene tree parsimony approaches seek the evolutionary scenario that implies the fewest gene duplications, duplications and losses, or deep coalescence (incomplete lineage sorting) events needed to reconcile a gene tree and a species tree. While a gene tree parsimony approach can be informative about genome evolution and phylogenetics, error in gene trees can profoundly bias the results.
Results
We introduce efficient algorithms that rapidly search local Subtree Prune and Regraft (SPR) or Tree Bisection and Reconnection (TBR) neighborhoods of a given gene tree to identify a topology that implies the fewest duplications, duplication and losses, or deep coalescence events. These algorithms improve on the current solutions by a factor of n for searching SPR neighborhoods and n^{2 }for searching TBR neighborhoods, where n is the number of taxa in the given gene tree. They provide a fast error correction protocol for ameliorating the effects of gene tree error by allowing small rearrangements in the topology to improve the reconciliation cost. We also demonstrate a simple protocol to use the gene rearrangement algorithm to improve gene tree parsimony phylogenetic analyses.
Conclusions
The new gene tree rearrangement algorithms provide a fast method to address gene tree error. They do not make assumptions about the underlying processes of genome evolution, and they are amenable to analyses of largescale genomic data sets. These algorithms are also easily incorporated into gene tree parsimony phylogenetic analyses, potentially producing more credible estimates of reconciliation cost.
Introduction
The availability of largescale genomic data from a wide variety of taxa has revealed much incongruence between gene trees and the phylogeny of the species in which the genes evolve. This incongruence may be caused by evolutionary processes such as gene duplication and loss, deep coalescence, or lateral gene transfer. The variation in gene tree topologies can be used to infer the processes of genome evolution. Gene tree  species tree (GTST) reconciliation methods seek to map the history of gene trees into the context of species evolution and thus potentially link processes of gene evolution to phenotypic changes and diversification. Yet these methods can be confounded by error in the gene trees, which also may cause incongruence between the gene and species topologies. We introduce efficient algorithms to correct gene tree topologies based on the gene duplication, duplication and loss, or deep coalescence cost models. The algorithms work by identifying the small rearrangements in the gene trees that reduce the reconciliation cost. They are extremely fast and thus amenable to analyses of enormous genomic data sets.
Perhaps the most commonly used and computationally feasible approach to GTST reconciliation is gene tree parsimony, which seeks to infer the fewest evolutionary events (e.g., duplication, loss, coalescence, or lateral gene transfer) needed to reconcile a gene tree and species tree topology [1]. This approach also can be extended to infer species phylogenies, finding the species tree that implies the fewest evolutionary events implied by the gene trees (e.g., [24]). However, the gene trees often are estimated using heuristic methods from short sequence alignments, and consequently, there is often much error in the estimated gene tree topologies. Error in the gene trees creates more GTST incongruence and can radically affect GTST reconciliation analyses, implying far more duplications, duplications and losses, or deep coalescence events than actually exist. For example, Rasmussen and Kellis [5] estimated that error in gene tree reconstruction can lead to 23 fold overestimates of gene duplications and losses. Gene tree error also can erroneously imply large numbers of duplications near the root of the species tree [6,7], and it can mislead gene tree parsimony phylogenetic analyses (e.g., [810]).
Several approaches have been proposed to address gene tree error in GTST reconciliation. First, questionable nodes in a gene tree or nodes with low support may be collapsed prior to gene tree reconciliation, and the resulting nonbinary gene trees may be reconciled with species trees [1113]. Similarly, GTST reconciliations can use a distribution of gene tree topologies, such as bootstrap gene trees, rather than a single gene tree estimate [6,14,15]. Both of these approaches may help account for stochastic error and uncertainty in gene tree topologies, but they do not explicitly confront gene tree error. Methods also exist to simultaneously infer the gene tree topology and the gene tree reconciliation with a known species tree [5,16]. While these sophisticated statistical approaches appear very promising, they are computationally intensive, and it is unclear if they will be tractable for largescale analyses. Another, perhaps a more computationally feasible, approach is to allow a limited number of local rearrangements in the gene tree topology if they reduced the reconciliation cost [17,18].
Previously [17,18] described a method to allow NNIbranch swaps on selected branches of a gene tree to reduce the reconciliation cost. Following [17,18], we address gene tree error in the reconciliation process by assuming that the correct gene tree can be found in a particular neighborhood of the given gene tree. We describe this approach for the gene duplication, duplication and loss, and deep coalescence models, which identify the fewest respective events implied from a given gene tree and given species tree. This neighborhood consists of all trees that are within one edit operation of the gene tree. While [17,18] use Nearest Neighbor Interchange (NNI) edit operations to define the neighborhood, we use the standard tree edit operations SPR [19,20] and TBR [19], which significantly extend upon the search space of the NNI neighborhood. The SPR and TBR local search problems find a tree in the SPR and TBR neighborhood of a given gene tree, respectively, that has the smallest reconciliation cost when reconciled with a given species tree. Using the algorithm from Zhang [21] the best known (naïve) runtimes are O(n^{3}) for the SPR local search problem and O(n^{4}) for the TBR local search problem, where n is the number of taxa in the given gene tree. These runtimes typically are prohibitively long for the computation of larger GTST reconciliations. We improve on these solutions by a factor of n for the SPR local search problem and a factor of n^{2 }for the TBR local search problem. This makes the local search under the TBR edit operation as efficient as under the SPR edit operation, and it provides a highspeed gene tree errorcorrection protocol that is computationally feasible for largescale genomic data sets.
We also evaluated the performance of our algorithms using the implementation of SPR based local search algorithms. Note, that the SPR neighborhood is properly contained in the TBR neighborhood for any given tree. Thus the performance of the SPR based program provides a conservative estimate of the performance of the TBR based program. We test our programs on a collection of 106 yeast gene trees, some of which contain hundreds of leaves, and we demonstrate how it can be easily incorporated into largescale gene tree parsimony phylogenetic analyses.
Basic notation and preliminaries
Throughout this paper, the term tree refers to a rooted full binary (phylogenetic) tree.
Let T be a tree. The leaf set of T is denoted by Le(T). The set of all vertices of T is denoted by V(T) and the set of all edges by E(T). The root of T is denoted by Ro(T). The set of internal vertices of T is I(T):= V(T)\Le(T).
Given a vertex v ∈ V(T), we denote the parent of v by Pa_{T}(v). Let u := Pa_{T}(v). The edge that connects v with u is written as (u, v). The first element in the pair is always the parent of the second element. The set of all children of v is denoted by Ch_{T}(v) and the children are called siblings. For w ∈ Ch_{T}(v), the sibling of w is denoted by Sb_{T}(w).
We define ≤_{T }to be the partial order on V(T) where x≤_{T }y if y is a vertex on the path between Ro(T) to x, and write x <_{T }y if x ≤ _{T }y and x ≠ y. The least common ancestor of a nonempty subset L ⊆ V(T), denoted as LCA_{T}(L), is the unique smallest upper bound of L under ≤_{T}. Given x, y ∈ V(T), d_{T}(x, y) denotes the number of edges on the unique path between x and y in T.
Given U ⊆ V(T), we denote by T(U) the unique rooted subtree of T that spans U with the minimum number of vertices. Furthermore, the restriction of T to U, denoted by T_{U}, is the rooted tree that is obtained from T(U) by suppressing all nonroot vertices of degree two. The subtree of T rooted at u ∈ V(T), denoted by T_{u}, is defined to be T_{U}, for U:= {v ∈ Le(T): v ≤_{T }u}. Two trees T_{1 }and T_{2 }are called isomorphic if there exists a bijection between the vertex sets of T_{1 }and T_{2 }which maps a vertex u_{1 }of T_{1 }to vertex u_{2 }of T_{2 }if the subtree rooted at u_{1 }in T_{1 }has the same leaf set as the subtree rooted at u_{2 }in T_{2}. If an isomorphism exists between T_{1 }and T_{2}, we write T_{1 }≃ T_{2}.
Given function f : A → B, we denote by f(A') where A' ⊆ A a set of images of each element in A' under f.
The reconciliation cost models
A species tree is a tree that depicts the branching pattern representing the divergence of a set of species, whereas a gene tree is a tree that depicts the evolutionary history among the sequences encoding one gene (or gene family) for a given set of species. We assume that each leaf of the gene tree is labeled with the species from which that gene was sampled. Let G be a gene tree and S a species tree. In order to compare G with S, we require a mapping from each gene g ∈ V(G) to the most recent species in S that could have contained g.
Definition 1 (Mapping). The leafmapping ℒ_{G,S }: Le(G) → Le(S) is a surjection that maps each leaf g∈ Le(G) to that unique leaf s ∈ Le(S) which has the same label as g. The extension ℳ_{G,S }: V(G) → V(S) is the mapping defined by ℳ_{G,S}(g):= LCA(ℒ_{G,S}(Le(G_{g}))). For convenience, we write ℳ(g) instead of ℳ_{G,S}(g) when G and S are clear from the context.
Definition 2 (Comparability). Given trees G and S, we say that G is comparable to S if a leafmapping ℒ_{G,S}(g) is well defined.
Throughout this paper we use the following terminology: G is a gene tree that is comparable to the species tree S through a leafmapping ℒ_{G,S}, and n is the number of taxa present in both input trees.
Now we define different reconciliation costs from G to S for a given mapping ℳ_{G,S }that extends ℒ_{G,S}. The reconciliation cost are based on the models of gene duplication [22,23], duplicationloss [21], and deep coalescence [21].
Definition 3 (Duplication cost).
• The duplication cost from g ∈ V(G) to S,
• The duplication cost from G to S, .
Definition 4 (Duplicationloss cost).
• The loss cost from g ∈ V(G) to S,
• The duplicationloss cost from G to S, .
Definition 5 (Deep coalescence cost).
• The number of lineages from g ∈ V(G) to h ∈ Ch(g) in S,
• The deep coalescence cost from G to S, .
The errorcorrection problems
Here we give definitions for tree rearrangement operations TBR [19] and SPR [19,20], and then formulate the ErrorCorrection problems that were motivated in the introduction.
Definition 6 (Tree Bisection and Reconnection (TBR)). Let T be a tree. For this definition, we regard the planted tree Pl(T) as the tree obtained from adding the root edge{r, Ro(T)} to E(T), where r ∉ V(T).
Let e := (u, v) ∈ E(T), and X and Y be the connected components that are obtained by removing edge e from T such that v ∈ X and u ∈ Y. We define TBR_{T }(v, x, y) for x ∈ X and y ∈ Y to be the tree that is obtained from Pl(T) by first deleting edge e, and then adjoining a new edge f between X and Y as follows:
1. If x ≠ Ro(X) then suppress Ro(X) and create a new root by subdividing edge (Pa(x), x).
2. Subdivide edges (Pa(y), y) by introducing a new vertex y'.
3. Reconnect components X and Y by adding edge f = (y', Ro(x).
4. Suppress the vertex u, and rename vertex y' as u.
5. Contract the root edge.
We say that, the tree TBR_{T }(v, x, y) is obtained from T by a tree bisection and reconnection (TBR) operation that bisects the tree T into the components X and Y , and reconnects them above the nodes x and y. (See Figure 1.) We define the following neighborhoods for the TBR operation:
1. TBR_{G}(v, x) := ∪ _{y∈Y }TBR_{G}(v, x, y)
2. TBR_{G}(v) := ∪_{x∈X }TBR_{G}(v, x)
3. TBR_{G }:= ∪_{(u, v)∈E(G) }TBR_{G}(v)
Figure 1. An TBR operation. Tree T' = TBR_{T}(v, x, y) results from T after performing single TBR operation.
Definition 7 (Subtree Prune and Regrafting (SPR)). The SPR operation is defined as a special case of the TBR operation. Let e := (u, v) ∈ E(T), and X and Y be the connected components that are obtained by removing edge e from T such that v ∈ X and u ∈ Y. We define SPR_{T }(v, y) for y ∈ Y to be TBR_{T }(v, v, y). We say that the tree SPR_{T }(v, y) is obtained from T by performing subtree prune and regraft (SPR) operation that prunes subtree T_{v }and regrafts it above y. (See Figure 2(a).)
Figure 2. The NNI adjacency graph. (a) The tree is obtained from G by pruning and regrafting subtree G_{v }to the root of G. The vertex x ∈ V(G) is suppressed, and the new vertex above root in is named x. (b) Two NNI operations NNI_{G'}(z') and NNI_{G'}(z) produce leftchild G'_{l }and rightchild G'_{r }of G' in .
We define the following neighborhoods for the SPR operation:
1. SPR_{G}(v) := ∪_{y∈Y }SPR_{G}(v, y)
2. SPR_{G }:= ∪_{(u, v)∈E(G) }SPR_{G}(v)
We now state the SPR based errorcorrection problems for duplication (D), duplicationloss (DL), and deep coalescence (DC). Let Γ ∈ {D, DL, DC}.
Problem 1 (SPR based errorcorrection for Γ (SECΓ))
Instance: A gene tree G and a species tree S.
Find: A gene tree G* ∈ SPR_{G }such that .
The TBR based errorcorrection for Γ (TECΓ) problems are defined analogously to the SPR based errorcorrection for Γ (SECΓ) problems.
Solving the SECΓ problems
In this section we study the SPR based errorcorrection problems, for duplication (D), duplicationloss (DL), and deep coalescence (DC), in more detail. Our efficient solution for these problems are based on solving restricted versions of these problems efficiently. For each Γ ∈ {D, DL, DC} we first define a restricted version of the SECΓ problem, which we call the restricted SPR based errorcorrection for the Γ (RSECΓ) problem.
Problem 2 (Restricted SPR based errorcorrection for (RSECΓ))
Instance: A gene tree G, a species tree S, and v ∈ V(G).
Find: A gene tree G* ∈ SPR_{G}(v) such that .
Observation 1. Let Γ ∈ {D, DL, DC}. Given a gene tree G and a species tree S, the SECΓ problem can be solved as follows: (i) solve the RSECΓ problem for every v ∈ V(G) where v ≠ Ro(G), (ii) under all solutions found return a minimum scoring gene tree G*.
Naïvely, the RSECΓ problem can be solved in Θ(n^{2}) time by computing the cost for each G' ∈ SPR_{G}(v). The cost for a given gene and species tree can be computed in Θ(n) time [21]. We introduce a novel algorithm for the RSECΓ problem that improves by a factor of n on the naïve solution. This speedup is achieved by semiordering the trees in SPR_{G}(v), for each v ∈ V(G), such that the scoredifference of any two consecutive trees in this order can be computed in constant time.
Ordering the trees in SPR_{G}(v)
Consider a graph on trees in SPR_{G}(v), in which every two adjacent trees are one NNI [19] operation apart. We show that such a graph is a rooted full binary tree, after providing necessary definitions.
Definition 8 (Nearest Neighbor Interchange (NNI)). We define the NNI operation as a special case of the SPR operation. Let e ∈ E(T) where e := (u, v), and X and Y be the connected components that are obtained by removing edge e from T such that v ∈ X and u ∈ Y. We define NNI_{T }(v) to be SPR_{T }(v, y) for y:= Pa(u), and say that NNI_{T }(v) is obtained from T by performing nearest neighbor interchange (NNI) operation that prunes subtree T_{v }and regrafts it above the parent of v's parent. (See Figure 2(b).)
Definition 9 (NNI distance). Let the NNIdistance, denoted as d_{NNI}(T_{1, }T_{2}), between two trees T_{1 }and T_{2 }over n taxa be the minimum number of NNI operations required to transform T_{1 }into T_{2}.
Definition 10 (NNIadjacency graph). The NNIadjacency graph, denoted as , is the graph where V = SPR_{G}(v) and {T_{1, }T_{2}} ∈ E if d_{NNI}(T_{1, }T_{2}) = 1.
Proof. We prove it by showing that there exists a unique path between every two vertices in . Let G', G'' ∈ V(), thus G', G'' ∈ SPR_{G}(v). Let G':= SPR_{G}(v, x_{1}) and G'':= SPR_{G}(v, x_{2}). We use induction on d_{G}(x_{1}, x_{2}). Let d_{G}(x_{1}, x_{2}) = 1 and assume without loss of generality that x_{2 }= Pa_{G}(x_{1}). Thus, G'' = NNI_{G'' }(Sb(x_{1})). So the hypothesis holds for d_{G}(x_{1}, x_{2}) = 1. Assume now that the hypothesis is true for d_{G}(x_{1}, x_{2}) ≤ k and suppose d_{G}(x_{1}, x_{2}) = k + 1. Since G is a tree, there must be a unique path between x_{1 }and x_{2}; let y be a vertex on this path. Let d_{G}(y, x_{1}) = 1, and G^{n }:= SPR_{G}(v, y). If y = Pa_{G}(x_{1}), then G^{n }= NNI_{G'}(v); otherwise G^{n }= NNI_{G'}(Sb(y)). Since d_{G}(y, x_{2}) = k, thus (by induction hypothesis) the hypothesis is valid for d_{G}(x_{1}, x_{2}) = k + 1. □
Theorem 1. is a rooted full binary tree.
Proof. In view of Lemma 1, it suffices to show that except a unique vertex of degree 2 all other vertices in are of degree 1 or 3. Let G' ∈ V(), thus G' = SPR_{G}(v, y) for some y ∈ V(G). There are three cases:
Case 1: y is a root. Let y_{1 }∈ Ch_{G}(y). Let G^{1 }:= SPR_{G}(v, y_{1}), thus ). Hence {G^{1}, G'} ∈ E(). Since Ch_{G}(y) = 2, G' must be a degree 2 vertex in .
Case 2: y is a leaf. Let y_{1 }= Pa_{G}(y). Let G^{1 }:= SPR_{G}(v, y_{1}), thus G^{1 }= NNI_{G'}(v). Hence {G^{1}, G'} ∈ E(), and consequently, G' is a degree 1 vertex in .
Case 3: y is an internal vertex. Let y_{1 }= Pa_{G}(y) and y_{2 }∈ Ch_{G}(y). Let G^{1 }:= SPR_{G}(v, y_{1}), thus G^{1 }= NNI_{G'}(v). Let G^{2 }:= SPR_{G}(v, y_{2}), thus G' = NNI_{G}^{2}(v). Since y has one parent and two children in G, thus G' is a degree 3 vertex in .
This completes the proof.
The score difference of consecutive trees in
To solve the RSECΓ problems we traverse tree . Two adjacent trees in V() are one NNI operation apart. We show that score of a tree can be computed in constant time from the LCA computation of its adjacent tree.
Let e := (G', G'') be an edge in . Let x := Pa(v), y := Sb(v), and z, z' ∈ Ch(y) in G' (see Figure 2(b)). Without loss of generality, let G'':= NNI_{G'}(z). (Observe, G'' is similar to of Figure 2(b).)
Lemma 2. ℳ_{G'',S}(y) = ℳ_{G',S}(x).
Proof. From NNI operation, v, z' ∈ Ch_{G''}(x) and z, x ∈ Ch_{G''}(y). Also, , so . Thus, . □
Lemma 3. ℳ_{G'',S}(w) = ℳ_{G',S}(w), for all w ∈ V(G')\{x, y}.
Proof. For , since , therefore ℳ_{G',S}(g) = ℳ_{G'',S}(g). Also, except for subtree , the rest of the tree remains the same in . Thus by Lemma 2, ℳ_{G',S}(Pa_{G' }(x)) = ℳ_{G'',S}(Pa_{G''}(y)). Inductively, ℳ_{G',S}(g) = ℳ_{G'',S}(g), for all g ∈ V(G')\V(G'_{x}). □
Lemma 4. ℳ_{G''},_{S}(x) = LCA(ℳ_{G',S}(v), ℳ_{G', S}(z')).
Proof. From Lemma 3, ℳ_{G'',S}(v) = ℳ_{G',S}(v) and ℳ_{G'',S}(z') = ℳ_{G',S}(z'). Thus, ℳ_{G'',S}(x) = LCA(ℳ_{G'',S}(v),ℳ_{G'',S}(z')) = LCA(ℳ_{G',S}(v),ℳ_{G',S}(z')). □
Lemma 5. , for all g ∈ V(G'')\{x, y} and Γ ∈ {D, DL, DC}.
Proof. The gene duplication and loss status of a vertex, and the number of lineages from a vertex to its children in G' can change in G'' if its mapping or mapping of any of its children changes in ℳ_{G'',S}. From Lemma 3, and also, since ℳ_{G''},_{S}(w) = ℳ_{G',S}(w), for w ∈ Ch(Pa_{G'}(x)), must have . Thus the Lemma follows. □
Let e := (G', G'') ∈ E() and Γ ∈ {D, DL, DC}. We define with respect to the given species tree S. Observe that this score can be negative too. We study how Γ_{e }can be computed efficiently for each edge e in .
□
Definition 11. Let , and let P_{G' }be a path from to G^{' }in . For G^{'}, we define the scoredifference as .
Theorem 3. For given S, G, and v ∈ V(G), the tree G' ∈ V() is the output of a RSECΓ problem iff .
Proof. Let . We prove that G' is the output of RSECΓ problem. Since , thus G' gives the minimum normalized score over all trees in V(). Hence, G' must be the output of the RSECΓ problem. The other direction follows similarly. □
The algorithm
We describe a general algorithm, called AlgoRSECΓ, to solve the RSECΓ problem for each Γ ∈ {D, DL, DC}. Initially AlgoRSECΓ computes (i) the root vertex of the NNIadjacency graph , which we call , by regrafting the subtree G_{v }above the root of G, (ii) the LCA mapping from to S, and (iii) the Γ score from to S. Then recursively AlgoRSECΓ computes the LCA mapping and Γ score for every vertex G' in when the LCA mapping and Γ score of its parent vertex in is known. Algorithm 1 details AlgoRSECΓ.
Algorithm 1  AlgoRSECΓ
Input: A gene tree G, a species tree S, and v ∈ V(G)
Output: A tree G* ∈ SPR_{G}(v) such that
01. Compute by pruning G_{v }and regrafting at Ro(G)
06. For each in preorder traversal of , do
07. If not backtracking, then
08. Set x = Pa_{G'}(v), y = Sb_{G}'(v)
09. Set G'' = NNI_{G'}(Sb_{G'}(k))
10. Set ℳ_{G''.S }= ℳ_{G',S}, ℳG'',_{S}(y) = ℳ_{G',S}(x)
11. ℳ_{G'',S}(x) = LCA(ℳ_{G',S}(k),ℳ_{G',S}(v))
16. Else,
18. Set G'' = NNI_{G'}(v)
24. Return BestTree
Algorithm 2  AlgoCompScore
Input: A gene tree G, a species tree S, and LCA mapping ℳ_{G,S}
01. score = 0
02. For each g ∈ I(G) in preorder traversal of G, do
04. If Γis DC, then
05. score = score  I(S)
06. Return score
Algorithm 3  AlgoGScore
Input: A gene tree G, a species tree S, LCA mapping ℳ_{G}_{,S}, and g ∈ I(G)
01. If Γ is D, then
02. Ifℳ(g) ∈ ℳ(Ch(g)), then
03. Return 1
04. ElseIf Γis DL, then
06. Ifℳ(g) ∈ ℳ(Ch(g)), then
07. Return ls + 1
08. Else
09. Return ls
10. Else //Γ is DC
Lemma 6. The RSECΓ problem is correctly solved by AlgoRSECΓ.
Proof. Lemma 15 and Theorem 13 directly imply that in order to prove the correctness of algorithm AlgoRSECΓ, it is sufficient to prove that it correctly returns G' of minimum among all G' ∈ V(). We will show that algorithm AlgoRSEC accounts each G' ∈ V(), correctly computes for Γ ∈ {D, DL, DC}, and returns the right G' as output.
From Definition 10, V() = SPR_{G}(v). In AlgoRSECΓ, step 1 prunes subtree G_{v }and regrafts it above the root of G to create . Step 5 sets G' to . The forloop in step 6 traverses subtree in preorder. For each traversed vertex , step 9 builds the tree G'':= SPR_{G}(v, k) by applying NNI operation on the last build G''. Each forloop iteration sets G' to the last build G'' in step 23. and G''s constitute all the trees in SPR_{G}(v).
For , step 2 computes the LCA mapping and step 5 sets to zero. Following Lemma 24 and Theorem 2, step 10 and 11 update the LCA of G'' and step 12 computes by calling algorithm AlgoGScore. Depending on Γ ∈ {D, DL, DC}, there are three cases:
Case 1: Γ is D. AlgoGScore returns 1, if the vertex g ∈ V(G'') maps to the same vertex in S as any of its children maps to, otherwise 0.
Case 2: Γ is DL. AlgoGScore computes losses by applying the formula of Definition 4. Further, it adds 1 if there is a duplication.
Case 3: Γ is DC. AlgoGScore, returns the number of lineages from g to each of its children h ∈ Ch(g) in S. For each h ∈ Ch(g), depth of ℳ(g) is subtracted from depth of ℳ(h) to count number of edges between ℳ(g) and ℳ(h).
In AlgoRSECΓ, step 13 computes by adding and . When backtracking, steps 1722 are executed to restore the right G' to compute the next unique G'' ∈ Ch_{}(G). This ensures that the correct is computed for each G' ∈ V().
In AlgoRSECΓ, step 4 sets as the BestTree and as BestScore. Every time a new G'' ∈ Ch_{}(G) is encountered, step 14 compares with BestScore, and updates BestTree with G'' of the minimum . After the forloop, step 24 returns the BestTree. □
Lemma 7. The RSECΓ and SECΓ problems can be solved in Θ(n) and Θ(n^{2}) time, respectively.
Proof. We will prove that the algorithm AlgoRSECΓsolves the restricted SPR based errorcorrection problems for each Γ ∈ {D, DL, DC} in Θ(n) time. In AlgoRSECΓ, step 1 takes constant time. Step 2 precomputes LCA values for species tree in O(n) time [24], and so, finds LCA mapping from to S in O(n) time in bottomup manner. Step 3 computes the duplication, duplicationloss or deep coalescence score of and S by calling AlgoCompScore. In AlgoCompScore, step 1 and step 2 runs for O(1) and O(n) time, respectively. Step 3 calls AlgoGScore in each iteration of forloop. AlgoGScore runs for O(1) time for Γ ∈ {D, DL, DC}.
When Γ is DC, steps 4 and 5 are further executed in AlgoCompScore for constant time. Thus in AlgoRSECΓ, step 3 runs for O(n) time. Further, steps 4 and 5 take constant time. The loop of step 6 runs for Θ(n) time. If condition of step 7 is true, steps 810 executes in constant time. With precomputed LCA values from step 2, step 11 executes in constant time. AlgoGScore runs for constant time for Γ ∈ {D, DL, DC}, and lets step 12 to execute in constant time. Further, steps 1315 execute for constant time too. If the condition in step 7 is false, then steps 1722 execute in constant time, similarly. Finally, step 23 runs for constant time, and hence, the RSECΓ problem can be solved in Θ(n) time. From Observation 1, AlgoRSECΓ is called Θ(n) time to solve SECΓ problem. Thus, the SECΓ problem can be solved in Θ(n^{2}) time. □
Solving the TECΓ problems
In this section we study the TBR based errorcorrection problems, for duplication (D), duplicationloss (DL), and deep coalescence (DC). More precisely, we extend our solution for the SECΓ problems to solve the TECΓ problems. A TBR operation can be viewed as an SPR operation, except that the pruned subtree can be rerooted before it is regrafted. Our speedup for the SECΓ problems is achieved by observing that the scores Γof any rerooted pruned subtree and its remaining pruned tree are independent of each other. We define the RTECΓ problems for the TECΓ problems, as we defined the RSECΓ problems for the SECΓ problems. We will show that the RTECΓ problems can be solved by solving two smaller problems separately and combining their solutions.
Definition 12. Let T be a tree and x ∈ V(T). RR(T, x) is defined to be the tree T , if x = Ro(T) or x ∈ Ch(Ro(T)). Otherwise, RR(T, x) is the tree obtained by suppressing Ro(T), and subdividing the edge (Pa(x), x) by the new root node.
Lemma 8. Given a tuple 〈G, S, v〉, and G'':= TBR_{G}(v, x, y), for x ∈ V(G_{v}), y ∈ V(G)\V(G_{v}). Then, and .
Proof. (⇒) Let G^{1 }:= TBR_{G}(v, x_{1, }y), for x_{1 }∈ V(G_{v}), and x_{1 }≠ x. Now observe that, Also, let G^{2 }:= TBR_{G}(v, x, y_{1}), for y_{1 }∈ V(G)\V(G_{v}), and y_{1 }≠ y. Observe that, . Thus, if G'' gives the minimum duplication, duplicationloss, or deep coalescence score among all trees in TBR_{G}(v), then the score contribution of vertices in V(G_{v}) and V(G)\V(G_{v}) is independent. Now looking at vertices of G, the best score is achieved when G_{v }is rooted at x, i.e. ; also the best score is achieved when RR(G_{v, }x) is regrafted at y, i.e., . This follows similarly. □
Lemma 8 implies that a tree in TBR_{G}(v) with the minimum duplication, duplicationloss, or deep coalescence cost can be obtained by optimizing the rooting for the pruned subtree, and the regraft location, separately. A best rooting for the pruned subtree is linear time computable [17,25], and the solution to the RSEC problem identifies a best regraft location in Θ(n) time. This allows to obtain a tree in TBR_{G}(v) with the minimum duplication, duplicationloss, or deep coalescence cost by evaluating only Θ(n) trees. Thus the RTECΓ problem can be solved in Θ(n) time. The TECΓ problem can be solved by calling the solution of RTECΓ problem Θ(n) times, and Theorem 4 follows.
Theorem 4. The TECΓ problem can be solved in Θ(n^{2}) time.
Experimental results
We tested the performance of the gene tree rearrangement algorithms on a set of 106 gene alignments containing sequences from 8 yeast taxa from Rokas et al. [26]. There is a well accepted phylogeny for the yeast species, and the data set has been used to test algorithms for gene tree parsimony based on the deep coalescence problem [27,28]. We constructed maximum likelihood gene trees for each gene using RAxMLVIHPC version 7.0.4 [29], the gene trees were rooted with the outgroup Candida albicans. We used the new error correction algorithms to examine how much a single SPR rearrangement in the gene tree reduces the reconciliation cost based on deep coalescence and also gene duplications and losses. Over all genes the SPR error correction reduced the total deep coalescence cost from 151 to 53 (Table 1) and the duplication and loss cost from 481 to 175 (Table 2). Both the algorithms took only seconds to run for all 106 genes on a standard laptop.
Table 1. Error correction based on deep coalescence model
Table 2. Error correction based on duplication and loss model
We also implemented a protocol to use the gene rearrangement algorithm to correct for gene tree error in gene tree parsimony phylogenetic analyses. We first took a collection of input gene trees and performed a SPR species tree search using Duptree [30], which seeks the species tree with the minimum gene duplication cost. We used the duplication only cost (instead of duplications and losses) because when there is no complete sampling of all existing genes, the loss estimates may be inflated by missing sequences. After finding the locally optimal species tree, we used our SPR gene tree rearrangement algorithm to find gene tree topologies with a lower duplication cost. We then performed another SPR species tree search using Duptree, starting from the locally optimal species tree and using the new gene tree topologies. This search strategy is similar to rerooting protocol in Duptree, which checks for better gene tree roots after a SPR species tree search [30,31]. We tested this protocol on data set of 6,084 genes (with a combined 81,525 leaves) from 14 seed plant taxa. This is the same data set used by [31], except that all gene tree clades containing sequences from a single species were collapsed to a single leaf. Our original SPR tree search found a species tree with 23,500 duplications. The SPR tree search after the gene tree rearrangements identified the same species tree, but the new gene trees had a reconciliation cost of only 18,213. This tree search protocol took just under 4 hours on a Mac Powerbook with a 2 GHz Intel Core 2 Duo processor and 2 GB memory.
Conclusion
GTST reconciliation provides a powerful approach to study the patterns and processes of gene and genome evolution. Yet it can be thwarted by the error that is an inherent part of gene tree inference. Any reliable method for GTST reconciliation must account for gene tree error; however, any useful method also must be computationally tractable for largescale genomic data. We introduce fast and effective algorithms to correct error in the gene trees. These algorithms, based on SPR and TBR rearrangements, greatly extend upon the range of possible errors in the gene tree from existing algorithms [17,18], while remaining fast enough to use on data sets with thousands of genes. These algorithms can correct trees based on a broad variety of evolutionary factors that can cause conflict between gene trees and species trees, including gene duplication, duplications and losses, and deep coalescence.
Our analysis on 106 yeast gene trees demonstrates that even a single SPR correction on the gene trees can radically improve upon the reconciliation cost. Our results in the yeast analysis are very similar to the 23 fold improvement in implied duplications and losses reported from the parametric gene tree estimation and reconciliation method of Rasmussen and Kellis [5]. However, in contrast, to this computationally complex method, the gene tree rearrangement algorithm is extremely fast, does not require assumption about the rates of duplication and loss, and is also amenable to corrections based on deep coalescence and duplications as well as duplications and losses. We do not claim that the gene correction algorithms produce a more accurate reconciliation than these parametric methods. However, they do present an extremely fast and flexible alternative.
We also demonstrated that this error correction protocol could easily be incorporated into a gene tree parsimony phylogenetic analysis. Previous studies have emphasized that gene tree parsimony is sensitive to the topology of the input trees. For example, the species tree may differ whether the gene trees are made using parsimony or maximum likelihood [8,10]. In our study, although the gene tree rearrangement did not affect the species tree inference, it did greatly reduce the gene duplication reconciliation cost.
While the results of the experiments are promising, they also suggest several directions for future research. First, further investigation is needed to characterize the effects of error on gene tree topologies. For example, it seems likely that gene tree errors may extend beyond a single SPR or TBR neighborhood. Yet, if we allow unlimited rearrangements, the gene trees will simply converge on the species tree topology. One simple improvement may be to weight the possible gene tree rearrangements based on support for different clades in the gene tree. Thus, wellsupported clades may be rarely or never be subject to rearrangement, while poorly supported clades may be subject to extensive rearrangements. Finally, these approaches implicitly assume that all differences between gene trees and species trees are due to either coalescence, duplications, or duplications and losses. Future work will seek to combine these objectives and also address lateral transfer.
Competing interests
The authors declare that they have no competing interests.
Author's contributions
RC was responsible for algorithm design and program implementation, and wrote major parts of the paper. JGB performed the experimental evaluation and the analysis of the results, and contributed to the writing of the paper. OE supervised the project, contributed to the algorithmic design and writing of the paper.
All authors read and approved the final manuscript.
Acknowledgements
The authors thank André Wehe for his generous support with the implementation. This work was conducted in parts with support from the Gene Tree Reconciliation Working Group at NIMBioS through NSF award EF0832858, with additional support from the University of Tennessee. R.C. and O.E. were supported in parts by NSF awards #0830012 and #10117189.
This article has been published as part of BMC Bioinformatics Volume 13 Supplement 10, 2012: "Selected articles from the 7th International Symposium on Bioinformatics Research and Applications (ISBRA'11)". The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/13/S10.
References

Maddison WP: Gene Trees in Species Trees.
Systematic Biology 1997, 46:523536. Publisher Full Text

Goodman M, Czelusniak J, Moore GW, RomeroHerrera AE, Matsuda G: Fitting the gene lineage into its species lineage. A parsimony strategy illustrated by cladograms constructed from globin sequences.
Systematic Zoology 1979, 28:132163. Publisher Full Text

Guigó R, Muchnik I, Smith TF: Reconstruction of Ancient Molecular Phylogeny.
Molecular Phylogenetics and Evolution 1996, 6(2):189213. PubMed Abstract  Publisher Full Text

Slowinski JB, Knight A, Rooney AP: Inferring Species Trees from Gene Trees: A Phylogenetic Analysis of the Elapidae (Serpentes) Based on the Amino Acid Sequences of Venom Proteins.
Molecular Phylogenetics and Evolution 1997, 8:349362. PubMed Abstract  Publisher Full Text

Rasmussen MD, Kellis M: A Bayesian approach for fast and accurate gene tree reconstruction.
Molecular Biology and Evolution 2011, 28:273290. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Burleigh JG, Bansal MS, Wehe A, Eulenstein O: Locating LargeScale Gene Duplication Events through Reconciled Trees: Implications for Identifying Ancient Polyploidy Events in Plants.
Journal of Computational Biology 2009, 16:10711083. PubMed Abstract  Publisher Full Text

Hahn MW: Bias in phylogenetic tree reconciliation methods: implications for vertebrate genome evolution.
Genome Biology 2007, 8:R141. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Burleigh JG, Bansal MS, Eulenstein O, Hartmann S, Wehe A, Vision TJ: Genomescale phylogenetics: inferring the plant tree of life from 18,896 discordant gene trees.
Systematic Biology 2011, 60(2):117125. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Huang H, Knowles LL: What Is the Danger of the Anomaly Zone for Empirical Phylogenetics?
Systematic Biology 2009, 58:527536. PubMed Abstract  Publisher Full Text

Sanderson MJ, McMahon MM: Inferring angiosperm phylogeny from EST data with widespread gene duplication.

BerglundSonnhammer A, Steffansson P, Betts MJ, Liberles DA: Optimal Gene Trees from Sequences and Species Trees Using a Soft Interpretation of Parsimony.
Journal of Molecular Evolution 2006, 63:240250. PubMed Abstract  Publisher Full Text

Vernot B, Stolzer M, Goldman A, Durand D: Reconciliation with nonbinary species trees.

Yu Y, Warnow T, Nakhleh L: Algorithms for MDCBased Multilocus Phylogeny Inference. In RECOMB, Volume 6577 of Lecture Notes in Computer Science. Edited by Bafna V, Sahinalp SC. Springer; 2011:531545.

Cotton JA, Page RDM: Going nuclear: gene family evolution and vertebrate phylogeny reconciled.
P Roy Soc Lond B Biol 2002, 269:15551561. Publisher Full Text

Joly S, Bruneau A: Measuring Branch Support in Species Trees Obtained by Gene Tree Parsimony.
Systematic Biology 2009, 58:100113. PubMed Abstract  Publisher Full Text

Arvestad L, Berglund A, Lagergren J, Sennblad B: Gene tree reconstruction and orthology analysis based on an integrated model for duplications and sequence evolution.

Chen K, Durand D, FarachColton M: Notung: a program for dating gene duplications and optimizing gene family trees.
Journal of Computational Biology 2000, 7:429447. PubMed Abstract  Publisher Full Text

Durand D, Halldórsson BV, Vernot B: A Hybrid MicroMacroevolutionary Approach to Gene Tree Reconstruction.
Journal of Computational Biology 2006, 13(2):320335. PubMed Abstract  Publisher Full Text

Allen BL, Steel M: Subtree transfer operations and their induced metrics on evolutionary trees.
Annals of Combinatorics 2001, 5:113. Publisher Full Text

Bordewich M, Semple C: On the computational complexity of the rooted subtree prune and regraft distance.

Zhang L: On a MirkinMuchnikSmith Conjecture for Comparing Molecular Phylogenies.
Journal of Computational Biology 1997, 4(2):177187. PubMed Abstract  Publisher Full Text

Page RDM: Maps between trees and cladistic analysis of historical associations among genes, organisms, and areas.

Eulenstein O: Predictions of geneduplications and their phylogenetic development.
PhD thesis.
University of Bonn, Germany 1998. [GMD Research Series No. 20/1998, ISSN: 14352699]

Górecki P, Tiuryn J: Inferring phylogeny from whole genomes.

Rokas A, Williams BL, King N, Carroll SB: Genomescale approaches to resolving incongruence in molecular phylogenies.
Nature 2003, 425(6960):798804. PubMed Abstract  Publisher Full Text

Than C, Nakhleh L: Species tree inference by minimizing deep coalescences.
PLoS Comput Biol 2009, 5(9):e1000501. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bansal MS, Burleigh JG, Eulenstein O: Efficient genomescale phylogenetic analysis under the duplicationloss and deep coalescence cost models.
BMC Bioinformatics 2010, 11(Suppl 1):S42. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Stamatakis A: RAxMLVIHPC: maximum likelihoodbased phylogenetic analyses with thousands of taxa and mixed models.
Bioinformatics 2006, 22(21):26882690. PubMed Abstract  Publisher Full Text

Wehe A, Bansal MS, Burleigh JG, Eulenstein O: DupTree: a program for largescale phylogenetic analyses using gene tree parsimony.

Chang W, Burleigh JG, FernándezBaca D, Eulenstein O: An ILP solution for the gene duplication problem.
BMC Bioinformatics 2011, 12(Suppl 1):S14. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text