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

  • Pawel Górecki1Email author and

    Affiliated with

    • Oliver Eulenstein2

      Affiliated with

      BMC Bioinformatics201213(Suppl 10):S14

      DOI: 10.1186/1471-2105-13-S10-S14

      Published: 25 June 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.
      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_Fig1_HTML.jpg
      Figure 1

      Rooted reconciliation. An lca-mapping M from the gene tree http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif into the species tree http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif and the corresponding embedding. M is shown for the internal nodes of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif .

      Related work

      Gene tree reconciliation is a well-studied method for resolving topological inconsistencies between a gene tree and a trusted species tree [27]. 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 [810], reconstructing domain decompositions [11], and species supertree construction [8, 1214].

      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, 1619].

      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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq3_HTML.gif be a species tree. http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif 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, http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq4_HTML.gif , 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 ab or ba 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, http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq4_HTML.gif siblings when a + b is a parent of a and b. For a, http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq4_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq5_HTML.gif be the least common ancestor (lca) mapping, from rooted http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif into http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif that preserves the labeling of the leaves. Formally, if v is a leaf in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif then M(v) is the node in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif labeled by the label of v. If v is internal node in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq6_HTML.gif which for all nodes http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq7_HTML.gif , http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq8_HTML.gif assigns a real ξ(v, a) representing a contribution to node a which comes from v when reconciling http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif with http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif . Having ξ we can define http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq9_HTML.gif to be a total contribution from v in the reconciliation of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif with http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif . We call κ a contribution function. Finally, http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq10_HTML.gif is the total cost of reconciliation of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif with http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif .

      Now we present examples of cost functions that are used in the duplication model. We assume that if v is an internal node in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif then w 1 and w 2 are its children. The Duplication cost function is defined as follows: ξ D (v, a) = 1 if http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq11_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq11_HTML.gif and aSb (M(w 1), M(w 2)), and ξ L (v, a) = 0 otherwise. It can be proved that if http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq11_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq11_HTML.gif and M(v) = M(w i ) = a for some i, and ξ D (v, a) = 0 otherwise. Loss cost function: ξ L (v, a) = 1 if http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq11_HTML.gif and aSb(M(w 1), M(w 2)), and ξ L (v, a) = 0 otherwise. It can be proved that if http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq11_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq12_HTML.gif 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) http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif with http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif . 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq13_HTML.gif is an unrooted gene tree. We define a rooting of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif by selecting an edge http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq14_HTML.gif on which the root is to be placed. Such a rooted tree will be denoted by http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq15_HTML.gif , where v * is a new node defining the root. To distinguish between rootings of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif , 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq15_HTML.gif is independent of e. Without loss of generality the following is assumed: (A1) http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif and http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif have at least one internal node and (A2) M e (v *)=⊤; that is, the root of every rooting is mapped into the root of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif (we may always consider the subtree of the species tree rooted in M e (v *) with no change of the cost).

      First, we transform http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif into a directed graph http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq16_HTML.gif where http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq17_HTML.gif . In other words each edge 〈v, w〉 in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif is replaced in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif by a pair of directed edges 〈v, w〉 and 〈w, v〉.

      Edges in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif are labeled by nodes of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif as follows. If http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq19_HTML.gif is a leaf labeled by a, then the edge http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq20_HTML.gif is labeled by a. When v is an internal node in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif we assume that 〈w 1, v〉 and 〈w 2, v〉 are labeled by b 1 and b 2, respectively. Then the edge http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq21_HTML.gif such that w 3w 1 and w 3w 2 is labeled by b 1 + b 2. Such labeling will be used to explore mappings of rootings of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif . An edge {v, w} in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif is called asymmetric if exactly one of the labels of 〈v, w〉 and 〈w, v〉 in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif is equal to ⊤, otherwise it is called symmetric.

      Every internal node v, and its neighbors in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif define a subtree of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq22_HTML.gif , 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.
      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_Fig2_HTML.jpg
      Figure 2

      Unrooted reconciliation. a) A star in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif . 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif be a species tree and http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif be unrooted gene tree. The set of optimal edges, that is, candidates for best rootings, is defined as follows: http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq23_HTML.gif , where http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq24_HTML.gif is the total cost for the weighted mutation cost defined by http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq25_HTML.gif , e is an edge in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif and α, β are two positive reals. Then (M1) if http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq26_HTML.gif , then http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq27_HTML.gif consists of all edges present in all stars of type S4 or S5, (M2) if http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq28_HTML.gif , then http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq27_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq27_HTML.gif . More precisely, the star edges outside http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq27_HTML.gif are asymmetric and share the same direction. Thus, to find an optimal edge it is sufficient to follow the direction of non ⊤ edges in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif .

      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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif ) can be computed in constant time after an initial preprocessing step requiring http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq29_HTML.gif time. Other structures like http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif with the labeling can be computed in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq30_HTML.gif time. The same complexity has the procedure of finding an optimal edge in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif . 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.
      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_Fig3_HTML.jpg
      Figure 3

      NNI. A single NNI on http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif and http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif . On the left e i and http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq31_HTML.gif (for i = 0, ... , 4) denote edges in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif and its NNI-neighbor http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq32_HTML.gif , respectively. On the right each node a i denote the labeling of edges in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif . Notation http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq33_HTML.gif denote the lca-mapping of complementary subtrees, for instance, http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq34_HTML.gif , etc. For brevity, we omit each subtree T i attached to w i 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif , (b) lca structure for http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif that allows to answer lca-queries in constant time, (c) an unrooted gene tree http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif , (d) http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif with the labeling of edges, (e) http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq27_HTML.gif - 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq35_HTML.gif . 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq36_HTML.gif into http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq37_HTML.gif , where T i -s are (rooted) subtrees of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif . The edge that connects the roots of (T 1, T 2) and (T 3, T 4) in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq38_HTML.gif and http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq39_HTML.gif in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq40_HTML.gif .

      An NNI operation is depicted in Figure 3 with the transformation of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif into http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq41_HTML.gif . The notation will be used from now on. Note that there is a second NNI operation, when http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif . Observe that the labels of the incoming and outgoing edges attached to each w i in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq18_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq42_HTML.gif . 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq43_HTML.gif .

      For convenience, assume that the NNI operation replaces e i with http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq44_HTML.gif as indicated in Figure 3. We call two disjoint edges from http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq45_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq89_HTML.gif 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:

      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq47_HTML.gif

      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq48_HTML.gif and each pair of semi-alternating edges contains at least one symmetric edge,

      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq90_HTML.gif consists of only the center edge,

      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq50_HTML.gif for some i >0 and the center is asymmetric after the NNI operation.

      Proof: (EQ1) All edges in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq51_HTML.gif are asymmetric (2 stars S1). Then, after the NNI operation http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq52_HTML.gif is asymmetric and ( http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq53_HTML.gif has 2 stars S1). (EQ2) http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq54_HTML.gif consists of 2 stars of type S4/S5 and at most two asymmetric edges. It follows from EQ2 that the asymmetric edges in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq53_HTML.gif cannot form a star of type other than S5. Together with M1 it follows that http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq53_HTML.gif is optimal. (EQ3) By M1 the center is symmetric in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif . It remains symmetric after NNI. From C1 and M2, http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq55_HTML.gif consists of the center edge. (EQ4) Note, that the type of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq56_HTML.gif is S1, S2 or S3.

      Lemma 3 (NE1). If http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq57_HTML.gif and there exists a pair {e i , e j } of asymmetric semi-alternating edges, then http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq58_HTML.gif .

      Proof: The type of http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq59_HTML.gif is S1 or S3 and the other star has type S4 or S5. By M2 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq60_HTML.gif and http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq61_HTML.gif are not optimal.

      Lemma 4 (NE2). If http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq62_HTML.gif for some i >0 and the center is symmetric after the NNI operation then http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq63_HTML.gif .

      Proof: In this case http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq64_HTML.gif has two arrows and http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq65_HTML.gif is of type S5.

      Lemma 5. Assume that http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq66_HTML.gif , where i ≠0,

      (NE3) If both e i and e j are symmetric then http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq67_HTML.gif ,

      (NE4) If e j is asymmetric and is symmetric then http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq68_HTML.gif .

      (NE5) If both e j and are asymmetric then http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq91_HTML.gif .

      Proof: Note that {e 0, e i , e j } must be a star in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq70_HTML.gif has type S4 or S5. After the transformation the two stars http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq71_HTML.gif and http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq72_HTML.gif have type S5. Both are optimal in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq73_HTML.gif has type S5. After the transformation http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq74_HTML.gif has type S5 and http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq75_HTML.gif has type S3. Only the first is optimal in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq76_HTML.gif has type S5 while the other star in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq51_HTML.gif has type S3. After the transformation only http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq77_HTML.gif remains symmetric in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq78_HTML.gif therefore it is the only optimal edge in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq79_HTML.gif .

      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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif and http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq80_HTML.gif 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: http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq81_HTML.gif . 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif and http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq80_HTML.gif . 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 ) 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq83_HTML.gif .

      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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq84_HTML.gif . 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]).
      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_Fig4_HTML.jpg
      Figure 4

      Algorithm. Optimal weighted cost for http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif and its k-NNI neighborhood.

      General reconstruction problems

      We present several approaches to problems of error correction and phylogeny reconstruction. Let us assume that http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq85_HTML.gif is the cost computed by algorithm from Figure 4, where α, β > 0, k ≥ 0, http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif is a rooted species tree and http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif is an unrooted gene tree.

      Problem 1 (kNNIC). Given a rooted species tree http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif and a set of unrooted gene trees, G compute the total cost http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq86_HTML.gif .

      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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif that minimizes the total cost http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq87_HTML.gif .

      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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif has length. We call an edge e in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq88_HTML.gif , where l is the number of weak edges in http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq1_HTML.gif .

      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.

      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.
      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_Fig5_HTML.jpg
      Figure 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).
      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_Fig6_HTML.jpg
      Figure 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.

      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_Fig7_HTML.jpg
      Figure 7

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

      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_Fig8_HTML.jpg
      Figure 8

      Branch lengths. Histogram of branch lengths.

      http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_Fig9_HTML.jpg
      Figure 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 http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif are given. Then we can state the following problem: find ω and μ such that http://static-content.springer.com/image/art%3A10.1186%2F1471-2105-13-S10-S14/MediaObjects/12859_2012_5219_IEq2_HTML.gif 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.

      Declarations

      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/​bmcbioinformatic​s/​supplements/​13/​S10.

      Authors’ Affiliations

      (1)
      Institute of Informatics, University of Warsaw
      (2)
      Department of Computer Science, Iowa State University

      References

      1. Graur D, Li WH: [http://​www.​amazon.​com/​exec/​obidos/​redirect?​tag=​citeulike07–20\&​path=​ASIN/​0878932666] Fundamentals of Molecular Evolution. 2nd edition. Sinauer Associates; 2000.
      2. Page RDM: Maps between trees and cladistic analysis of historical associations among genes, organisms, and areas. Systematic Biology 1994, 43:58–77.
      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.View Article
      4. Eulenstein O, Mirkin B, Vingron M: Duplication-Based Measures of Difference Between Gene and Species Trees. J Comput Biol 1998, 5:135–148.PubMedView Article
      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.View Article
      6. Mirkin B, Muchnik IB, Smith TF: A Biologically Consistent Model for Comparing Molecular Phylogenies. J Comput Biol 1995,2(4):493–507.PubMedView Article
      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] BMC Evolutionary Biology 2007.,7(Suppl 1):
      8. Bansal MS, Eulenstein O: The multiple gene duplication problem revisited. Bioinformatics 2008,24(13):i132–8.PubMedView Article
      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.
      10. Guigó R, Muchnik IB, Smith TF: Reconstruction of ancient molecular phylogeny. Molecular Phylogenetics and Evolution 1996,6(2):189–213.PubMedView Article
      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.View Article
      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.
      13. Ma B, Li M, Zhang L: From Gene Trees to Species Trees. SIAM Journal on Computing 2000,30(3):729–752.View Article
      14. Page RDM: GeneTree: comparing gene and species phylogenies using reconciled trees. Bioinformatics 1998,14(9):819–820.PubMedView Article
      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] Genome biology 2007,8(7):R141.PubMedView Article
      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.PubMedView Article
      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] J Comput Biol 2006,13(2):320–335.PubMedView Article
      18. Górecki P, Tiuryn J: Inferring phylogeny from whole genomes. Bioinformatics 2007,23(2):e116–22.PubMedView Article
      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.PubMedView Article
      20. Eulenstein O, Huzurbazar S, Liberles D: Reconciling phylogenetic trees. In Evolution After Gene Duplication. Dittmar, Liberles, Wiley; 2010.
      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.
      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.
      23. Górecki P, Tiuryn J: URec: a system for unrooted reconciliation. Bioinformatics 2007,23(4):511–512.PubMedView Article
      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] Nucleic Acids Research 2009,37(suppl 1):D550-D554.PubMedView Article
      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] J Mol Biol 2000, 302:205–217.PubMedView Article

      Copyright

      © Górecki and Eulenstein; licensee BioMed Central Ltd. 2012

      Advertisement