An ILP solution for the gene duplication problem
- Wen-Chieh Chang^{1},
- Gordon J Burleigh^{2},
- David F Fernández-Baca^{1} and
- Oliver Eulenstein^{1}Email author
https://doi.org/10.1186/1471-2105-12-S1-S14
© Chang et al; licensee BioMed Central Ltd. 2011
Published: 15 February 2011
Abstract
Background
The gene duplication (GD) problem seeks a species tree that implies the fewest gene duplication events across a given collection of gene trees. Solving this problem makes it possible to use large gene families with complex histories of duplication and loss to infer phylogenetic trees. However, the GD problem is NP-hard, and therefore, most analyses use heuristics that lack any performance guarantee.
Results
We describe the first integer linear programming (ILP) formulation to solve instances of the gene duplication problem exactly. With simulations, we demonstrate that the ILP solution can solve problem instances with up to 14 taxa. Furthermore, we apply the new ILP solution to solve the gene duplication problem for the seed plant phylogeny using a 12-taxon, 6, 084-gene data set. The unique, optimal solution, which places Gnetales sister to the conifers, represents a new, large-scale genomic perspective on one of the most puzzling questions in plant systematics.
Conclusions
Although the GD problem is NP-hard, our novel ILP solution for it can solve instances with data sets consisting of as many as 14 taxa and 1, 000 genes in a few hours. These are the largest instances that have been solved to optimally to date. Thus, this work can provide large-scale genomic perspectives on phylogenetic questions that previously could only be addressed by heuristic estimates.
Keywords
Background
With recent advances in DNA sequencing technology, there is much interest in using genomic data sets to infer phylogenetic trees. However, evolutionary events such as gene duplication and loss, incomplete lineage sorting (deep coalescence), and lateral gene transfer can produce discordance between gene trees and the phylogeny of the species in which the genes evolve (e.g., [1]). The gene tree parsimony (GTP) problem [1–4] provides a direct approach to infer a species phylogeny from discordant gene trees. Given a collection of gene trees, this problem seeks a species tree that implies the minimum reconciliation cost, i.e., the fewest number of evolutionary events that can explain discordance in the gene phylogenies.
One of the most widely studied variants of the GTP problems is the gene duplication (GD) problem, in which the reconciliation cost is based on gene duplication events. The GD problem is W[2]-hard when parameterized by the number of gene duplications events and hard to approximate better than a logarithmic factor [5]. One way to cope with this intractability in practice is using heuristics [6, 7]. Although these heuristics do not guarantee optimal solutions or any non-trivial theoretical bound, in many cases they appear to have produced credible estimates [8–11]. However, the lack of performance guarantees makes the pursuit of exact solutions for the GD problem desirable.
Exact solutions can be found by exhaustive search for every NP-complete problem, but run times typically become prohibitively large for even rather small sized instances. However, exact algorithms that are substantially faster than exhaustive search have been progressively developed (e.g. [12, 13]). Unfortunately, little work has focused on such algorithms for the GD problem [14]. Here, we describe an ILP formulation solving the GD problem exactly and demonstrate its performance on both simulated and empirical data sets.
Related work
Exact solutions to the GD problem were obtained by exhaustively searching all possible species trees in data sets with up to 8 taxa [15, 16]. More recently, a branch-and-bound algorithm to identify exact solutions for the GD problem was introduced [14]. This algorithm was applied to a data-set consisting of 1, 111 gene trees with 29-taxa, but it did not run to completion. However, the branch-and-bound algorithm was able to solve this instance on reduced search spaces that resulted from providing some of the relationships in the species tree. Although constraining the search space for a species tree can help solving difficult instances of the GD problem, there are no theoretical guarantees to support this approach.
ILP formulations have provided an effective strategy to solve moderately sized instances of several NP-hard phylogenetic problems (e.g. [17–22]). Most similar to the GD problem, ILP formulations have been introduced for the deep coalescence problem, the variant of the GTP problem in which the reconciliation cost is based on the deep coalescence events [23]. These formulations solved instances with up to 8 taxa. However, perhaps due to the difficulty of directly expressing gene duplications in terms of linear equations, there have been no ILP formulations for the DP problem.
Our contributions
We introduce a novel approach to solve the GD problem exactly by describing the first ILP formulation for this problem. This solution is made possible by revealing an alternate description of the GD problem, called the triple inconsistency problem, which expresses gene duplications in terms of rooted triples. Rooted triples are rooted full binary trees with three leaves, and are the smallest unit of phylogenetic information. They, together with an equivalent presentation of species trees through cluster hierarchies, provide the fundamental elements of our ILP solution.
We demonstrate that our ILP formulation can solve non-trivial instances with up to 14 taxa and 1,000 gene trees. This greatly improves upon the largest (unconstrained) instances of the GD problem that have been solved exactly to date. Finally, we use the ILP formulation to address the relationships among the major seed plant lineages.Our ILP formulation was able to solve the GD problem exactly for a 12-taxon data set using 6,084 gene trees.
Methods
Preliminaries
Basic definitions
A rooted tree T is a connected and acyclic graph consisting of a vertex set V(T), an edge set E(T), and that has exactly one distinguished vertex called root, which we denote by Rt(T). Let T be a rooted tree. We define ≤_{ T } to be the partial order on V(T), where u ≤_{ T } v if v is a vertex on the path between Rt(T) and u. Moreover, we write u <>_{ T } v if neither u ≤_{ T } v nor v ≤_{ T } u is true. The set of minima under ≤_{ T } is denoted by L(T) and its elements are called leaves. We call u a child of v if u ≤ v and {u,v} ∈ V(E). The set of all children of v is denoted by Ch_{ T }(v). For a vertex v ∈ V(T) we denote by T(v) the subtree of T that consists of all vertices u ≤_{ T } v. The least common ancestor of a non-empty subset X ⊆ V(T), denoted as LCA_{ T }(X), is the unique smallest upper bound of X under ≤_{ T }. T is called full binary if every vertex has either two or zero children. Throughout this work, the term tree refers to a full and rooted binary tree.
Gene duplication (GD) problem
The terms species tree and gene tree refer to trees that represent the evolutionary history of a gene family or species respectively.
To compare a gene tree with a species tree, a mapping from each gene in the gene tree to the most recent species in the species tree that could have contained the gene is required.
Definition 1 (Mapping). Let G be a gene tree and S a species tree. A leaf-mapping from G to S is a function L_{ G,S } : L(G) → L(S). The extension M_{ G }_{,}_{ S }: V(G) → V(S) of the leaf-mapping L_{ G }_{,}_{ S } is the mapping defined by M_{ G,S }(u) := LCA_{ S }(L_{ G,S }(G(u)).
To simplify the exposition we shall assume that leaf-mappings are injections, and w.l.o.g. we identify the genes with the species from which they were sampled. After describing our ILP solution for identity leaf-mappings, we extend this formulation to cover non-injective leaf-mappings.
Definition 2 (Comparable). Let S be a species tree. A gene tree G is comparable to S, denoted by G ⊢ S, if L_{ G,S } exists. A set of gene trees is comparable to S, denoted by G ⊢ S, if G ⊢ S for each gene tree G ∈ G .
We shall adopt the following notation: we use S for a species tree, G for a set of gene trees that is comparable to S, and G for an gene tree in G .
Definition 3 (Duplication). A node g ∈ V(G) is a duplication (w.r.t. S) if M_{ G }_{,}_{ S }(g) ∈ M_{ G }_{,}_{ S }(Ch_{ G }(g)).
For consistency we follow the common practice to call what is stated above a definition, even though it is actually a theorem [24] that follows from the gene duplication model [2].
Definition 4 (Duplication cost). We define the following duplication costs:
1. Dup(G, S) := |{g ∈ V(G): g is a duplication}| is the duplication cost from G to S.
2. Dup(G , S) := ∑_{G∈}_{ G } Dup(G, S) is the duplication cost from G to S.
3. Dup(G) := min_{ G }_{⊢}_{ T } Dup(G, T) is the duplication cost of G.
Problem 1 (Gene-Duplication (GD)).
Instance: A set of gene trees G .
Find: The duplication cost Dup(G), and a species tree S* such that Dup(G, S*) = Dup(G).
The Triple-Inconsistency problem and its equivalence to the GD problem
A rooted triple is a tree with three leaves. The rooted triple with leaves x, y, and z is denoted xy|z if the path between x and y does not intersect with the path between z and the root. A rooted triple is displayed by a tree T if LCA_{ T }(x, y) ≤_{ T } LCA_{ T }(x, z) (= LCA_{ T }(y, z)). The set of rooted triples xy|z displayed by tree T that are rooted at vertex u ∈ V(T), (i.e., u = LCA_{ T }({x, y, z})) is denoted by Trip_{ T }(v), and the set of all triples displayed by T is denoted by Trip(T).
Definition 5 (T(riple)-inconsistency). A rooted triple t ∈ Trip(G) is said to be inconsistent with S if t ∉ Trip(S). A vertex v ∈ V(G) is called t(riple)-inconsistent with S if there is a rooted triple in Trip_{G}(v) that is inconsistent with S.
Definition 6 (T-inconsistency cost). We define the following t-inconsistency costs:
1. Tin(G,S) := |{v ∈ V(G): v is t-inconsistent with S}| is the t-inconsistency cost from G to S.
2. Tin(G, S) := ∑_{ G }_{∈}_{ G } Tin(G, S) is the t-inconsistency cost from G to S.
3. Tin(G) := min_{ G }_{⊢T} Tin(G,T) is the t-inconsistency cost of G .
Problem 2 (T(riple)-inconsistency).
Instance: A set of gene trees G .
Find: The t-inconsistency cost Tin(G), and species tree S* such that Tin(G, S*) = Tin(G).
Theorem 1 (Equivalence between duplication and t-inconsistency). Let u ∈ (G). Then u is a duplication w.r.t. S if and only if u is t-inconsistent with S.
Proof. Let x := M_{ G,S }(u).
Suppose u is a not a duplication. Let ab|c ∈ Trip_{G}(u). We will show that ab|c ∈ Trip(S). By the definition of ab|c ∈ Trip_{G}(u) we know that LCA_{G}({α, b,c}) = u, and together with our assumption that G is fully binary it follows that u has two children v and w, where w.l.o.g. a, b ∈ L(G(v)) and c ∈ L(G(w)). Let v′ := M_{ G,S }(v) and w′ := M_{ G,S }(W). From a, b ∈ L(G(v)) and c ∈ L(G(w)) follows that a, b ∈ L(S(v′)) and c ∈ L(S(w′)) respectively. Now, since u is not a duplication we have v′ <>_{ S } w′. Otherwise, we would have w′ ≤_{ S } v′ or v′ ≤_{ S } w′ from which x = v′ or x = w′ would follow respectively; contradicting that v is not a duplication. Hence, from v′ <>_{ S } w′ and a, b ∈ L(S(v′)) and c ∈ L(S(w′)) follows that ab|c ∈ Trip(S).
Suppose u is a duplication, and thus we have x = M_{ G,S }(v) for a child v ∈ Ch(u). So u is not a leaf in G, and since G is fully binary it follows that there are two distinct vertices a, b ∈ L(G(u)) such that LCA_{ S }({a,b}) = x. Therefore, x has two children y and z such that a ≤_{ S } y and b ≤_{ S } z. Now we distinguish different cases for the vertices a and b based on their possible order relation to the children of u. Since G is fully binary and v is a child of u, there exists a child w ∈ Ch(u) where w ≠ v. Now, we have the following cases.
Case 1: a ≤_{ G } v, b ≤_{ G } v: Let c ≤_{ G } w. Then ab|c ∈ Trip_{ G }(u). Further c ≤_{ S } y or c ≤_{ S }z and with a ≤_{ S }y and b ≤_{ S } z, it follows that either ac|b ∈ Trip_{s}(x) or bc|a ∈ Trip_{ S }(x). Hence, u is t-inconsistent as desired.
Case 2: a ≤_{ G } v, b ≤_{ G } w: We know that x has two children y and z and that M_{ G,S }(v) = x. Therefore there exist c ≤_{ S } y and d ≤_{ S } z such that LCA_{s}(c, d) = M(v) where c,d ∈ L(G(v)). From the order relations a ≤_{ S } y, d ≤_{ S } z and d ≤_{ G } v, b ≤_{ G } w and a ≤_{ G } v, b ≤_{ G } w, it follows that a, b and d are pairwise different. Therefore the rooted triples ad|b ∈ Trip_{ G }(u) and bd|a ∈ Trip_{ S }(x) are well defined, from which follows that the vertex u is t-inconsistent.
Case 3: a ≤_{ G } W, b ≤_{ G } w or b ≤_{ G } v, a ≤_{ G } w: Similarly to the previous cases it follows that u is t-inconsistent
From Theorem 1, the next corollary follows.
Corollary 1 (Equivalence between the GD problem and the T-Inconsistency problem). The t-inconsistency problem is a mathematical equivalent formulation of the duplication problem (i.e. Dup(G,S) = Tin(G,S)).
An ILP solution for the T-Inconsistency problem
Table 1
Notation | Definition |
---|---|
M(i, j) | Taxon-cluster representation of (the) species tree: M(i, j) = 1 iff taxon i is in the cluster j. Additional constraints on M require the cluster set to form a binary hierarchy (tree). |
C(p, q, xy) | Compatibility: C(p, q, xy) = 1 exactly if the cluster pair (p, q) has the gamete xy ∈ {01, 10, 11}. |
T(a, b, c, xyz) | Rooted triple: T(a, b, c, xyz) = 1 exactly if the rooted triple with leaf set {a, b, c} and topology xyz is displayed in M. Topologies for xyz are 011, 101, and 110 and refer to the rooted triples bc|a, ac|b and ab|c respectively. |
D(g) | t-inconsistency: D(g) = 1 if the gene vertex g is t-inconsistent w.r.t. a tree represented by matrix M. |
Let X := ∪_{ G }_{∈}_{ G } L(G) be the taxon set, n := |X|, m := |∪_{ G }_{∈}_{ G }Trip(G)|, and k := | G |. It follows that Σ_{ G }_{∈}_{ G }|G| = O(kn).
Formulating candidate species trees in terms of cluster hierarchies
Here we formulate constraints that describe all species trees that are possible candidates for solving the t-inconsistency problem, that is, all trees to which the given gene tree set G is compatible. Based on our assumption that the leaf label function is the identity function, these are all trees with the leaf set X. Our ILP formulation is based on an alternative way of describing trees by specifying their clusters through a hierarchy of subsets of X.
Definition 7 (Clusters). Let T be a tree. For each vertex v ∈ V(T) we define the cluster at v as {u ∈ L(T) : u ≤_{ T } v}, i.e., L(T(u)). We shall denote the set of all clusters of T by H(T).
Definition 8 ((Full) Binary hierarchy). Let F be a finite set. We call a set H of non-empty subsets of F a (full) binary hierarchy on F if the following properties are satisfied:
1. Trivial set property: F ∈ H and {v} ∈ H for each v ∈ F
2. Compatibility property: every pair of sets A and B in H is compatible; that is A ∩ B ∈ {A, B, Ø}.
3. Cardinality property: | H |= 2|F| –1
A hierarchy is defined as a binary hierarchy without requiring the cardinality property. There is a well-known and fundamental equivalence between hierarchies and trees that are not necessarily binary (e.g. [25]). The next result follows from this equivalence and the fact that a binary tree over l leaves has exactly 2l – 1 clusters.
Theorem 2 (Equivalence between binary hierarchies and binary trees). Let H be a set of non-empty subsets of a set F. Then there is a binary tree T such that H = H(F) if and only if H is a binary hierarchy on F.
Since trees and binary hierarchies are equivalent, we use these terms interchangeably from now on. Now we formulate constraints that describe the hierarchies on X using the binary matrix presentation.
Binary matrix. We describe 2n – 1 subsets of a hierarchy on X using a binary matrix M with a row for each species in X and 2n – 1 columns, where each column p represents the set {a ∈ X: M(a,p) = 1}.
- 2
≤ Σ_{ a }_{∈}_{ X }M(a, p) ≤ (n – 1).
Compatibility. Incompatibility can be tested directly by using the three-gamete condition (e.g., [26]). An incompatibility occurs for two columns p and q in M if and only if there exist three rows a, b and c in M that contain the gametes (0,1), (1,0), and (1,1) in p and q respectively (i.e. (M(a,p),M(a, q)) = (0,1), (M(b, p),M(b,q)) = (1,0), and (M(c, p),M(c, q)) = (1,1)). To identify if a certain gamete (x, y) ∈ {(0,1), (1,0), (1,1)} exists for p and q, we define a set of binary variables C(p, q, xy) under the following constraints over all rows a in M.
C(p, q, 01) ≥ –M(a, p) + M(a, q),
C(p, q, 10) ≥ M(a, p) – M(a, q),
C(p, q, 11) ≥ M(a, p) + M(a, q) – 1.
These constraints capture that C(p, q, xy) = 1 as long as M(a, p) = x and M(a, q) = y is satisfied for a gamete (x, y) in a certain row a in M. However, the reverse condition does not necessarily hold true without adding further constraints. To guarantee that clusters p, q are compatible, we require the following constraints
C(p, q, 01) + C(p, q, 11) + C(p, q, 10) = 2.
Number of variables and constraints. There are O(n^{ 2 }) variables for the matrix M, and O(n^{2}) variables of the type C(p, q, xy). O(n) constraints are needed to exclude trivial sets and to guarantee uniqueness, and O(n^{ 3 }) constraints guarantee compatibility. In summary, there are O(n^{2}) variables and O(n^{3}) constraints to describe all candidates for the species tree.
Formulating the T-lnconsistency problem. To formulate the t-inconsistency problem, we first describe variables T(a, b, c, xyz) that detect whether a rooted triple is displayed by the tree presented by M. Then we describe variables D(g) that detect if g is t-inconsistent by using the variables T(a, b, c, xyz). Finally, we formulate the objective of the t-inconsistency problem based on the variables D(g).
Variables T ( a, b, c, xyz ) . We describe the binary variables T(a, b, c, xyz) that are 1 exactly if a rooted triple over the leaf set {a, b, c} with topology (x, y, z) ∈ {(0,1,1), (1,0,1), (1,1,0)} is displayed by the tree that is presented by M. The parameters a, b, c are rows in M, and the settings 011, 101, and 110 of (x, y, z) refer to the rooted triples bc|a, ac|b and ab|c respectively. For each column p in M, we introduce the following constraints.
T(a, b, c, 011) ≥ –M(a, p) + M(b, p) + M(c, p) – 1 ;
T(a, b, c, 101) ≥ M(a, p) – M(b, p) + M(c, p) – 1 ;
T(a, b, c, 110) ≥ M(a, p) + M(b, p) - M(c, p) – 1 ;
T(a, b, c, 011) + T(a, b, c, 101) + T(a, b, c, 110) = 1 ,
since a rooted triple is uniquely resolved in a tree.
Variables T(a, b, c, 011), T(a, b, c, 101), and T(a, b, c, 110) are constructed for every triple {a, b, c} for which a rooted triple is displayed by a gene tree in . Thus, there are O(m) variables of this type. For each variable we have O(n) constraints, which results in O(nm) constraints overall.
Variables D ( g ) . We express the t-inconsistency of each vertex g ∈ V(G) where G ∈ G by the binary variable D(g). The variable is 1 if g is t-inconsistent with the tree described by matrix M, given the following constraints
D(g) ≥ 1 – T(a, b, c, xyz) ,
where the rooted triple over the leaf set {a, b, c} and topology xyz is an element in Trip_{ G }(g).
Variables D(g) are constructed for each internal vertex of a gene tree in G , which results in O(kn) variables. Intuitively, a constraint is constructed for each rooted triple that is displayed by a gene tree in G, which yields O(km) constraints. However, the following observation reduces the number of such constraints to O(kn^{2}).
Let u ∈ V(G) such that Trip_{G}(u) ≠ Ø, {v,w} = Ch(u), a ∈ L(G(u)) and b ∈ L(G(v)). A rooted triple xy|z is in Trip_{G}(u) if and only if all ax|b, ay|b, and bz|a are in Trip_{G}(u). Therefore, instead of enumerating all rooted triples in Trip_{G}(u) (which sums up to O(m) in each gene tree G), we only need to enumerate a number of O(n) rooted triples to represent Trip_{ G }(u) while detecting if u is t-inconsistent (hence O(kn^{ 2 }) constraints over all).
T-lnconsistency objective. This objective is expressed by the following expression.
min Σ_{g∈}_{ V }_{(}_{ G }_{)}D(g).
Once the optimal objective cost is found, a unique tree corresponding to the cost can be constructed from M. It is worth noting that an instance of unique optimal tree does not ensure an unique optimal solution to the corresponding ILP due to relaxed constraints for variables C. Although this can be addressed by adding additional constraints, the correctness of the objective value and the resulting tree is not affected.
Number of variables and constraints. In summary, there are O(n^{2} + m + kn) variables, and the number of constraints is O(n^{3}+ mn + kn^{2}).
Handling non-injective leaf mappings
A leaf mapping L_{ G,S } is non-injective if and only if there is a vertex u ∈ V(G) with distinct children v and w such that L_{ G,S }(L(G(v)))ØL_{ G,S }(L(G(w))) ≠ Ø; and if the latter holds true, it follows that u is a duplication. Therefore, it can be determined if u is a gene-duplication regardless of the topology of S. By pre-processing all such determined duplication vertices, the leaf-mapping over the remaining internal vertices of G can be made injective. Hence, the existing ILP formulation solves input gene trees with non-injective leaf mappings. Since the input gene tree size can be arbitrary, under the non-injective leaf mapping assumption, the ILP formulation has O(n^{2}+ m + l) variables and O(n^{3}+ mn + In) constraints where Σ_{ G }_{∈}_{ G }|G| = l.
Generating optimal species trees
The species tree corresponding to a feasible solution of an ILP instance can be constructed in O(n^{2}) time [27]. Furthermore, a gene node g is identified as a duplication if and only if D(g) = l.
Implementation
We implemented an ILP generator in Python that, given a set of gene trees, outputs the ILP described in the preceding section. We tested our formulation with both simulated and empirical gene tree data sets (described below). All analyses were on a GNU/Linux based PC with an Intel Core2 Quad 2.4 GHz CPU. We choose Gurobi 2.0 [28] to solve the ILP directly and CPLEX 12.1 [29] to enumerate optimal solutions when necessary.
Simulation experiments
We first evaluated the performance of our ILP solution with simulated gene tree data sets. Our simulation protocol included the following steps: (1) a species tree S of n taxa was randomly generated as the template of a gene tree; (2) a depth-first-search walk starting from Rt(S) simulated at most one evolutionary event at each vertex based on given probabilities for each event. These events could be a duplication (duplicating the whole current subtree) or a loss (cutting the current subtree). If there is neither a duplication nor a loss, the process proceeds to the next vertex. We used the same species tree to generate k gene trees.
In our simulation experiments, we used a duplication rate of 0.25 duplications per gene at each spe-ciation vertex and a loss rate of 0.3. These events produced a similar tree size distribution and optimal duplication cost to the gene trees used by Sanderson and McMahon [16]. We varied the number of taxa in the species tree from 6 to 14 and the number of input gene trees from 10 to 1000. We performed 10 simulation replicates for each different combination of species and gene tree number. For each simulated data set, we also compared the ILP score to results from DupTree [7], a fast hill-climbing heuristic implementation for the problem, to determine if the heuristic finds the optimal solution.
Seed plant analysis
Next, we tested the ability of the ILP formulation to solve the seed plant phylogeny problem using a large-scale genomic data set. First, to build the gene trees, amino acid alignments for gene families were selected from Phytome v. 2, an online comparative genomics database based on publicly available sequence data from 136 plant species [30]. To ensure positional homology throughout the alignments, columns and sequences of questionable certainty were masked using default settings of the program REAP [30, 31].We sampled sequences from the nine gymnosperm taxa represented in Phytome with the most data, including cycad taxon Cycas rumphii, Gnetales taxa Gnetum gnemon and Welwitschia mirabilis, and, from the conifers, Cryptomeria japonica from Cupres-saceae, and Pseudotsuga menziesii, Picea glauca, Picea sitchensis, Pinus pinaster, and Pinus taeda from Pinaceae. We also sampled sequences from two representative angiosperm taxa, Arabidopsis thaliana and Oryza sativa, and the non-seed plant, Physcomitrella patens.
We selected all the 6,084 masked amino acid alignments from gene families in Phytome that had at least 4 sequences and had sequences from at least 3 of the selected taxa. All species were found in at least 376 gene families. To build the gene trees, we performed ML phylogenetic analyses on each of the gene alignments using RAxML-VI-HPC version 2.2.3 [32]. The ML analyses used the JTT amino acid substitution model [33] with rate variation among sites (the “PROTMIX” model; see [32]). The trees were then rooted using mid-point rooting, as implemented in the Phylip program retree [34]. We applied the ILP formulation to solve the GD problem using all 6, 084 gene trees.
Results and discussion
Simulations
Table 2
n = 6 | n = 8 | n = 10 | n = 12 | n = 14 | ||||||
---|---|---|---|---|---|---|---|---|---|---|
k | time | Dup | time | Dup | time | Dup | time | Dup | time | Dup |
10 | 0.06 | 34.80 | 0.34 | 49.70 | 22.98 | 60.10 | 200.53 | 68.80 | 12597.21 | 78.40 |
50 | 0.03 | 189.50 | 1.26 | 265.00 | 8.74 | 280.00 | 159.26 | 346.40 | 2953.62 | 393.10 |
100 | 0.06 | 382.80 | 0.63 | 523.30 | 9.64 | 598.50 | 117.38 | 701.60 | 2191.65 | 825.70 |
200 | 0.05 | 788.20 | 0.54 | 994.90 | 11.03 | 1217.30 | 168.85 | 1372.50 | 2709.91 | 1627.70 |
500 | 0.25 | 1910.30 | 0.79 | 2458.60 | 13.92 | 2987.00 | 220.17 | 3678.80 | 4270.05 | 4001.70 |
1000 | 0.57 | 3842.60 | 0.96 | 5283.10 | 23.54 | 6140.90 | 330.34 | 7026.40 | 5014.61 | 8258.80 |
Seed plant analysis
The relationships among the major lineages of seed plants has long been a major question in plant systematics, especially with regard to the position of Gnetales, a clade of three genera (Gnetum, Ephedra, and Welwitschia) that lack obvious morphological links to other extant seed plants (e.g.,, [35, 36, 39]). Cladistic analyses of morphological characters generally have placed Gnetales sister to the angiosperms, or flowering plants [36, 40–44]; however, early analyses of molecular characters rarely supported this placement [35, 37, 39, 45]. Most recently, maximum likelihood (ML) and maximum parsimony (MP) analysis of 15-17 plastid loci placed Gnetales sister to the other seed plants [46]. However, a loss of plastid ndh genes appears to link Gne-tales with Pinaceae [47]. An MP analysis of EST sequences from 43 nuclear genes similarly linked Gnetales with the conifers [48]. Yet later MP and ML analyses of EST sequences from over 1,200 nuclear loci placed Gnetales sister to the other gym-nosperms [49]. All of these molecular analyses of the seed plant phylogeny have been limited to puta-tively orthologous genes. However, the GD problem provides a way to incorporate large gene families into the phylogenetic inference of seed plants.
Conclusions
Our ILP formulation provides exact solutions to the largest instances of the GD problem analyzed to date. Thus, it can provide a large-scale genomic perspective on important phylogenetic questions that previously could only be addressed by heuristics. Furthermore, our simulation experiments demonstrate that these heuristic estimates can be misled with as few as 10 taxa. Even when heuristics identify an optimal solution they cannot, unlike ILP, determine if the solution is unique. In future research the ILP implementation will be useful, not only for solving empirical data sets, but for assessing the performance of different heuristics by comparing their estimates to the exact ILP solution. Ultimately, it also will be useful to expand the scale of solvable instances beyond 14 taxa. While this challenge may be addressed by improved ILP formulations, investigations into other algorithm concepts might also be effective (e.g., [14, 23]).
Authors contributions
WCC was responsible for developing the solution, running experiments, and writing of the manuscript. JGB performed the experimental evaluation and the analysis of the results, and contributed to the writing of the manuscript. OE and DFB supervised the project and contributed to the writing of the paper. All authors read and approved the final manuscript.
Declarations
Acknowledgements
We thank Mukul Bansal for discussions and reviewers for helpful comments. This work was supported in part by NSF awards #0830012 and #1017189.
This article has been published as part of BMC Bioinformatics Volume 12 Supplement 1, 2011: Selected articles from the Ninth Asia Pacific Bioinformatics Conference (APBC 2011). The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2105/12?issue=S1.
Authors’ Affiliations
References
- Maddison WP: Gene trees in species trees. Syst. Biol 1997, 46: 523–536. 10.1093/sysbio/46.3.523View ArticleGoogle Scholar
- 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. Syst. Zool 1979, 28: 132–163. 10.2307/2412519View ArticleGoogle Scholar
- Guigó R, Muchnik I, Smith TF: Reconstruction of Ancient Molecular Phylogeny. Mol. Phylogenet. Evol. 1996, 6(2):189–213.View ArticlePubMedGoogle Scholar
- Slowinski JB, Knight A, Rooney AP: Inferring Species Trees from Gene Trees: A Phylogenetic Analysis of the Elapidae (Serpentes) Based on the Amino Acid Sequences of Venom Proteins. Mol. Phylogenet. Evol. 1997, 8(3):349–362. 10.1006/mpev.1997.0434View ArticlePubMedGoogle Scholar
- Bansal MS, Shamir R: A Note on the Fixed Parameter Tractability of the Gene-Duplication Problem. IEEE/ACM Trans. Comput. Biol. Bioinf. 2010.Google Scholar
- Bansal MS, Burleigh JG, Eulenstein O, Wehe A: Heuristics for the Gene-Duplication Problem: A Θ(n) Speed-Up for the Local Search. RECOMB, Volume 4453 of LNCS 2007, 238–252.Google Scholar
- Wehe A, Bansal MS, Burleigh JG, Eulenstein O: Dup-Tree: a program for large-scale phylogenetic analyses using gene tree parsimony. Bioinformatics 2008, 24(13):1540–1541. 10.1093/bioinformatics/btn230View ArticlePubMedGoogle Scholar
- Page RDM: Extracting Species Trees From Complex Gene Trees: Reconciled Trees And Vertebrate Phylogeny. Mol. Phylogenet. Evol. 2000, 14: 89–106. 10.1006/mpev.1999.0676View ArticlePubMedGoogle Scholar
- Cotton JA, Page RDM: Going Nuclear: Gene Family Evolution And Vertebrate Phylogeny Reconciled. Proc Biol Sci 2002, 269: 1555–1561. 10.1098/rspb.2002.2074PubMed CentralView ArticlePubMedGoogle Scholar
- Martin AP, Burg TM: Perils of Paralogy: Using HSP70 Genes for Inferring Organismal Phylogenies. Syst. Biol. 2002, 51(4):570–587. 10.1080/10635150290069995View ArticlePubMedGoogle Scholar
- McGowen MR, Clark C, Gatesy J: The Vestigial Olfactory Receptor Subgenome of Odontocete Whales: Phylogenetic Congruence between Gene-Tree Reconciliation and Supermatrix Methods. Syst. Biol. 2008, 57(4):574–590. 10.1080/10635150802304787View ArticlePubMedGoogle Scholar
- Applegate DL, Bixby RE, Chvatal V, Cook WJ: The Traveling Salesman Problem: A Computational Study (Princeton Series in Applied Mathematics). Princeton University Press; 2007.Google Scholar
- Woeginger GJ: Exact algorithms for NP-hard problems: A survey. Combinatorial OptimizationÂ–Eureka, You Shrink! 2003, 2570/2003: 185–207.View ArticleGoogle Scholar
- Doyon JP, Chauve C: Branch-and-Bound Approach for Parsimonious Inference of a Species Tree From a Set of Gene Family Trees. In Tech. rep.. LIRMM; 2010.Google Scholar
- Burleigh JG, Bansal MS, Eulenstein O, Vision TJ: Inferring Species Trees From Gene Duplication Episodes. Proc. ACM-BCB 2010, 198–203.Google Scholar
- Sanderson MJ, McMahon M: Inferring angiosperm phylogeny from EST data with widespread gene duplication. BMC Evol. Biol. 2007, 7(Suppl 1):S3. 10.1186/1471-2148-7-S1-S3PubMed CentralView ArticlePubMedGoogle Scholar
- Brown DG, Harrower IM: Integer Programming Approaches to Haplotype Inference by Pure Parsimony. IEEE/ACM Trans. Comput. Biol. Bioinf. 2006, 3(2):141–154. 10.1109/TCBB.2006.24View ArticleGoogle Scholar
- Dong J, Fernández-Baca D, McMorris FR: Constructing majority-rule supertrees. Algorithms for Molecular Biology 2010, 5: 2. 10.1186/1748-7188-5-2PubMed CentralView ArticlePubMedGoogle Scholar
- Gusfield D: The Multi-State Perfect Phylogeny Problem with Missing and Removable Data: Solutions via Integer-Programming and Chordal Graph Theory. RECOMB 2009, 236–252.Google Scholar
- Gusfield D, Frid Y, Brown DG: Integer Programming Formulations and Computations Solving Phylogenetic and Population Genetic Problems with Missing or Genotypic Data. COCOON 2007, 51–64.Google Scholar
- Sridhar S, Lam F, Blelloch GE, Ravi R, Schwartz R: Efficiently finding the most parsimonious phylogenetic tree via linear programming. Int. J. Bioinf. Res. Appl. 2007, 4463: 37–48. full_textView ArticleGoogle Scholar
- Chimani M, Rahmann S, Sebastian B: Exact ILP Solutions for Phylogenetic Minimum Flip Problems. Proc. ACM BCB 2010, 147–153.Google Scholar
- Than C, Nakhleh L: Species Tree Inference by Minimizing Deep Coalescences. PLoS Comput. Biol. 2009, 5(9):e1000501. 10.1371/journal.pcbi.1000501PubMed CentralView ArticlePubMedGoogle Scholar
- Eulenstein O: Vorhersage von Genduplikationen und deren Entwicklung in der Evolution. In PhD dissertation. University of Bonn; 1998.Google Scholar
- Semple C, Steel MA: Phylogenetics. Oxford University Press; 2003.Google Scholar
- Gusfield D: Algorithms on Strings, Trees and Sequences: Computer Science and Computational Biology. Cambridge University Press; 1997.View ArticleGoogle Scholar
- Gusfield D: Efficient algorithms for inferring evolutionary trees. Networks 1991, 21: 19–28. 10.1002/net.3230210104View ArticleGoogle Scholar
- Gurobi Optimization, Inc: Gurobi Optimization 2.0.2.2010. [http://www.gurobi.com/]Google Scholar
- IBM, Inc: IBM ILOG CPLEX 12.1.2009. [http://www.ibm.com/software/integration/optimization/cplex/]Google Scholar
- Hartmann S, Lu D, Phillips J, Vision TJ: Phytome: a platform for plant comparative genomics. Nucleic Acids Res 2006, 34(Database issue):D724-D730. 10.1093/nar/gkj045PubMed CentralView ArticlePubMedGoogle Scholar
- Hartmann S, Vision TJ: Using ESTs for phylogenomics: Can one accurately infer a phylogenetic tree from a gappy alignment? BMC Evol. Biol. 2008, 8: 95. 10.1186/1471-2148-8-95PubMed CentralView ArticlePubMedGoogle Scholar
- Stamatakis A: RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics 2006, 22(21):2688–2690. 10.1093/bioinformatics/btl446View ArticlePubMedGoogle Scholar
- Jones DT, Taylor WR, Thornton JM: The rapid generation of mutation data matrices from protein sequences. Comput. Appl. Biosci. 1992, 8(3):275–282.PubMedGoogle Scholar
- Felsenstein J: PHYLIP (Phylogeny Inference Package) version 3.6. Distributed by the author 2005.Google Scholar
- Burleigh JG, Mathews S: Phylogenetic signal in nucleotide data from seed plants: Implications for resolving the seed plant tree of life. Am. J. Bot. 2004, 91(10):1599–1613. 10.3732/ajb.91.10.1599View ArticlePubMedGoogle Scholar
- Donoghue MJ, Doyle JA: Seed plant phylogeny: Demise of the anthophyte hypothesis? Current Biology 2000, 10(3):R106-R109. 10.1016/S0960-9822(00)00304-3View ArticlePubMedGoogle Scholar
- Magallón S, Sanderson MJ: Relationships among Seed Plants Inferred from Highly Conserved Genes: Sorting Conflicting Phylogenetic Signals among Ancient Lineages. Am. J. Bot. 2002, 89(12):1991–2006.View ArticlePubMedGoogle Scholar
- Mathews S: Phylogenetic relationships among seed plants: Persistent questions and the limits of molecular data. Am. J. Bot. 2009, 96: 228–236. 10.3732/ajb.0800178View ArticlePubMedGoogle Scholar
- Soltis DE, Soltis PS, Zanis MJ: Phylogeny of Seed Plants Based on Evidence from Eight Genes. Am. J. Bot. 2002, 89(10):1670–1681. 10.3732/ajb.89.10.1670View ArticlePubMedGoogle Scholar
- Crane PR: Phylogenetic Analysis of Seed Plants and the Origin of Angiosperms. Annals of the Missouri Botanical Garden 1985, 72: 716–793. 10.2307/2399221View ArticleGoogle Scholar
- Doyle JA: Seed Ferns and the Origin of Angiosperms. The Journal of the Torrey Botanical Society 2006, 133: 169–209. 10.3159/1095-5674(2006)133[169:SFATOO]2.0.CO;2View ArticleGoogle Scholar
- Doyle JA, Donoghue MJ: Seed plant phylogeny and the origin of angiosperms: An experimental cladistic approach. The Botanical Review 1986, 52(4):321–431. 10.1007/BF02861082View ArticleGoogle Scholar
- Hilton J, Bateman RM: Pteridosperms are the backbone of seed-plant phylogeny. The Journal of the Torrey Botanical Society 2006, 133: 119–168. 10.3159/1095-5674(2006)133[119:PATBOS]2.0.CO;2View ArticleGoogle Scholar
- Nixon KC, Crepet WL, Stevenson DW, Friis EM: A Reevaluation of Seed Plant Phylogeny. Annals of the Missouri Botanical Garden 1994, 81(3):484–533. 10.2307/2399901View ArticleGoogle Scholar
- Rydin C, Kallersjo M, Friis EM: Seed Plant Relationships and the Systematic Position of Gnetales Based on Nuclear and Chloroplast DNA: Conflicting Data, Rooting Problems, and the Monophyly of Conifers. Int. J. Plant Sci. 2002, 163(2):197–214. 10.1086/338321View ArticleGoogle Scholar
- Rai HS, Reeves PA, Peakall R, Olmstead RG, Graham SW: Inference of higher-order conifer relationships from a multi-locus plastid data set. Botany 2008, 86: 658–669. 10.1139/B08-062View ArticleGoogle Scholar
- Braukmann TWA, Kuzmina M, Stefanovic S: Loss of all plastid ndh genes in Gnetales and conifers: extent and evolutionary significance for the seed plant phylogeny. Current Genetics 2009, 55(3):323–337. 10.1007/s00294-009-0249-7View ArticlePubMedGoogle Scholar
- de La Torre-Bárcena JE, Egan M, Katari MS, Brenner ED, Stevenson DW, Coruzzi GM, DeSalle R: ESTimating plant phylogeny: lessons from partitioning. BMC Evol. Biol. 2006, 6: 48.View ArticleGoogle Scholar
- de La Torre-Bárcena JE, Kolokotronis SO, Lee EK, Stevenson DW, Brenner ED, Katari MS, Coruzzi GM, DeSalle R: The Impact of Outgroup Choice and Missing Data on Major Seed Plant Phylogenetics Using Genome-Wide EST Data. PLoS ONE 2009, 4(6):e5764.PubMed CentralView ArticlePubMedGoogle Scholar
- Burleigh JG, Mathews S: Assessing systematic error in the inference of seed plant phylogeny. Int. J. Plant Sci. 2007, 168(2):125–135. 10.1086/509588View ArticleGoogle Scholar
- Wu CS, Wang YN, Liu SM, Chaw SM: Chloroplast Genome (cpDNA) of Cycas taitungensis and 56 Cp Protein-coding Genes of Gnetum parvifolium: Insights into CpDNA Evolution and Phylogeny of Extant Seed Plants. Mol. Biol. Evol. 2007, 24: 1366–1379. 10.1093/molbev/msm059View ArticlePubMedGoogle Scholar
- Holland BR, Penny D, Hendy MD: Outgroup Misplacement and Phylogenetic Inaccuracy under a Molecular Clock: A Simulation Study. Syst. Biol. 2003, 52(2):229–238. 10.1080/10635150390192771View ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.