Efficient error correction algorithms for gene tree reconciliation based on duplication, duplication and loss, and deep coalescence

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 n2 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 large-scale 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 large-scale 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 (GT-ST) 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 GT-ST 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., [2][3][4]). 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 GT-ST incongruence and can radically affect GT-ST 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 2-3 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., [8][9][10]).
Several approaches have been proposed to address gene tree error in GT-ST reconciliation. First, questionable nodes in a gene tree or nodes with low support may be collapsed prior to gene tree reconciliation, and the resulting non-binary gene trees may be reconciled with species trees [11][12][13]. Similarly, GT-ST 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 large-scale 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 GT-ST 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 error-correction protocol that is computationally feasible for large-scale 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 large-scale 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 non-empty 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 non-root 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 leaf-mapping L 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(L 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 leaf-mapping L 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 leaf-mapping L 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 L G,S . The reconciliation cost are based on the models of gene duplication [22,23], duplication-loss [21], and deep coalescence [21].
• The duplication cost from g V(G) to S, Definition 4 (Duplication-loss cost).
• The loss cost from g V(G) to S, • The duplication-loss 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 error-correction problems
Here we give definitions for tree rearrangement operations TBR [19] and SPR [19,20], and then formulate the Error-Correction 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. Re-connect 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: 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).) We define the following neighborhoods for the SPR operation: We now state the SPR based error-correction problems for duplication (D), duplication-loss (DL), and deep coalescence (DC). Let Γ {D, DL, DC}.

Problem 1 (SPR based error-correction for Γ (SEC-Γ))
Instance: A gene tree G and a species tree S.
Find: A gene tree G* SPR G such that  The TBR based error-correction 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 error-correction problems, for duplication (D), duplication-loss (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 error-correction for the Γ (R-SEC-Γ) problem.

Problem 2 (Restricted SPR based error-correction for (R-SEC-Γ))
Instance: A gene tree G, a species tree S, and v V(G). Find Given a gene tree G and a species tree S, the SEC-Γ problem can be solved as follows: (i) solve the R-SEC-Γ 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 R-SEC-Γ problem can be solved in Θ(n 2 ) time by computing the cost C G , S 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 R-SEC-Γ problem that improves by a factor of n on the naïve solution. This speedup is achieved by semi-ordering the trees in SPR G (v), for each v V(G), such that the score-difference 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 NNI-distance, 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 .
Lemma 1. X is a tree. Proof. We prove it by showing that there exists a unique path between every two vertices in X . Let G', G'' V( X ), 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. X 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 X are of degree 1 or 3. Let G' V( X ), 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 G = NNI G 1 (v)). Hence {G 1 , G'} E( X ). Since |Ch G (y)| = 2, G' must be a degree 2 vertex in X .
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 X .
This completes the proof.
The score difference of consecutive trees in X To solve the R-SEC-Γ problems we traverse tree X . Two adjacent trees in V( X ) are one NNI operation apart. We show that C 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 X . 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 G r of Figure 2 Proof. From NNI operation, v, z' Ch G'' (x) and z, x Ch G'' (y). Also, Thus, M G ,S (x) = LCA(L G ,S (Le(G x ))) = LCA(L G ,S (Le(G y ))) = M G ,S y . □ Lemma 3. ℳ G'',S (w) = ℳ G',S (w), for all w V(G')\{x, y}.

The algorithm
We describe a general algorithm, called Algo-R-SEC-Γ, to solve the R-SEC-Γ problem for each Γ {D, DL, DC}. Initially Algo-R-SEC-Γ computes (i) the root vertex of the NNI-adjacency graph X , which we call G , by regrafting the subtree G v above the root of G, (ii) the LCA mapping from G to S, and (iii) the Γ score from G to S. Then recursively Algo-R-SEC-Γ computes the LCA mapping and Γ score for every vertex G' in X when the LCA mapping and Γ score of its parent vertex in X is known. Algorithm 1 details Algo-R-SEC-Γ. .
Proof. Lemma 1-5 and Theorem 1-3 directly imply that in order to prove the correctness of algorithm Algo-R-SEC-Γ, it is sufficient to prove that it correctly returns G' of minimum Ḡ ,G among all G' V( X ). We will show that algorithm Algo-R-SEC-accounts each G' V( X ), correctly computes Ḡ ,G for Γ {D, DL, DC}, and returns the right G' as output.
From Definition 10, V( X ) = SPR G (v). In Algo-R-SEC-Γ, step 1 prunes subtree G v and regrafts it above the root of G to create G .
Step 5 sets G' to G . The forloop in step 6 traverses subtreeḠ Sb(v) in preorder. For each traversed vertex k = Ro(Ḡ Sb(v) ), step 9 builds the tree G'':= SPR G (v, k) by applying NNI operation on the last build G''. Each for-loop iteration sets G' to the last build G'' in step 23. G and G''s constitute all the trees in SPR G (v).
For G , step 2 computes the LCA mapping and step 5 sets Ḡ ,G to zero. Following Lemma 2-4 and Theorem 2, step 10 and 11 update the LCA of G'' and step 12 computes {G ,G } by calling algorithm Algo-G-Score. Depending on Γ {D, DL, DC}, there are three cases: Case 1: Γ is D. Algo-G-Score 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. Algo-G-Score computes losses by applying the formula of Definition 4. Further, it adds 1 if there is a duplication.
Case 3: Γ is DC. Algo-G-Score, 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 Algo-R-SEC-Γ, step 13 computes G,G by adding G,G and {G ,G } . When backtracking, steps 17-22 are executed to restore the right G' to compute the next unique G'' Ch X (G). This ensures that the correct Ḡ ,G is computed for each G' V( X ). In Algo-R-SEC-Γ, step 4 sets G as the BestTree and Ḡ ,G = 0 as BestScore. Every time a new G'' Ch X (G) is encountered, step 14 compares Ḡ , G with BestScore, and updates BestTree with G'' of the minimum Ḡ , G . After the for-loop, step 24 returns the BestTree. □ Lemma 7. The R-SEC-Γ and SEC-Γ problems can be solved in Θ(n) and Θ(n 2 ) time, respectively.
Proof. We will prove that the algorithm Algo-R-SEC-Γsolves the restricted SPR based error-correction problems for each Γ {D, DL, DC} in Θ(n) time. In Algo-R-SEC-Γ, step 1 takes constant time.
Step 2 precomputes LCA values for species tree in O(n) time [24], and so, finds LCA mapping from G to S in O(n) time in bottom-up manner.
Step 3 computes the duplication, duplication-loss or deep coalescence score of G and S by calling Algo-Comp-Score. In Algo-Comp-Score, step 1 and step 2 runs for O(1) and O

Solving the TEC-Γ problems
In this section we study the TBR based error-correction problems, for duplication (D), duplication-loss (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 speed-up for the SEC-Γ problems is achieved by observing that the scores Γof any re-rooted pruned subtree and its remaining pruned tree are independent of each other. We define the R-TEC-Γ problems for the TEC-Γ problems, as we defined the R-SEC-Γ problems for the SEC-Γ problems. We will show that the R-TEC-Γ 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, S, g) Also, let G 2 := TBR G (v, x, y 1 ), for y 1 V(G)\V(G v ), and y 1 ≠ y. Observe that, ∀g ∈ V(G v ), C (G , S, g) = C (G 2 , S, g). 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.
. This follows similarly. □ Lemma 8 implies that a tree in TBR G (v) with the minimum duplication, duplication-loss, 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 R-SEC problem identifies a best regraft location in Θ(n) time. This allows to obtain a tree in TBR G (v) with the minimum duplication, duplication-loss, or deep coalescence cost by evaluating only Θ(n) trees. Thus the R-TEC-Γ problem can be solved in Θ(n) time. The TEC-Γ problem can be solved by calling the solution of R-TEC-Γ 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 RAxML-VI-HPC 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.
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 re-rooting 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

Conclusion
GT-ST 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 GT-ST reconciliation must account for gene tree error; however, any useful method also must be computationally tractable for large-scale 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 2-3 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, well-supported 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.