Here, we propose a reconciliation model based on DTL parsimony that distinguishes between regions of the species tree where ILS is likely, and those where only gene duplication and transfer need be considered. These differences are specified using a non-binary species tree: at binary nodes, we assume that ILS is so rare that incongruence is always evidence of gene duplication or transfer. At polytomies, ILS is considered, and gene duplication and transfer are invoked only if topological disagreement cannot be explained by ILS. This model can be invoked for both non-binary species trees and for binary species trees with short branches where ILS is suspected: even when the binary branching order of the species tree is known, the user can collapse edges in the species tree to indicate in which lineages ILS should be considered as an alternate hypothesis.
A key aspect of our model is that even when ILS is allowed, it is not possible to explain all incongruence in terms of ILS, even in a uniquely labeled gene tree. Let
g be a node in
TG and let
s ![[set membership]](/corehtml/pmc/pmcents/x2208.gif)
VS be the associated node in the species tree. We wish to determine whether the divergence at
g is consistent with a co-divergence at
s or whether it can only be explained by events that give rise to a new locus; i.e. duplication and transfer. If the branch point at
g arose through a co-divergence with
s, then each species lineage descending from
s should inherit at most one descendant of
g. The presence of more than one descendant of
g indicates that the divergence at
g must be due to acquisition of an additional locus by duplication or transfer. An operational test for detecting more than one descendant on a branch results from the observation that any branching pattern that is consistent with a binary resolution of the polytomy can be explained by lineage sorting.
For example, the gene tree in represents a valid, binary resolution of the species tree, consistent with ILS. The embedding of the gene tree in the species tree shows that each species tree lineage inherits exactly one descendant of x1 and at most one descendant of x2. Both x1 and x2 can be interpreted as deep coalescences. In contrast, there is no binary resolution of the species tree that corresponds to the gene tree in . The embedding of this gene tree requires two descendants of y2 in the lineage from e to f, a violation of model constraints. The only way to explain two descendants of y2 on the branch from e to f is by inferring a duplication () or a transfer ().
2.1 The DTLI algorithm
In our DTLI model, divergence in a gene tree arises through one of four events: duplication (
), transfer (
), speciation (
) and deep coalescence (
). The score of a reconciliation under this model is the weighted sum of the number of duplications (
), losses (
), and transfers (
):
where
δ,
λ and
τ, respectively, are the costs of duplication, loss and transfer. Speciation and deep coalescence represent co-divergence with binary nodes and polytomies, respectively, in the species tree and have zero cost. We refer to the cost of event
ε {
,
,
,
} as
κ (
ε).
A rooted, binary gene tree
TG; a rooted, arbitrary species tree
TS; a mapping
ML :
L(
VG) →
L(
VS) from contemporary genes to the species from which they were sampled and a set of permitted events are given as input. The reconciliation of
TG with
TS results in an annotated tree,
RGS = (
VG,
EG), in which every internal node,
g, is annotated with the species
s VS that contained gene
g, designated
M(
g), and the event that caused the divergence at
g, designated
. In addition, every
g VG \ {
ρG} is annotated with
, the genes lost on the edge from
P(
g) to
g. Each loss is labeled with the species in which the loss occurred. We say (
u,
v)
EG is a transfer edge if
and
M (
u)
![[not less, greater]](/corehtml/pmc/pmcents/x2278.gif)
S
M(
v) and define Λ(
RGS)
![[subset or is implied by]](/corehtml/pmc/pmcents/x2282.gif)
EG to be the set of transfer edges in
RGS. If (
u,
v)
Λ(
RGS), a transfer occurred from donor species
d =
M(
u) to recipient species
r =
M(
v).
Here, we present the DTLI event inference problem under the constraint that a deep coalescent is inferred at g if each lineage descending from M(g) inherits at most one descendant of g:
The DTLI event inference problem
- Input: A rooted non-binary species tree, TS; a rooted, binary gene tree, TG; the leaf mapping, ML.
- Output: All reconciliation histories RGS that minimize π and satisfy the model constraints.
Algorithms for the DTLI event model must address several issues that do not arise when only a subset of the events is considered: (1) there may be more than one combination of duplications, transfers and losses that gives rise to the same pattern of tree incongruence (i.e. there may be more than one optimal solution,
RGS). (2) The value of
M(
g) is not uniquely determined by the children of
g and multiple possible values of
M(
g) must be considered because transfers cause genes to jump to distant locations in the species tree. (3) An optimal reconciliation at the root may entail a suboptimal reconciliation at an internal node,
g. Inferring a more costly event at
g may change the values of
M(·) in nodes ancestral to
g such that the overall score is reduced. Therefore, the values of
M(
g) and
required for an optimal solution cannot be determined using only local information, and more than one optimal solution may result.
To accommodate these requirements, it is necessary to enumerate all possible assignments of
M(
g) and
, for each node
g VG. At each
g, the associated information is stored in two tables,
and
. For each candidate assignment
s VS, the score that minimizes the cost of reconciling
TG(
g) with
TS(
s), is stored in
. The associated events and other information needed to reconstruct the history at
g are stored in
.
Optimal reconciliations are calculated by a two-pass algorithm. The first pass (Algorithm 2.1.1) is a dynamic program that populates each
and
in a post-order traversal of
TG. It returns the optimal reconciliation score, the values of
M(
ρg) and
corresponding to that score and the number of optimal histories. The second pass (
Supplementary Algorithm S1.0.1) is a traceback algorithm that reads information from each
to construct an optimal solution. Each optimal history is generated by traversing, in pre-order of
TG, each unique path that leads to the optimal label(s) in
. Appropriate values of
M(
g) and
at each node
g are selected from
. Each candidate optimal history is then tested for temporal feasibility, as described in the next section. Only those histories that are temporally feasible are reported.
A key calculation in the dynamic program of firstPass is determination of the possible events at
g for a given candidate species assignment,
M(
g)=
s. These events, in turn, depend on
M(
c1) =
s1 and
M(
c2) =
s2, where
c1,
c2 C(
g). The basis for determining candidate events that are consistent with
s,
s1 and
s2 is the following observation: if a duplication occurred at
g, then the species that inherit the descendants of
c1 and the species that inherit the descendants of
c2 will not be disjoint.
We define a test, based on this observation, for distinguishing duplication from other events:
where
is the set of species that vertically inherit descendants of
P(
g). If
and
are disjoint, than one of the other three events (
,
or
) must have occurred. These events can be distinguished from one another using
,
M(
g) and
M(
c1) and
M(
c2), as seen in costCalc in Algorithm 2.1.1. Note that
Equation (2) is different from the standard least common ancestor (lca) test; however, when
M(
g)=
s is binary, the descendants of
s are partitioned into two sets, the left and right descendants of
s, if there is no duplication. Therefore,
Equation 2 is equivalent to lca reconciliation (
Vernot et al., 2008).
Because
only consists of elements that were vertically inherited, we must exclude transfer edges in the calculation. For this purpose, we define
the set of leaves of
TG(
g) that were acquired through HGT. Formally, we define
to be a mapping from
VG to sets of nodes in
VS, where
VS+ is the powerset of
VS.
is the set of children of
M(
P(
g)) such that
; otherwise,
One more piece of machinery is needed: to determine
, we must know the children of
M(
P(
g)), but we do not have that information until we visit
P(
g). Therefore, we define a similar set mapping,
, to aid in the calculation of
.
is the, set of children of
M(
g) that vertically inherit a descendant of
g. Formally, if
M(
g)
L(
TS),
; otherwise,
Algorithm 2.1.1 traverses
TG in post-order calling calcCost at each
g VG. The challenge in the DTLI model is to determine the sets of species that inherit the descendants of
c1 and
c2 when
M(
g) =
s is a polytomy; i.e. how to calculate
and
. When
s is binary, the descendants of
s are easily partitioned into two sets; when
s is a polytomy, all possible ways to partition the descendants must be considered. Each child of
g can be retained in any subset of the children of
s, ranging from size 1 to |
C(
s)| − 1. Our DTLI algorithm addresses this by considering all ways of partitioning
C(
s) into two non-empty subsets.
At each internal node
g, the algorithm assesses all possible values for
M(
g) and
by looping through all (
s1,
s2)
VS ×
VS and all
. Considering all power sets corresponds to considering all the ways to partition
C(
s1) and
C(
s2). The optimal event and child mapping under
s and
is determined by minimizing the cost of the candidate solution at
g:
where
, the number of losses on edge (
g,
ci), is calculated using the loss heuristic in (
Vernot et al., 2008). Note that for each
s, the local cost and history tables are also indexed by all possible values of
, which are in
C(
s)
+.
2.2 Temporal infeasibility
Because the donor and recipient species of any transfer must have coexisted, each transfer implies a temporal constraint. A reconciliation is temporally feasible if an ordering of species exists that satisfies the constraints of all inferred transfers. Because reconciliations inferred by Algorithm 2.1.1 are not guaranteed to be feasible, each candidate optimal solution is tested for feasibility post hoc.
To determine whether a reconciliation
RGS is temporally feasible, we construct a directed timing graph
Gt = (
Vt,
Et) that encodes all temporal constraints on species in
TS. Only species that are the donor,
d, or recipient,
r, of a transfer edge in Λ(
RGS) must be considered. Thus, the vertex set is defined as
Vt = {
v VS|
![[exists]](/corehtml/pmc/pmcents/x2203.gif)
(
g,
h)
Λ(
TG)
![[contains]](/corehtml/pmc/pmcents/x220B.gif)
v =
M(
g)
![[logical or]](/corehtml/pmc/pmcents/x2228.gif)
v =
M(
h)}.
The edges in
Et represent three types of temporal constraints:
- If species si is an ancestor of species sj in TS, then si predates sj: for every (si, sj) in Vt × Vt, add (si, sj) to Et if si ≥S
sj.
- Let (g, h) and (g′, h′) be transfers in Λ(RGS), such that g ≥G
g′. Then d = M(g) and r = M(h) must have occurred no later than both d′ = M(g′) and r′ = M(h′). We add (P(d), d′), (P(d), r′), (P(r), d′) and (P(r), r′) to Et .
- Given a transfer (g, h) Λ(RGS), species M(g) and M(h) must be contemporaneous. Furthermore, any species that predates M(g) must also predate M(h) and vice versa. For every (si, sj) Vt × Vt, add (si, sj) to Et if sk Vt such that si ≥S
sk and sk and sj are the donor and recipient, or vice versa, of some transfer (g, h) Λ(RGS).
We test each candidate optimal history for temporal feasibility by verifying that the associated timing graph
Gt is acyclic, using a modified topological sorting algorithm in Θ (|
Vt |+|
Et |) (
Cormen et al., 1990). Temporally infeasible histories are not reported. Note that it is not the case that if one optimal history is infeasible, all optimal histories are infeasible. Finding the optimal, temporally feasible reconciliation is NP-complete (
Tofigh et al., 2011); we leave the problem of obtaining an optimal, feasible solution when all candidate solutions have infeasible timing constraints for future work.
2.3 Complexity and running time
Our algorithm is fixed-parameter tractable with polynomial complexity when the size of the largest polytomy,
k*, is fixed. In practical data analyses,
k* is likely to be small. Recent genome-scale analyses of ILS have focused on species trees with
k* = 3 (
Ebersberger et al., 2007;
Pollard et al., 2006). In general, event inference will not yield informative results when the species tree is highly unresolved.
Theorem 2.1. Given a binary gene tree TG and a non-binary species tree TS, firstPass takes O(|VG|(|VS|+ nk 2k*)2(hS + k*)) time.
P
roof. firstPass visits each
g VG in post order. At each
g, costCalc is called once for every (
s1,
s2)
VS ×
VS and
, resulting in a total of
O(|
VG|(|
![[union or logical sum]](/corehtml/pmc/pmcents/x222A.gif)
sVs C(
s)
+|)
2) calls to costCalc. Because |
C(
s)
+| = 2
|C(s)| is
O(1) when
s is binary, |
sVs C(
s)
+| is bounded above by |
VS|–
nk +
nk2
k* and the number of calls to costCalc is
O(|
VG|(|
VS|+
nk2
k*)
2). We precalculate lca(
s1,
s2) and test whether
s1s2, for all species pairs, in
O(|
VS|
2) time. Therefore, the complexity of costCalc is dominated by the calculations of
for
l and
r,
and
. These values can be computed in
O(
hS),
O(log(
k*)) and
O(
k*) time, respectively. Thus, each call to costCalc has complexity
O(
hS +
k*). Once the post-order traversal is completed, we extract the minimum score in
, and all values of
M(
ρG) and
corresponding to that score. Since
, a linear search accomplishes this in
O(|
VS|+
nk2
k*) time. Thus, the total complexity is
O(|
VG|(|
VS|+
nk2
k*)
2(
hS +
k*)).
Theorem 2.2. secondPass returns each optimal reconciliation in O(|VG|(hS + k*)).
P
roof. secondPass starts from the
M(
ρG) and
found in firstPass. It then constructs an optimal solution by visiting each subsequent
g VG, assigning mappings and events by looking up values in
in constant time. Losses are inferred in
O(
k* +
hS) time (see
Vernot et al., 2008). Thus, the complexity for returning each optimal history is
O(|
VG|(
hS +
k*)).
When TS is binary, firstPass is completed in O(hS|VG||VS|2) time, and secondPass reports each optimal solution in O(hS|VG|) time.
Our Notung implementation is efficient in practice. We measured the time required to reconcile 1128 cyanobacterial gene trees with a species tree of size |VS| ≤ 21 for all the parameter settings given in . To assess the effect of polytomy size, we also collapsed edges in the species tree to create a polytomy ranging in size from 2 to 6. The maximum average running time observed on a single AMD Opteron 2.3 ghz, 64-bit processor was ~ 0.05s. per solution.
| Table 1.Event counts for the cyanobacteria dataset, with δ = 3 and λ = 2 |
This article has been ![[not greater-than-or-equal]](/corehtml/pmc/pmcents/x2271.gif)
i u and u 
