Email updates

Keep up to date with the latest news and content from BMC Bioinformatics and BioMed Central.

This article is part of the supplement: Selected articles from the 7th International Symposium on Bioinformatics Research and Applications (ISBRA'11)

Open Access Proceedings

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

Pawel Górecki1* and Oliver Eulenstein2

Author Affiliations

1 Institute of Informatics, University of Warsaw, Warsaw, 02-097, Poland

2 Department of Computer Science, Iowa State University, Ames, 50011, USA

For all author emails, please log on.

BMC Bioinformatics 2012, 13(Suppl 10):S14  doi:10.1186/1471-2105-13-S10-S14

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2105/13/S10/S14


Published:25 June 2012

© 2012 Górecki and Eulenstein; licensee 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.

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.

thumbnailFigure 1. Rooted reconciliation. An lca-mapping M from the gene tree into the species tree and the corresponding embedding. M is shown for the internal nodes of .

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(lk + 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M2">View MathML</a> 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, <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M3">View MathML</a>, 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, <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M3">View MathML</a>siblings when a + b is a parent of a and b. For a, <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M3">View MathML</a> 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M4">View MathML</a> 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M6">View MathML</a> which for all nodes <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M7">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M8">View MathML</a> assigns a real ξ(v, a) representing a contribution to node a which comes from v when reconciling with . Having ξ we can define <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M9">View MathML</a> to be a total contribution from v in the reconciliation of with . We call κ a contribution function. Finally, <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M10">View MathML</a> 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 w1 and w2 are its children. The Duplication cost function is defined as follows: ξD(v, a) = 1 if <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M11">View MathML</a> and M(v) = M(wi) = a for some i, and ξD(v, a) = 0 otherwise. The Loss cost function: ξL(v, a) = 1 if <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M11">View MathML</a> and a Sb (M(w1), M(w2)), and ξL(v, a) = 0 otherwise. It can be proved that if <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M11">View MathML</a> then κD(v) = D(M(w1), M(w2)) and κL (v) = L(M(w1), M(w2)) (in both cases 0 if v is a leaf).

The Duplication cost function is defined as follows: ξD(v, a) = 1 if <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M11">View MathML</a> and M(v) = M(wi) = a for some i, and ξD(v, a) = 0 otherwise. Loss cost function: ξL(v, a) = 1 if <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M11">View MathML</a> and a Sb(M(w1), M(w2)), and ξL(v, a) = 0 otherwise. It can be proved that if <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M11">View MathML</a> then κD(v) = D(M(w1), M(w2)) and κL(v) = L(M(w1), M(w2)) (in both cases 0 if v is a leaf).

Observe that a node <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M12">View MathML</a> 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M87">View MathML</a> is an unrooted gene tree. We define a rooting of by selecting an edge <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M13">View MathML</a> on which the root is to be placed. Such a rooted tree will be denoted by <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M14">View MathML</a>, 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M14">View MathML</a> is independent of e. Without loss of generality the following is assumed: (A1) and have at least one internal node and (A2) Me(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 Me(v*) with no change of the cost).

First, we transform into a directed graph <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M15">View MathML</a> where <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M16">View MathML</a>. In other words each edge 〈v, w〉 in is replaced in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a> by a pair of directed edges 〈v, w〉 and 〈w, v〉.

Edges in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a> are labeled by nodes of as follows. If <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M18">View MathML</a> is a leaf labeled by a, then the edge <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M19">View MathML</a> is labeled by a. When v is an internal node in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a> we assume that 〈w1, v〉 and 〈w2, v〉 are labeled by b1 and b2, respectively. Then the edge <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M20">View MathML</a> such that w3 w1 and w3 w2 is labeled by b1 + b2. 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a> is equal to ⊤, otherwise it is called symmetric.

Every internal node v, and its neighbors in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a> define a subtree of <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M21">View MathML</a>, called a star with a center v, as depicted in Figure 2. The edges 〈v, wi〉 are called outgoing, while the edges 〈wi, v〉 are called incoming. We will refer to the undirected edge {v, wi} as ei, for i = 1, 2, 3.

thumbnailFigure 2. Unrooted reconciliation. a) A star in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a>. b) Types of edges. c) All possible types of stars. We use simplified notation instead of the full topology.

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: <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M22">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M23">View MathML</a>is the total cost for the weighted mutation cost defined by <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M24">View MathML</a>, e is an edge in and α, β are two positive reals. Then (M1) if <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M25">View MathML</a>, then <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M88">View MathML</a> consists of all edges present in all stars of type S4 or S5, (M2) if <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M26">View MathML</a>, then <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M88">View MathML</a> 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M88">View MathML</a>. More precisely, the star edges outside <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M88">View MathML</a> are asymmetric and share the same direction. Thus, to find an optimal edge it is sufficient to follow the direction of non ⊤ edges in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a>.

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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M27">View MathML</a> time. Other structures like <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a> with the labeling can be computed in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M28">View MathML</a> 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.

thumbnailFigure 3. NNI. A single NNI on and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a>. On the left ei and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M83','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M83">View MathML</a> (for i = 0, ... , 4) denote edges in and its NNI-neighbor <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M84','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M84">View MathML</a>, respectively. On the right each node ai denote the labeling of edges in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a>. Notation <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M85','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M85">View MathML</a> denote the lca-mapping of complementary subtrees, for instance, <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M86">View MathML</a>, etc. For brevity, we omit each subtree Ti attached to wi in the left diagram.

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) <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a> with the labeling of edges, (e) <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M88">View MathML</a> - 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M29">View MathML</a>. 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M30">View MathML</a> into <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M31">View MathML</a>, where Ti-s are (rooted) subtrees of . The edge that connects the roots of (T1, T2) and (T3, T4) in is denoted by e0 and called the center edge. For each i = 1, 2, 3, 4 we assume the following: wi is the root of Ti, ei is the edge connecting wi with e0 and ai is the lca-mapping of Ti. Similarly, we define the center edge <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M32">View MathML</a>and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M33">View MathML</a> in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M34">View MathML</a>.

An NNI operation is depicted in Figure 3 with the transformation of <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a> into <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M35">View MathML</a>. The notation will be used from now on. Note that there is a second NNI operation, when is replaced with ((T1, T3), (T2, T4)). 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a>. Observe that the labels of the incoming and outgoing edges attached to each wi in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a> 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M17">View MathML</a> 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M36">View MathML</a>. 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M37">View MathML</a>.

For convenience, assume that the NNI operation replaces ei with <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M38">View MathML</a> as indicated in Figure 3. We call two disjoint edges from <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M39">View MathML</a>semi-alternating if they share a common node after the NNI operation. In Figure 3 {e1, e4} and {e2, e3} 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 ei is replaced by <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M40">View MathML</a> 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:

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M41">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M42">View MathML</a>and each pair of semi-alternating edges contains at least one symmetric edge,

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M43">View MathML</a>consists of only the center edge,

<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M44">View MathML</a>for some i >0 and the center is asymmetric after the NNI operation.

Proof: (EQ1) All edges in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M70">View MathML</a> are asymmetric (2 stars S1). Then, after the NNI operation <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M45">View MathML</a> is asymmetric and (<a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M46">View MathML</a> has 2 stars S1). (EQ2) <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M47">View MathML</a> consists of 2 stars of type S4/S5 and at most two asymmetric edges. It follows from EQ2 that the asymmetric edges in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M46">View MathML</a> cannot form a star of type other than S5. Together with M1 it follows that <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M46">View MathML</a> is optimal. (EQ3) By M1 the center is symmetric in . It remains symmetric after NNI. From C1 and M2, <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M48">View MathML</a> consists of the center edge. (EQ4) Note, that the type of <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M49">View MathML</a> is S1, S2 or S3.

Lemma 3 (NE1). If <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M50">View MathML</a>and there exists a pair {ei , ej} of asymmetric semi-alternating edges, then <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M51">View MathML</a>.

Proof: The type of <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M52">View MathML</a> is S1 or S3 and the other star has type S4 or S5. By M2 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M53">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M54">View MathML</a> are not optimal.

Lemma 4 (NE2). If <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M55">View MathML</a> for some i >0 and the center is symmetric after the NNI operation then <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M56">View MathML</a>.

Proof: In this case <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M57">View MathML</a> has two arrows and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M58">View MathML</a> is of type S5.

Lemma 5. Assume that <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M59">View MathML</a> , where i ≠ 0,

(NE3) If both ei and ej are symmetric then <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M60">View MathML</a>,

(NE4) If ej is asymmetric and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M57">View MathML</a> is symmetric then <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M61">View MathML</a>.

(NE5) If both ej and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M57">View MathML</a> are asymmetric then <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M62">View MathML</a>.

Proof: Note that {e0, ei , ej} must be a star in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M63">View MathML</a> has type S4 or S5. After the transformation the two stars <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M64">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M65">View MathML</a> have type S5. Both are optimal in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M66">View MathML</a> has type S5. After the transformation <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M67">View MathML</a> has type S5 and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M68">View MathML</a> has type S3. Only the first is optimal in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M69">View MathML</a> has type S5 while the other star in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M70">View MathML</a> has type S3. After the transformation only <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M71">View MathML</a> remains symmetric in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M72">View MathML</a> therefore it is the only optimal edge in <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M73">View MathML</a>.

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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M74">View MathML</a> 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: <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M75">View MathML</a>. 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(a1, a2, a3, a4), where ai (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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M74">View MathML</a>. 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 Ti ∪ {ei} that remains optimal after the NNI operation (under assumption that ei is replaced by <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M76">View MathML</a>) then the cost difference is Δ3(a4, a3, a2) if i = 1, Δ3(a3, a4, a1) if i = 2, Δ3(a2, a1, a4) if i = 3 and Δ3(a1, a2, a3) if i = 4.

Similarly to Lemma 6 we can prove Lemma 7 by comparing the rootings of ei and <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M77">View MathML</a>.

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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M78">View MathML</a>. 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]).

thumbnailFigure 4. Algorithm. Optimal weighted cost for and its k-NNI neighborhood.

General reconstruction problems

We present several approaches to problems of error correction and phylogeny reconstruction. Let us assume that <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M79">View MathML</a> 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 (kNNIC). Given a rooted species tree and a set of unrooted gene trees, G compute the total cost <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M80">View MathML</a>.

The kNNIC 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 (kNNIST). Given a set of unrooted gene trees G find the species tree that minimizes the total cost <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M81','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M81">View MathML</a>.

The complexity of the kNNIST problem is unknown. However, similar problems for the duplication model are NP-hard [13]. Therefore we developed heuristics for the kNNIST 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 kNNIC and kNNIST denoted by ω-kNNIC and ω-kNNIST, 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 <a onClick="popup('http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M82','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/13/S10/S14/mathml/M82">View MathML</a>, 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 kNNIST and ω-kNNIST.

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

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.

thumbnailFigure 5. Yeasts phylogeny. Species tree topologies. G3 - original phylogeny of Genolevures 3 data set [24]. 1NNIEC - optimal rooted species tree inferred from gene trees with all possible 1-NNI error corrections. NOEC - optimal species tree for the yeast gene trees with no NNI operations (cost 64413, no corrections). Rank denotes a position of a tree on the sorted list of the best trees. The trees below are inferred from other ω-kNNIST (see next figures). Please note that NOEC, G3, a1 and a2 are rooted variants of the same unrooted tree. Similar property holds for 1NNIEC, b1 and b2.

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 ω-kNNIST 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 ω-kNNIST 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 ω-kNNIST 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).

thumbnailFigure 6. ω-1NNIST and ω-2NNIST experiments. A summary of ω-1NNIST (top) and ω-2NNIST experiments (bottom) for ω = 0, 0.02, 0.04, ... , 0.98, μ = 2, 3, ... , 20. Optimal species trees found by the heuristics. Please note that in some cases two optimal trees were found.

thumbnailFigure 7. ω-3NNIST experiments. A summary of ω-3NNIST experiments for ω = 0, 0.02, 004, ... , 0.48 and μ = 2, 3, ... , 20.

thumbnailFigure 8. Branch lengths. Histogram of branch lengths.

thumbnailFigure 9. Rejected gene trees. The number of rejected trees as a function of μ and ω.

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.

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.

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.

References

  1. Graur D, Li WH: [http:/ / www.amazon.com/ exec/ obidos/ redirect?tag=citeulike07-20\&path=A SIN/ 0878932666] webcite

    Fundamentals of Molecular Evolution. 2nd edition. Sinauer Associates; 2000. OpenURL

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

    Systematic Biology 1994, 43:58-77. OpenURL

  3. 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. Publisher Full Text OpenURL

  4. Eulenstein O, Mirkin B, Vingron M: Duplication-Based Measures of Difference Between Gene and Species Trees.

    J Comput Biol 1998, 5:135-148. PubMed Abstract | Publisher Full Text OpenURL

  5. 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. Publisher Full Text OpenURL

  6. Mirkin B, Muchnik IB, Smith TF: A Biologically Consistent Model for Comparing Molecular Phylogenies.

    J Comput Biol 1995, 2(4):493-507. PubMed Abstract | Publisher Full Text OpenURL

  7. Sanderson M, McMahon M: Inferring angiosperm phylogeny from EST data with widespread gene duplication. [http://dx.doi.org/10.1186/1471-2148-7-S1-S3] webcite

    BMC Evolutionary Biology 2007., 7(Suppl 1) OpenURL

  8. Bansal MS, Eulenstein O: The multiple gene duplication problem revisited.

    Bioinformatics 2008, 24(13):i132-8. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  9. Fellows MR, Hallett MT, Stege U: On the Multiple Gene Duplication Problem.

    In ISAAC, Volume 1533 of LNCS Edited by Chwa KY, Ibarra OH, Springer. 1998, 347-356. OpenURL

  10. Guigó R, Muchnik IB, Smith TF: Reconstruction of ancient molecular phylogeny.

    Molecular Phylogenetics and Evolution 1996, 6(2):189-213. PubMed Abstract | Publisher Full Text OpenURL

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

  12. Bansal MS, Burleigh GJ, Eulenstein O, Wehe A: Heuristics for the Gene-Duplication Problem: A Θ(n) Speed-Up for the Local Search. In RECOMB, Volume 4453 of LNCS. Springer; 2007:238-252. OpenURL

  13. Ma B, Li M, Zhang L: From Gene Trees to Species Trees.

    SIAM Journal on Computing 2000, 30(3):729-752. Publisher Full Text OpenURL

  14. Page RDM: GeneTree: comparing gene and species phylogenies using reconciled trees.

    Bioinformatics 1998, 14(9):819-820. PubMed Abstract | Publisher Full Text OpenURL

  15. Hahn MW: Bias in phylogenetic tree reconciliation methods: implications for vertebrate genome evolution. [http://dx.doi.org/10.1186/gb-2007-8-7-r141] webcite

    Genome biology 2007, 8(7):R141. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  16. 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. PubMed Abstract | Publisher Full Text OpenURL

  17. Durand D, Halldorsson BV, Vernot B: A Hybrid Micro-Macroevolutionary Approach to Gene Tree Reconstruction. [http://dx.doi.org/10.1089/cmb.2006.13.320] webcite

    J Comput Biol 2006, 13(2):320-335. PubMed Abstract | Publisher Full Text OpenURL

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

    Bioinformatics 2007, 23(2):e116-22. PubMed Abstract | Publisher Full Text OpenURL

  19. 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. PubMed Abstract | Publisher Full Text OpenURL

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

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

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

  23. Górecki P, Tiuryn J: URec: a system for unrooted reconciliation.

    Bioinformatics 2007, 23(4):511-512. PubMed Abstract | Publisher Full Text OpenURL

  24. 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. [http://nar.oxfordjournals.org/content/37/suppl_1/D550.abstract] webcite

    Nucleic Acids Research 2009, 37(suppl 1):D550-D554. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  25. Notredame C, Higgins DG, Jaap H: T-coffee: a novel method for fast and accurate multiple sequence alignment. [http://dx.doi.org/10.1006/jmbi.2000.4042] webcite

    J Mol Biol 2000, 302:205-217. PubMed Abstract | Publisher Full Text OpenURL