- Proceedings
- Open access
- Published:

# Algorithms: simultaneous error-correction and rooting for gene tree reconciliation and the gene duplication problem

*BMC Bioinformatics*
**volume 13**, Article number: S14 (2012)

## Abstract

### Background

Evolutionary methods are increasingly challenged by the wealth of fast growing resources of genomic sequence information. Evolutionary events, like gene duplication, loss, and deep coalescence, account more then ever for incongruence between gene trees and the actual species tree. Gene tree reconciliation is addressing this fundamental problem by invoking the minimum number of gene duplication and losses that reconcile a rooted gene tree with a rooted species tree. However, the reconciliation process is highly sensitive to topological error or wrong rooting of the gene tree, a condition that is not met by most gene trees in practice. Thus, despite the promises of gene tree reconciliation, its applicability in practice is severely limited.

### Results

We introduce the problem of reconciling unrooted and erroneous gene trees by simultaneously rooting and error-correcting them, and describe an efficient algorithm for this problem. Moreover, we introduce an error-corrected version of the gene duplication problem, a standard application of gene tree reconciliation. We introduce an effective heuristic for our error-corrected version of the gene duplication problem, given that the original version of this problem is NP-hard. Our experimental results suggest that our error-correcting approaches for unrooted input trees can significantly improve on the accuracy of gene tree reconciliation, and the species tree inference under the gene duplication problem. Furthermore, the efficiency of our algorithm for error-correcting reconciliation is capable of handling truly large-scale phylogenetic studies.

### Conclusions

Our presented error-correction approach is a crucial step towards making gene tree reconciliation more robust, and thus to improve on the accuracy of applications that fundamentally rely on gene tree reconciliation, like the inference of gene-duplication supertrees.

## Background

The wealth of newly sequenced genomes has provided us with an unprecedented resource of information for phylogenetic studies that will have extensive implications for a host of issues in biology, ecology, and medicine, and promise even more. Yet, before such phylogenies can be reliably inferred, challenging problems that came along with the newly sequenced genomes have to be overcome. Evolutionary biologists have long realized that gene-duplication and subsequent loss, a fundamental evolutionary process [1], can largely obfuscate phylogenetic inference [2]. Gene-duplication can form complex evolutionary histories of genes, called gene trees, whose topologies are traditionally used to derive species trees. This approach relies on the assumption that the topologies from gene trees are consistent with the topology of the species tree. However, frequently genes that evolve from different copies of ancestral gene-duplications can become extinct and result in gene trees with correct topologies that are inconsistent with the topology of the actual species tree (see Figure 1). In many such cases phylogenetic information from the gene trees is indispensable and may still be recovered using gene tree reconciliation.

### Related work

Gene tree reconciliation is a well-studied method for resolving topological inconsistencies between a gene tree and a trusted species tree [2–7]. Inconsistencies are resolved by invoking gene-duplication and loss events that reconcile the gene tree to be consistent with the actual species tree. Such events do not only reconcile gene trees, but also lay foundation for a variety of evolutionary applications including ortholog/paralog annotation of genes, locating episodes of gene-duplications in species trees [8–10], reconstructing domain decompositions [11], and species supertree construction [8, 12–14].

A major problem in the application of gene tree reconciliation is its high sensitivity to error-prone gene trees. Even seemingly insignificant errors can largely mislead the reconciliation process and, typically undetected, infer incorrect phylogenies (e.g., [7, 15]). Errors in gene trees are often topological errors and rooting errors. Topological error results in an incorrect topology of the gene tree that can be caused by the inference process (e.g. noise in the underlying sequence data) or the inference method itself (e.g. heuristic results). This problem has been addressed for rooted gene trees by 'correcting the error'; that is, editing the given tree such that the number of invoked gene-duplications and losses is minimized [16, 17]. However, most inference methods used in practice return only unrooted gene trees (e.g. parsimony and maximum likelihood based methods) that have to be rooted for the gene tree reconciliation process. Rooting error is a wrongly chosen root in an unrooted gene tree. Whereas rooting can be typically achieved in species trees by outgroup analysis, this approach may not be possible for gene trees if there is a history of gene duplication and loss [7]. Other rooting approaches like midpoint rooting or molecular clock rooting assume a constant rate of evolution that is often unrealistic. However, rooting problems can be bypassed by identifying roots that minimize the invoked number of gene duplications and losses [7, 16–19].

In summary, even small topological error or a slightly misplaced root can incorrectly identify enormous numbers of gene duplications and losses, and therefore largely mislead the reconciliation process. Therefore, gene tree reconciliation requires gene trees that are free of error and correctly rooted at the same time [5]. However, as previous work has incorporated topological error-correction only separately from correctly rooting gene trees into the reconciliation process [16, 18], this process can still be misled.

### Our contribution

We address the problem of reconciling erroneous and unrooted gene trees by error-correcting and rooting them at the same time. Solving this problem efficiently is a crucial step towards making gene tree reconciliation more robust, and thus to improve on the accuracy of applications that rely on gene tree reconciliation like the construction of gene-duplication supertrees. We introduce the problem and design an efficient algorithm that facilitates a much more precise gene tree reconciliation, even for large-scale data sets. Our algorithm detects and corrects errors in unrooted gene trees, and thus we avoid the biologists' difficulty and uncertainty of handling erroneous gene trees and correctly rooting them. The presented experimental results suggest that our novel reconciliation algorithms can identify and correct topological error in unrooted input gene trees, and at the same time root them optimally.

Our algorithm is designed to search for the correct and rooted tree of a given unrooted tree in local search neighborhoods of the given tree. The size of these neighborhoods is described by a positive integer *k* that allows to fine-tune the search. While in theory *k* can be large it is assumed that gene trees have only small topological error, which typically can be captured by small values of *k*. For a fixed but freely choosable integer *k* the runtime of our algorithm is *O*(*l*^{k} + max(*n, m*)), where *n* and *m* is the size of the gene tree and species tree respectively, and *l* is the number of edges in the gene tree that potentially contain an error (such edges will be called *weak*). Thus, for a small error, which is expressed by *k* = 1, our algorithm runs in linear time. Our experiments show that error-correction runs of the algorithm for *k* = 3 are still possible even for trees with large number of weak edges (e.g., *l* = 200) on a standard workstation configuration.

Further, we address the problem of constructing rooted supertrees by reconciling unrooted and erroneous gene trees with assigned weak edges, a key problem in illuminating the role and effect of gene duplication and loss in shaping the evolution of organisms. We introduce the problem and develop an effective local search heuristic that makes the construction of more accurate supertrees possible and allows a much better postulation of gene duplication histories. Our experimental results demonstrate that our approach is effective in identifying gene duplication histories given erroneous gene trees and producing more accurate supertrees under gene tree reconciliation.

### Duplication-loss model

We introduce the fundamentals of the classical duplication-loss model. Our definitions are mostly adopted from [18]. For a more detailed introduction to the duplication-loss model we refer the interested reader to [2, 5, 10, 20].

Let ℐ be the set of species consisting of *N* > 0 elements. The *unrooted gene tree* is an undirected acyclic graph in which each node has degree 3 (internal nodes) or 1 (leaves), and the leaves are labeled by the elements from ℐ. A *species tree* is a rooted binary tree with *N* leaves uniquely labeled by the elements from ℐ. In some cases, a node of a tree will be referred by "cluster" of labels of its subtree leaves. For instance, a species tree (*a*, (*b, c*)) has 5 nodes denoted by: *a, b, c, bc* and *abc*. *A rooted gene tree* is a rooted binary tree with leaves labeled by the elements from ℐ. The internal nodes of a tree *T* we denote by int(*T*).

Let \mathcal{S}=\u27e8{V}_{\mathcal{S}},{E}_{\mathcal{S}}\u27e9 be a *species tree*. can be viewed as an upper semilattice with + a binary least upper bound operation and ⊤ the top element, that is, the root. In particular for *a*, b\in {V}_{\mathcal{S}}, *a* <*b* means that *a* and *b* are on the same path from the root, with *b* being closer to the root than *a*. We define the *comparability predicate D*(*a, b*) = 1, if *a* ≤ *b* or *b* ≤ *a* and *D*(*a, b*) = 0, when *a* and *b* are incomparable. The *distance function ρ*(*a, b*) is used to denote the number of edges on the unique (non-directed) path connecting *a* and *b*.

We call distinct nodes *a*, b\in {V}_{\mathcal{S}}*siblings* when *a* + *b* is a parent of *a* and *b*. For *a*, b\in {V}_{\mathcal{S}} let **Sb**(*a, b*) be the set of nodes defined by the following recurrent rule: **(i) Sb**(*a, b*) = ∅ if *a* = *b* or *a* and *b* are siblings, **(ii) Sb**(*a, b*) = {*c*} ∪ **Sb**(*a* + *c, b*), if *a* <*b* or *a* + *c* <*a* + *b*; here *c* is the sibling of *a*, and **(iii) Sb**(*a, b*) = **Sb**(*b, a*) otherwise.

By *L*(*a, b*) we denote the number of elements in **Sb**(*a, b*). Observe that *L*(*a, b*) = *ρ*(*a, b*) - 2 · (1 - *D*(*a, b*)). Let M:{V}_{\mathcal{G}}\to {V}_{\mathcal{S}} be the *least common ancestor (lca) mapping*, from rooted into that preserves the labeling of the leaves. Formally, if *v* is a leaf in then *M*(*v*) is the node in labeled by the label of *v*. If *v* is internal node in with two children *a, b*, then *M*(*v*) = *M*(*a*) + *M*(*b*). An example is depicted in Figure 1.

In this general setting let us assume that we are given a *cost function* \xi :{V}_{\mathcal{G}}\times {V}_{\mathcal{S}}\to \mathbf{R} which for all nodes v\in {V}_{\mathcal{G}}, a\in {V}_{\mathcal{S}} assigns a real *ξ*(*v, a*) representing a contribution to node *a* which comes from *v* when reconciling with . Having *ξ* we can define k\left(v\right)={\sum}_{a}\xi \left(v,a\right) to be a total contribution from *v* in the reconciliation of with . We call *κ* a *contribution* function. Finally, \sigma ={\sum}_{v}k\left(v\right) is the total cost of reconciliation of with .

Now we present examples of cost functions that are used in the duplication model. We assume that if *v* is an internal node in then *w*_{1} and *w*_{2} are its children. The *Duplication cost* function is defined as follows: *ξ*^{D}(*v, a*) = 1 if v\in \mathsf{\text{int}}\phantom{\rule{1em}{0ex}}\left(\mathcal{G}\right) and *M*(*v*) = *M*(*w*_{
i
}) = *a* for some *i*, and *ξ*^{D}(*v, a*) = 0 otherwise. The *Loss cost* function: *ξ*^{L}(*v, a*) = 1 if v\in \mathsf{\text{int}}\phantom{\rule{1em}{0ex}}\left(\mathcal{G}\right) and *a* ∈ **Sb** (*M*(*w*_{1}), *M*(*w*_{2})), and *ξ*^{L}(*v, a*) = 0 otherwise. It can be proved that if v\in \mathsf{\text{int}}\phantom{\rule{1em}{0ex}}\left(\mathcal{G}\right) then *κ*^{D}(*v*) = *D*(*M*(*w*_{1}), *M*(*w*_{2})) and *κ*^{L} (*v*) = *L*(*M*(*w*_{1}), *M*(*w*_{2})) (in both cases 0 if *v* is a leaf).

The *Duplication cost* function is defined as follows: *ξ*^{D}(*v, a*) = 1 if v\in \mathsf{\text{int}}\phantom{\rule{1em}{0ex}}\left(\mathcal{G}\right) and *M*(*v*) = *M*(*w*_{
i
}) = *a* for some *i*, and *ξ*^{D}(*v, a*) = 0 otherwise. Loss cost function: *ξ*^{L}(*v, a*) = 1 if v\in \mathsf{\text{int}}\phantom{\rule{1em}{0ex}}\left(\mathcal{G}\right) and *a* ∈ **Sb**(*M*(*w*_{1}), *M*(*w*_{2})), and *ξ*^{L}(*v, a*) = 0 otherwise. It can be proved that if v\in \mathsf{\text{int}}\phantom{\rule{1em}{0ex}}\left(\mathcal{G}\right) then *κ*^{D}(*v*) = *D*(*M*(*w*_{1}), *M*(*w*_{2})) and *κ*^{L}(*v*) = *L*(*M*(*w*_{1}), *M*(*w*_{2})) (in both cases 0 if *v* is a leaf).

Observe that a node v\in {V}_{\mathcal{G}} is called a duplication [4, 13] if *κ*^{D}(*v*) = 1. Moreover, *κ*^{L}(*v*) = *l*(*v*), where *l*(*v*) is the number of gene losses associated to *v*. It can be proved that *σ*^{D} and *σ*^{L} are the minimal number of gene duplications and gene losses (respectively) required to reconcile (or to embed) with . Please refer to [18] for more details. The example of an embedding is depicted in Figure 1.

### Introduction to unrooted reconciliation

Here we highlight some results from [18] that are used for the design of our algorithm. From now on, we assume that \mathcal{G}=\u3008{V}_{\mathcal{G}},{E}_{\mathcal{G}}\u3009 is an unrooted gene tree. We define a rooting of by selecting an edge e\in {E}_{\mathcal{G}} on which the root is to be placed. Such a rooted tree will be denoted by {\mathcal{G}}_{e}, where *v*_{*} is a new node defining the root. To distinguish between rootings of , the symbols defined in previous section for rooted gene trees will be extended by inserting index *e*. Please observe, that the mapping of the root of {\mathcal{G}}_{e} is independent of *e*. Without loss of generality the following is assumed: **(A1)** and have at least one internal node and **(A2)** *M*_{
e
}(*v*_{*})=⊤; that is, the root of every rooting is mapped into the root of (we may always consider the subtree of the species tree rooted in *M*_{
e
}(*v*_{*}) with no change of the cost).

First, we transform into a directed graph \hat{\mathcal{G}}=\u27e8{V}_{\mathcal{G}},\hat{{E}_{\mathcal{G}}}\u27e9 where \hat{{E}_{\mathcal{G}}}=\left\{\u27e8v,\phantom{\rule{2.77695pt}{0ex}}w\u27e9|\phantom{\rule{2.77695pt}{0ex}}\left\{v,\phantom{\rule{2.77695pt}{0ex}}w\right\}\in {E}_{\mathcal{G}}\right\}. In other words each edge 〈*v, w*〉 in is replaced in \hat{\mathcal{G}} by a pair of directed edges 〈*v, w*〉 and 〈*w, v*〉.

Edges in \hat{\mathcal{G}} are labeled by nodes of as follows. If v\in {V}_{\mathcal{G}} is a leaf labeled by *a*, then the edge \u27e8v,w\u27e9\in \hat{{E}_{\mathcal{G}}} is labeled by *a*. When *v* is an internal node in \hat{\mathcal{G}} we assume that 〈*w*_{1}, *v*〉 and 〈*w*_{2}, *v*〉 are labeled by *b*_{1} and *b*_{2}, respectively. Then the edge \u27e8v,{w}_{3}\u27e9\in \hat{{E}_{\mathcal{G}}}, such that *w*_{3} ≠ *w*_{1} and *w*_{3} ≠ *w*_{2} is labeled by *b*_{1} + *b*_{2}. Such labeling will be used to explore mappings of rootings of . An edge {*v, w*} in is called *asymmetric* if exactly one of the labels of 〈*v, w*〉 and 〈*w, v*〉 in \hat{\mathcal{G}} is equal to ⊤, otherwise it is called *symmetric*.

Every internal node *v*, and its neighbors in \hat{\mathcal{G}} define a subtree of \hat{{E}_{\mathcal{G}}}, called a *star* with a center *v*, as depicted in Figure 2. The edges 〈*v, w*_{
i
}〉 are called *outgoing*, while the edges 〈*w*_{
i
}*, v*〉 are called *incoming*. We will refer to the undirected edge {*v, w*_{
i
}} as *e*_{
i
}, for *i* = 1, 2, 3.

The are several types of possible star topologies based on the labeling (for proofs and details see [18]): (S1) a star has one incoming edge labeled by ⊤ and two outgoing edges labeled ⊤ and these edges are connected to the three siblings of the center, (S2) a star has exactly two outgoing edges labeled by ⊤, (S3) a star has all outgoing edges and exactly one incoming edgd labeled by ⊤, (S4) a star has all edges labelled by *top*, and (S5) a star has all outgoing edges and exactly two incoming edges labeled by ⊤. Figure 2 illustrates the star topologies.

In summary stars are basic 'puzzle-like' units that can be used to assemble them into unrooted gene trees. However, not all star compositions represent a gene tree. For instance, there is no gene tree with 3 stars of type S2. It follows from [18] (see Lemma 4) that we need the following additional condition: (C1) if a gene tree has two stars of type S2 then they share a common edge.

Now we overview the main result of [18] (see Theorem 1 for more details). Let be a species tree and be unrooted gene tree. The set of optimal edges, that is, candidates for best rootings, is defined as follows: \mathbf{M}i{\mathbf{n}}_{\mathcal{G}}=\left\{e\in {E}_{\mathcal{G}}|{\sigma}_{e}^{{M}_{\alpha ,\beta}}\mathsf{\text{isminimal}}\right\}, where {\sigma}_{e}^{{M}_{\alpha ,\beta}}is the total cost for the weighted mutation cost defined by {\xi}_{e}^{{M}_{\alpha ,\beta}}\left(v,\phantom{\rule{2.77695pt}{0ex}}a\right)=\alpha \phantom{\rule{2.77695pt}{0ex}}\cdot \phantom{\rule{2.77695pt}{0ex}}{\xi}_{e}^{D}\left(v,\phantom{\rule{2.77695pt}{0ex}}a\right)+\beta \phantom{\rule{2.77695pt}{0ex}}\cdot \phantom{\rule{2.77695pt}{0ex}}{\xi}_{e}^{L}\left(v,\phantom{\rule{2.77695pt}{0ex}}a\right), *e* is an edge in and *α, β* are two positive reals. Then **(M1)** if \left|\mathbf{M}i{\mathbf{n}}_{\mathcal{G}}\right|>1, then \mathbf{M}i{\mathbf{n}}_{\mathcal{G}} consists of all edges present in all stars of type S4 or S5, **(M2)** if \left|\mathbf{M}i{\mathbf{n}}_{\mathcal{G}}\right|=1, then \mathbf{M}i{\mathbf{n}}_{\mathcal{G}} contains exactly one symmetric edge that is present in star of type S2 or S3. From the above statements, (C1) and star topologies we can easily determine \mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}}. More precisely, the star edges outside \mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}} are asymmetric and share the same direction. Thus, to find an optimal edge it is sufficient to follow the direction of non ⊤ edges in \hat{\mathcal{G}}.

Now we summarize the time complexity of this procedure. It follows from [21] that a single lca-query (that, is *a* + *b* for nodes *a* and *b* in ) can be computed in constant time after an initial preprocessing step requiring O\left(\left|\mathcal{S}\right|\right) time. Other structures like \hat{\mathcal{G}} with the labeling can be computed in O\left(\left|\mathcal{G}\right|\right) time. The same complexity has the procedure of finding an optimal edge in . In summary an optimal edge/rooting and the minimal cost can be computed in linear time. See [18] for more details and other properties.

## Methods

First we describe our algorithm for computing the optimal cost and the set of optimal edges after one nearest neighbor interchange (NNI) operation performed on an unrooted gene tree, and then extend it to a general case with *k* NNI operations. For the definition of NNI please refer to Def. 1 and Figure 3.

### Algorithm

Now we show that a single NNI operation can be completed in constant time if all structures required for computing the optimal rootings are already constructed. First, let us assume that the following is given: (a) two positive reals *α* and *β*, a species tree , (b) lca structure for that allows to answer lca-queries in constant time, (c) an unrooted gene tree , (d) \hat{\mathcal{G}} with the labeling of edges, (e) \mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}} - the set of optimal edges, and (f) *σ* - the minimal total weighted mutation cost. As observed in the previous section (b),(d)-(f) can be computed in O\left(\text{max}\left(\left|\mathcal{S}\right|,\left|\mathcal{G}\right|\right)\right). Now we show that (c)-(f) can be computed in constant time after a single NNI operation.

*NNI operation (c) and the update of lca-mappings (d)*.

**Definition 1**. *(Single NNI operation) An NNI operation transforms a gene tree* \mathcal{G}=\left(\left({T}_{1},{T}_{2}\right),\left({T}_{3},{T}_{4}\right)\right) *into*{\mathcal{G}}^{\prime}=\left(\left({T}_{2},{T}_{3}\right),\left({T}_{1},{T}_{4}\right)\right)*, where T*_{
i
}*-s are (rooted) subtrees of* *. The edge that connects the roots of* (*T*_{1}, *T*_{2}) *and* (*T*_{3}, *T*_{4}) *in* *is denoted by e*_{0} *and called the* center edge. *For each i* = 1, 2, 3, 4 *we assume the following: w*_{
i
} *is the root of T*_{
i
}*, e*_{
i
} *is the edge connecting w*_{
i
} *with e*_{0} *and a*_{
i
} *is the lca-mapping of T*_{
i
}. *Similarly, we define the* center edge {e}_{0}^{\prime}*and* {e}_{i}^{\prime} *in*{\mathcal{G}}^{\prime}.

An NNI operation is depicted in Figure 3 with the transformation of \hat{\mathcal{G}} into \hat{{\mathcal{G}}^{\prime}}. The notation will be used from now on. Note that there is a second NNI operation, when is replaced with ((*T*_{1}, *T*_{3}), (*T*_{2}, *T*_{4})). However, it can be easily defined and therefore it is omitted here. Observe that the NNI operation (without updating of lca-mappings) can be performed in constant time for both trees.

The right part of Figure 3 depicts the transformation of \hat{\mathcal{G}}. Observe that the labels of the incoming and outgoing edges attached to each *w*_{
i
} in \hat{\mathcal{G}} do not change during this operation. Lemma 1 follows directly from this observation.

**Lemma 1**. *An NNI operation changes only the labels of the center edge*.

We conclude that updating \hat{\mathcal{G}} requires only two lca-queries, and therefore can be performed in constant time.

*Reconstruction of optimal edges (e)*. We analyze the changes of the optimal set of edges Mi{n}_{\mathcal{G}}. To this end we consider a number of cases depending on the relation between the optimal set of edges and the set of edges, incident to the nodes of the center edge. Let {C}_{\mathcal{G}}={\left\{{e}_{i}\right\}}_{i=0,\dots ,4}.

For convenience, assume that the NNI operation replaces *e*_{
i
} with {e}_{i}^{\prime} as indicated in Figure 3. We call two disjoint edges from {C}_{\mathcal{G}}*semi-alternating* if they share a common node after the NNI operation. In Figure 3 {*e*_{1}, *e*_{4}} and {*e*_{2}, *e*_{3}} are semi-alternating. For two edges *a* and *b* that are incident to the same node let ⋆(*a, b*) be the set of three edges defining the unique star that contains *a* and *b*.

**Lemma 2**. *Assuming that e*_{
i
} *is replaced by* {e}_{i}^{\prime} *after the NNI operation the set of optimal edges does not require additional changes if and only if one of the following conditions is satisfied:*

\left(\mathbf{E}Q1\right)\phantom{\rule{2.77695pt}{0ex}}\mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}}\cap {C}_{\mathcal{G}}=\varnothing ,

\left(\mathbf{E}Q2\right)\phantom{\rule{2.77695pt}{0ex}}\mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}}\supseteq {C}_{\mathcal{G}}*and each pair of semi-alternating edges contains at least one symmetric edge*,

\left(\mathbf{E}Q3\right)\phantom{\rule{2.77695pt}{0ex}}\mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}}*consists of only the center edge*,

\left(\mathbf{E}Q4\right)\phantom{\rule{2.77695pt}{0ex}}\mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}}\cap {C}_{\mathcal{G}}=\left\{{e}_{i}\right\}*for some i >* 0 *and the center is asymmetric after the NNI operation*.

Proof: (EQ1) All edges in {C}_{\mathcal{G}} are asymmetric (2 stars S1). Then, after the NNI operation {e}_{0}^{\prime} is asymmetric and ({C}_{{\mathcal{G}}^{\prime}} has 2 stars S1). (EQ2) {C}_{\mathcal{G}} consists of 2 stars of type S4/S5 and at most two asymmetric edges. It follows from EQ2 that the asymmetric edges in {C}_{{\mathcal{G}}^{\prime}} cannot form a star of type other than S5. Together with M1 it follows that {C}_{{\mathcal{G}}^{\prime}} is optimal. (EQ3) By M1 the center is symmetric in . It remains symmetric after NNI. From C1 and M2, \mathbf{M}\mathbf{i}{\mathbf{n}}_{{\mathcal{G}}^{\prime}} consists of the center edge. (EQ4) Note, that the type of \star \left({e}_{i}^{\prime},{e}_{0}^{\prime}\right) is *S* 1, *S* 2 or *S* 3.

**Lemma 3** (NE1). *If* \mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}}\supseteq {C}_{\mathcal{G}}*and there exists a pair* {*e*_{
i
}*, e*_{
j
}} *of asymmetric semi-alternating edges, then* \mathbf{M}\mathbf{i}{{n}^{\prime}}_{\mathcal{G}}=\mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}}\backslash {C}_{\mathcal{G}}\cup \left({C}_{{\mathcal{G}}^{\prime}}\backslash \left\{{e}_{i}^{\prime},\phantom{\rule{2.77695pt}{0ex}}{e}_{j}^{\prime}\right\}\right).

Proof: The type of \star \left({e}_{i}^{\prime},{e}_{j}^{\prime}\right) is S1 or S3 and the other star has type S4 or S5. By M2 {e}_{i}^{\prime} and {e}_{j}^{\prime} are not optimal.

**Lemma 4** (NE2). *If* \mathbf{M}i{\mathbf{n}}_{\mathcal{G}}\cap {C}_{\mathcal{G}}=\left\{{e}_{i}\right\} *for some i >* 0 *and the center is symmetric after the NNI operation then* \mathbf{M}i{{n}^{\prime}}_{\mathcal{G}}=\mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}}\backslash \left\{{e}_{i}\right\}\cup \star \left({e}_{0}^{\prime},\phantom{\rule{2.77695pt}{0ex}}{e}_{j}^{\prime}\right).

Proof: In this case {e}_{0}^{\prime} has two arrows and \star \left({e}_{0}^{\prime},{e}_{i}^{\prime}\right) is of type S5.

**Lemma 5**. *Assume that* \mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}}\cap {C}_{\mathcal{G}}=\left\{{e}_{0},{e}_{i},{e}_{j}\right\}*, where i ≠* 0,

**(NE3)** *If both e*_{
i
} *and e*_{
j
} *are symmetric then* \mathbf{M}i{\mathbf{n}}_{{\mathcal{G}}^{\prime}}=\mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}}\backslash {C}_{\mathcal{G}}\cup {C}_{{\mathcal{G}}^{\prime}},

**(NE4)** *If e*_{
j
} *is asymmetric and* {e}_{0}^{\prime} *is symmetric then* \mathbf{M}\mathbf{i}{\mathbf{n}}_{{\mathcal{G}}^{\prime}}=\mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}}\backslash {C}_{\mathcal{G}}\cup \star \left({e}_{0}^{\prime},\phantom{\rule{2.77695pt}{0ex}}{e}_{i}^{\prime}\right).

**(NE5)** *If both e*_{
j
} *and* {e}_{0}^{\prime} *are asymmetric then* \mathbf{M}\mathbf{i}{\mathbf{n}}_{{\mathcal{G}}^{\prime}}=\mathbf{M}\mathbf{i}{\mathbf{n}}_{\mathcal{G}}\backslash {C}_{\mathcal{G}}\cup \left\{{e}_{i}^{\prime}\right\}.

Proof: Note that {*e*_{0}, *e*_{
i
}*, e*_{
j
}} must be a star in \mathcal{G}\cdot \left(\mathsf{\text{NE}}3\right)\star \left({e}_{i},{e}_{j}\right) has type S4 or S5. After the transformation the two stars \star \left({e}_{0}^{\prime},{e}_{i}^{\prime}\right) and \star \left({e}_{0}^{\prime},{e}_{j}^{\prime}\right) have type S5. Both are optimal in {\mathcal{G}}^{\prime}\cdot \left(\mathsf{\text{NE4}}\right)\star \left({e}_{i},{e}_{j}\right) has type S5. After the transformation \star \left({e}_{0}^{\prime},{e}_{i}^{\prime}\right) has type S5 and \star \left({e}_{0}^{\prime},{e}_{j}^{\prime}\right) has type S3. Only the first is optimal in {\mathcal{G}}^{\prime}\cdot \left(\mathsf{\text{NE5}}\right)\star \left({e}_{i},{e}_{j}\right) has type S5 while the other star in {C}_{\mathcal{G}} has type S3. After the transformation only {e}_{i}^{\prime} remains symmetric in {C}_{{\mathcal{G}}^{\prime}} therefore it is the only optimal edge in {C}_{{\mathcal{G}}^{\prime}}.

*Computing the optimal cost (f)*. Observe that from Lemmas 2-5 at least one optimal edge remains optimal after the NNI operation. Therefore, to compute the difference in costs between optimal rootings of and {\mathcal{G}}^{\prime} we start with the cost analysis for the rootings of such edge.

First, we introduce a function for computing the cost differences. Consider three nodes *x, y, z* of some rooted gene tree such that *x* and *y* are siblings and the parent of them (denoted by *xy*), is a sibling of *z*. In other words we can denote this subtree by ((*x, y*), *z*). Then, the partial contribution of ((*x, y*), *z*) to the total weighted mutation cost can be described as follows: {\sum}_{a\in \mathcal{S}}\alpha *\left({\xi}^{D}\left(xy,\phantom{\rule{2.77695pt}{0ex}}a\right)+{\xi}^{D}\left(xyz,\phantom{\rule{2.77695pt}{0ex}}a\right)\right)+\beta *\left({\xi}^{L}\left(xy,\phantom{\rule{2.77695pt}{0ex}}a\right)+{\xi}^{L}\left(xyz,\phantom{\rule{2.77695pt}{0ex}}a\right)\right). Assume that *x, y* and *z* are mapped into *a, b* and *c* (from the species tree), respectively. It can be proved from the defnition of *ξ*^{D} and *ξ*^{L} that the above contribution equals: *ϕ*(*a, b, c*) = *α* * (*D*(*a, b*) + *D*(*a* + *b, c*)) + *β* * (*L*(*a, b*) + *L*(*a* + *b, c*)). Now, assume that a single NNI operation changes ((*x, y*), *z*)) into (*x*, (*y, z*)). It should be clear that the cost difference is given by: Δ_{3}(*a, b, c*) = *ϕ*(*c, b, a*) - *ϕ*(*a, b, c*). Similarly, we can define a cost difference when a single NNI operation changes ((*x, y*), (*z, v*)) into ((*x, v*), (*y, z*)). Assume, that *v* is mapped into *d*. Then, the cost contribution of the first subtree is *ϕ*'(*a, b, c, d*) = *ϕ*(*a, b, c* + *d*) + *α* * (*D*(*c, d*) + *β* * *L*(*c, d*). The cost difference is given by: Δ_{4}(*a, b, c, d*) = *ϕ*'(*a, d, b, c*) - *ϕ*'(*a, b, c, d*).

**Lemma 6**. *If the center edge is optimal and remains optimal after the NNI operation then the cost difference equals* Δ_{4}(*a*_{1}, *a*_{2}, *a*_{3}, *a*_{4})*, where a*_{
i
} *(for i* = 1, 2, 3, 4*) is the mapping as indicated in* Figure 3.

As mentioned the above lemma can be proved by comparing the rootings placed on the center edges in and {\mathcal{G}}^{\prime}. Lemma 6 gives a solution for cases: EQ2, EQ3, NE1 and NE3. The next lemma gives a solution for the remaining cases.

**Lemma 7**. *If for some i >* 0 *there exists an optimal edge in T*_{
i
} ∪ {*e*_{
i
}} *that remains optimal after the NNI operation (under assumption that e*_{
i
} *is replaced by* {e}_{i}^{\prime}*) then the cost difference is* Δ_{3}(*a*_{4}, *a*_{3}, *a*_{2}) *if i* = 1, Δ_{3}(*a*_{3}, *a*_{4,} *a*_{1}) *if i* = 2, Δ_{3}(*a*_{2}, *a*_{1}, *a*_{4}) *if i* = 3 *and* Δ_{3}(*a*_{1}, *a*_{2}, *a*_{3}) *if i* = 4.

Similarly to Lemma 6 we can prove Lemma 7 by comparing the rootings of *e*_{
i
} and {e}_{i}^{\prime}.

*Error correction algorithm*. Finally, we can present the algorithm for computing the optimal weighted mutation cost for a given gene tree and its *k*-NNI neighborhood. See Figure 4 for details. It should be clear that the complexity of this algorithm is O\left(|\mathcal{G}{|}^{k}+\text{max}\left(\left|\mathcal{G}\right|,\left|\mathcal{S}\right|\right)\right). We write that a gene tree has *errors* if the optimal cost is computed for one of its NNI variants. Otherwise, we write that a gene tree *does not require corrections*. Please note that it for a special case of *k* = 1, this algorithm is linear in time (see also our preliminary article [22]).

### General reconstruction problems

We present several approaches to problems of error correction and phylogeny reconstruction. Let us assume that {\sigma}_{\alpha ,\beta ,k}\left(\mathcal{S},\phantom{\rule{2.77695pt}{0ex}}\mathcal{G}\right) is the cost computed by algorithm from Figure 4, where *α, β* > 0, *k* ≥ 0, is a rooted species tree and is an unrooted gene tree.

**Problem 1** (*k* NNIC). *Given a rooted species tree* *and a set of unrooted gene trees, G compute the total cost* {\sum}_{\mathcal{G}\in G}{\sigma}_{\alpha ,\beta ,k}\left(\mathcal{S},\phantom{\rule{2.77695pt}{0ex}}\mathcal{G}\right).

The *k* NNIC problem can be solved in polynomial time by an iterative application of our algorithm. Additionally, we can reconstruct the optimal rootings as well as the correct topology of each gene tree. Please note that for *k* = 0 (no error correction), we have the cost inference problem for the reconciliation of an unrooted gene tree with a rooted species tree [18].

**Problem 2** (*k* NNIST). *Given a set of unrooted gene trees G find the species tree* *that minimizes the total cost*{\sum}_{\mathcal{G}\in G}{\sigma}_{\alpha ,\beta ,k}\left(S,\phantom{\rule{2.77695pt}{0ex}}\mathcal{G}\right).

The complexity of the *k* NNIST problem is unknown. However, similar problems for the duplication model are NP-hard [13]. Therefore we developed heuristics for the *k* NNIST problem to use them in our experiments.

In applications there is typically no need to search over all NNI variants of a gene tree. For instance, a good candidate for an NNI operation is *a weak edge*. A weak edge is usually defined on the basis of its length, where short length indicates weakness. To formalize this property, let us assume that each edge in a gene tree has length. We call an edge *e* in *weak* if the length of *e* is smaller than *ω*, where *ω* is a non-negative real. Now we can define variants of *k* NNIC and *k* NNIST denoted by *ω*-*k* NNIC and *ω*-*k* NNIST, respectively, where the NNI operations are performed on weak edges only. These straighforward definitions are omitted. Please note that the time complexity of the algorithm with NNIs limited to weak edges is O\left({l}^{k}+\text{max}\left(\left|\mathcal{G}\right|,\left|\mathcal{S}\right|\right)\right), where *l* is the number of weak edges in .

### Software

The unrooted reconciliation algorithm [18] and its data structures are implemented in program URec [23]. Our algorithm partially depends on theses data structures and therefore was implemented as a significantly extended version of URec. Additionally, we implemented a hill climbing heuristic to solve *k* NNIST and *ω*-kNNIST.

Software and datasets from our experiments are made freely available through http://bioputer.mimuw.edu.pl/~gorecki/ec.

### Experimental results and discussion

#### Data preparation

First, we inferred 4133 unrooted gene trees with branch lengths from nine yeast genomes contained in the Genolevures 3 data set [24], which contains protein sequences from the following nine yeast species: *C. glabrata* (4957 protein sequences, abbreviation CAGL), *S. cerevisiae* (5396, SACE), *Z. rouxii* (4840, ZYRO), *S. kluyveri* (5074, SAKL), *K. thermotolerans* (4933, KLTH), *K. lactis* (4851, KLLA), *Y. lipolytica* (4781, YALI), *D. hansenii* (5006, DEHA) and *E. gossypii* (4527, ERGO).

We aligned the protein sequences of each gene family by using the program TCoffee [25] using the default parameter setting. Then maximum likelihood (unrooted) gene trees were computed from the alignments by using proml from the phylip software package. The original species tree of these yeasts [24], here denoted by G3, is shown in Figure 5.

#### Inferring optimal species trees

The optimal species tree reconstructed with error corrections (1NNIST optimization problem) is depicted in Figure 5 and denoted by 1NNIEC. This tree differs from G3 in the rooting and in the middle clade with KLLA and ERGO. Additionally, we inferred by the heuristic an optimal species tree, denoted here by NOEC, with no error corrections (0NNIST optimization). All the trees from this figure are highly scored in each of the optimization schemas.

#### From weak edges to species trees

In the previous experiment, the NNI operations were performed on almost every gene tree in the optimal solution and with no restrictions on the edges. In order to reconstruct the trees more accurately, we performed experiments for *ω*-*k* NNIST optimization with various *ω* parameters and subsets of gene trees. The filtering of gene trees was determined by an integer *μ* > 0 that defines the maximum number of allowed weak edges in a single gene tree. Each gene tree that did not satisfy such condition was rejected.

Figures 6 and 7 depict a summary of error correction experiments for weak edges. For each *ω* and *μ* we performed 20 runs of the *ω*-*k* NNIST heuristic for finding the optimal species tree in the set of gene trees filtered by *μ*. The optimal species trees are depicted in the diagram, where each cell represents the result of a single *ω*-*k* NNIST experiment. We observed that G3, 1NNIEC and NOEC are significantly well represented in the set of optimal species trees in *ω*-1NNIST experiments, while in *ω*-2NNIST and *ω*-3NNIST experiments only G3 and NOEC were detected. Note that the original yeast phylogeny (G3, black squares in Figures 6 and 7) is inferred for *ω* = 0.1-0.2 (in other words approx. 30-40% of edges are weak, see Figure 8) and *μ* ≥ 10 in most experiments. In particular for *ω* = 0.15 and *μ* = 10, 364 gene trees were rejected (see Figure 9). These results significantly support the G3 phylogeny. Please note that the results for the standard unrooted reconciliation algorithms without error correction are located in the first column of diagrams (*ω* = 0).

#### From trusted species tree to weak edges in gene trees - automated and manual curation

Assume that the set of unrooted gene trees and the rooted (trusted) species tree are given. Then we can state the following problem: find *ω* and *μ* such that is the optimal species tree in *ω*-NNIST problem for the set of gene trees filtered by *μ*. For instance in our dataset, if we assume that G3 is a given correct phylogeny of yeasts, then from the diagrams (Figure 6 and 7) one can determine appropriate values of *ω* and *μ* that yield G3 as optimal. In other words we can automatically determine weak edges by *ω* and filter gene trees by *μ*. This approach can be applied in tree curation procedures to correct errors in an automated way as well as to find candidates (rejected trees) for further manual curation. For instance, in the previous case, when *ω* = 0.1 and *μ* = 10, we have 3164 trees that can be corrected and rooted by our algorithm, while the 364 rejected trees could be candidates for further manual correction.

## Discussion

We present novel theoretical and practical results on the problem of error correction and phylogeny reconstruction. In particular, we describe a polynomial time and space algorithm that simultaneously solves the problem of correction topological errors in unrooted gene trees and the problem of rooting unrooted gene trees. The algorithm allows us to perform efficiently experiments on truly large-scale datasets available for yeast genomes. Our experiments suggest that our algorithm can be used to (i) detect errors, (ii) to infer a correct phylogeny of species under the presence of weak edges in gene trees, and (iii) to help in tree curation procedures.

## Conclusion

We introduced a novel polynomial time algorithm for error-corrected and unrooted gene tree reconciliation. Experiments on yeast genomes suggests that an implementation of our algorithm can greatly improve on the accuracy of gene tree reconciliation, and thus, curate error-prone gene trees. Moreover, we use our error-corrected reconciliation to make the gene duplication problem, a standard application of gene tree reconciliation, more robust. We conjecture that the error-corrected gene duplication problem is intrinsically hard to solve, since the gene duplication problem is already NP-hard. Therefore, we introduced an effective heuristic for error-corrected gene duplication problem. Our experimental results for a wide range of error-correction tests on yeasts phylogeny show that our error-corrected reconciliations result in improved predictions of invoked gene duplication and loss events that then allow to infer more accurate phylogenies.

The presented error correction is based on gene-species tree reconciliation using gene duplication and loss. However, there are other major evolutionary mechanism that infer gene tree topologies that are inconsistent with the actual species tree topology, like horizontal gene transfer and deep coalescence. Gene tree reconciliation using these mechanisms is highly sensitive to topological error, similar to gene tree reconciliation under gene duplication and loss. Future work will focus on the development of algorithms that can also reconcile unrooted and erroneous gene trees using horizontal gene transfer and deep coalescence.

## References

Graur D, Li WH: Fundamentals of Molecular Evolution. 2000, Sinauer Associates, [http://www.amazon.com/exec/obidos/redirect?tag=citeulike07-20\&path=ASIN/0878932666]2

Page RDM: Maps between trees and cladistic analysis of historical associations among genes, organisms, and areas. Systematic Biology. 1994, 43: 58-77.

Bonizzoni P, Della Vedova G, Dondi R: Reconciling a gene tree to a species tree under the duplication cost model. Theoretical Computer Science. 2005, 347 (1-2): 36-53. 10.1016/j.tcs.2005.05.016.

Eulenstein O, Mirkin B, Vingron M: Duplication-Based Measures of Difference Between Gene and Species Trees. J Comput Biol. 1998, 5: 135-148. 10.1089/cmb.1998.5.135.

Goodman M, Czelusniak J, Moore GW, Romero-Herrera 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 (2): 132-163. 10.2307/2412519.

Mirkin B, Muchnik IB, Smith TF: A Biologically Consistent Model for Comparing Molecular Phylogenies. J Comput Biol. 1995, 2 (4): 493-507. 10.1089/cmb.1995.2.493.

Sanderson M, McMahon M: Inferring angiosperm phylogeny from EST data with widespread gene duplication. BMC Evolutionary Biology. 2007, 7 (Suppl 1): [http://dx.doi.org/10.1186/1471-2148-7-S1-S3]

Bansal MS, Eulenstein O: The multiple gene duplication problem revisited. Bioinformatics. 2008, 24 (13): i132-8. 10.1093/bioinformatics/btn150.

Fellows MR, Hallett MT, Stege U: On the Multiple Gene Duplication Problem. ISAAC, Volume 1533 of LNCS. Edited by: Chwa KY, Ibarra OH, Springer. 1998, 347-356.

Guigó R, Muchnik IB, Smith TF: Reconstruction of ancient molecular phylogeny. Molecular Phylogenetics and Evolution. 1996, 6 (2): 189-213. 10.1006/mpev.1996.0071.

Behzadi B, Vingron M: Reconstructing Domain Compositions of Ancestral Multi-domain Proteins. Comparative Genomics, Volume 4205 of LNCS. Edited by: Bourque G, El-Mabrouk N. 2006, Springer, 1-10.

Bansal MS, Burleigh GJ, Eulenstein O, Wehe A: Heuristics for the Gene-Duplication Problem: A Θ(

*n*) Speed-Up for the Local Search. RECOMB, Volume 4453 of LNCS. 2007, Springer, 238-252.Ma B, Li M, Zhang L: From Gene Trees to Species Trees. SIAM Journal on Computing. 2000, 30 (3): 729-752. 10.1137/S0097539798343362.

Page RDM: GeneTree: comparing gene and species phylogenies using reconciled trees. Bioinformatics. 1998, 14 (9): 819-820. 10.1093/bioinformatics/14.9.819.

Hahn MW: Bias in phylogenetic tree reconciliation methods: implications for vertebrate genome evolution. Genome biology. 2007, 8 (7): R141-10.1186/gb-2007-8-7-r141. [http://dx.doi.org/10.1186/gb-2007-8-7-r141]

Chen K, Durand D, Farach-Colton M: NOTUNG: a program for dating gene duplications and optimizing gene family trees. J Comput Biol. 2000, 7 (3-4): 429-447. 10.1089/106652700750050871.

Durand D, Halldorsson BV, Vernot B: A Hybrid Micro-Macroevolutionary Approach to Gene Tree Reconstruction. J Comput Biol. 2006, 13 (2): 320-335. 10.1089/cmb.2006.13.320. [http://dx.doi.org/10.1089/cmb.2006.13.320]

Górecki P, Tiuryn J: Inferring phylogeny from whole genomes. Bioinformatics. 2007, 23 (2): e116-22. 10.1093/bioinformatics/btl296.

Wehe A, Bansal MS, Burleigh GJ, Eulenstein O: Dup-Tree: a program for large-scale phylogenetic analyses using gene tree parsimony. Bioinformatics. 2008, 24 (13): 1540-1541. 10.1093/bioinformatics/btn230.

Eulenstein O, Huzurbazar S, Liberles D: Reconciling phylogenetic trees. Evolution After Gene Duplication. 2010, Dittmar, Liberles, Wiley

Bender MA, Farach-Colton M: The LCA Problem Revisited LATIN, Volume 1776 of LNCS. Edited by: Gonnet GH, Panario D, Viola A. 2000, Springer, 88-94.

Górecki P, Eulenstein O: A Linear Time Algorithm for Error-Corrected Reconciliation of Unrooted Gene Trees. Bioinformatics Research and Applications, Volume 6674 of Lecture Notes in Computer Science. Edited by: Chen J, Wang J, Zelikovsky A. 2011, Springer Berlin/Heidelberg, 148-159.

Górecki P, Tiuryn J: URec: a system for unrooted reconciliation. Bioinformatics. 2007, 23 (4): 511-512. 10.1093/bioinformatics/btl634.

Sherman DJ, Martin T, Nikolski M, Cayla C, Souciet JL, Durrens P: Gènolevures: protein families and synteny among complete hemiascomycetous yeast proteomes and genomes. Nucleic Acids Research. 2009, 37 (suppl 1): D550-D554. [http://nar.oxfordjournals.org/content/37/suppl_1/D550.abstract]

Notredame C, Higgins DG, Jaap H: T-coffee: a novel method for fast and accurate multiple sequence alignment. J Mol Biol. 2000, 302: 205-217. 10.1006/jmbi.2000.4042. [http://dx.doi.org/10.1006/jmbi.2000.4042]

## Acknowledgements

The reviewers have provided several valuable comments that have improved the presentation. This work was conducted in parts with support from the Gene Tree Reconciliation Working Group at NIMBioS through NSF award #EF-0832858. PG was partially supported by the grant of MNiSW (N N301 065236) and OE was 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.

## Author information

### Authors and Affiliations

### Corresponding author

## Additional information

### Competing interests

The authors declare that they have no competing interests.

### Authors' contributions

PG and OE were responsible for algorithm design and writing the paper. PG implemented the programs, and performed the experimental evaluation and the analysis of the results. Both authors read and approved the final manuscript.

## Rights and permissions

This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## About this article

### Cite this article

Górecki, P., Eulenstein, O. Algorithms: simultaneous error-correction and rooting for gene tree reconciliation and the gene duplication problem.
*BMC Bioinformatics* **13**
(Suppl 10), S14 (2012). https://doi.org/10.1186/1471-2105-13-S10-S14

Published:

DOI: https://doi.org/10.1186/1471-2105-13-S10-S14