PhyloNet: a software package for analyzing and reconstructing reticulate evolutionary relationships

Tree-based comparison

Let $\mathcal{T}$(N₁) and $\mathcal{T}$(N₂) be the two sets of all trees induced by the two networks, and let d(·,·) be a distance metric on trees (see [26] for examples of such metrics). The idea is to compare the two networks based on how similar their corresponding sets of trees are. We formalize this as follows. Construct a weighted complete bipartite graph G(U₁, U₂, E), where |U_i| = |$\mathcal{T}$(N_i)|, and there are two bijections f_i : U_i → $\mathcal{T}$(N_i) for i = 1, 2. The weight of an edge e = (u, v) ∈ E for u ∈ U₁ and v ∈ U₂, w(e) = d(f₁(u), f₂(v)). Then, the tree-based measure m^tree(N₁, N₂) is defined as the weight of a minimum-weight edge cover of G divided by, the number of the edges in the cover. In its current implementation, PhyloNet uses the Robinson-Foulds distance measure [27] for d. The tree-based measure was first introduced by Nakhleh et al. [28]. An illustration of tree-based comparison of the two networks N₁ and N₂ in Figure 2 is given in Figure 6. Shown on the left of Figure 6 is the bipartite graph G built from the sets $\mathcal{T}$(N₁) and $\mathcal{T}$(N₂) of trees induced by the two networks; these two sets are shown in Figure 5. The weight of each edge connecting two nodes in G is the RF distance between the two trees corresponding to these two nodes. These weights can be normalized by the number of internal edges in the trees. Since each of the eight trees has six internal edges, the weight of each edge in G can be divided by six to normalize it.

Shown on the right of Figure 6 is the minimum-weight edge cover of G, which is the set of edges that satisfies two conditions: (1) each node in G must be the endpoint of at least one edge in the set, and (2) the sum of the weights of the edges in the set is minimum among all sets of edges satisfying condition 1. In this case, the four edges shown are a cover, since each node in G is "covered" by at least one edge (here, each node is covered by exactly one edge). Further, it is of minimum weight, which equals 2, since a simple inspection yields that every other cover has a weight larger than 2. Since the cover has four edges in it, we have m^tree(N₁, N₂) = (0 + 0 + 1/6 + 1/6)/4 = 1/12. If we use the raw RF values, then m^tree(N₁, N₂) = (0 + 0 + 1 + 1)/4 = 1/2.

Cluster-based comparison

Let C₁ = $\mathcal{C}$(N₁) and C₂ = $\mathcal{C}$(N₂) be the two sets of all clusters induced by the two networks. We define the measure based on these two sets to be

The rationale behind this measure is that it is the sum of the ratios of clusters present in one but not both networks. The cluster-based measure was first introduced by Nakhleh et al. [29]. The sets C₁ = $\mathcal{C}$(N₁) and C₂ = $\mathcal{C}$(N₂) of the two networks N₁ and N₂ in Figure 2 are listed in Table 1, with |C₁| = |C₂| = 14. Since |C₁ - C₂| = |C₂ - C₁| = 2 (the two highlighted clusters in Table 1), we have m^cluster(N₁, N₂) = 1/7. A similar weighting scheme to that described in the previous section can be used to incorporate the fraction of trees in which a cluster appears into the measure calculation.

Tripartition-based comparison

Let θ₁ = θ(N₁) and θ₂ = θ(N₂) be the two sets of all tripartitions induced by the two networks. We define the measure based on these two sets to be

This measure views the two networks in terms of the sets of edges they define (where an edge is in a 1-1 correspondence with a tripartition) and computes the sum of the ratios of edges present in one but not both networks. The tripartition-based measure was devised by Moret et al. [19]. The sets θ₁ = θ(N₁) and θ₂ = θ(N₂) of the two networks N₁ and N₂ in Figure 2 are listed in Table 2, with |θ₁| = |θ₂| = 12. Since |θ₁ - θ₂| = |θ₂ - θ₁| = 1 (the highlighted tripartition in Table 2), we have m^tripartition(N₁, N₂) = 1/12.

Which measure to use?

Several distance measures, such as the Robinson-Fould measure and the Subtree Prune and Regraft (SPR) distance, have been introduced over the years to quantify the difference between the topologies of a pair of phylogenetic trees; e.g., see [25,26] for description of many of these measures. Even though these measures may compute different distance values on the same pair of trees, there has been no consensus as to which measure should be used in general [30]. It may be the case that the Robinson-Foulds measure is more commonly used than the others, but this may be a mere reflection of its very low time requirements as compared to the other, more compute-intensive, measures.

Regarding the three measures for comparing networks, a scenario analogous to that in phylogenetic trees arises here: each measure gives a different quantification of the dissimilarity between two networks based on one of the three ways to characterize a given network. As shown in the examples above, some or all of these measures may compute the same value for a given pair of networks, but that may not always be the case. Tree-based comparison of networks can be viewed as a method to quantify how similar, or dissimilar, two networks are in terms of their quality as a summary of a collection of trees. In some cases, even though two networks "look different," they may be identical in terms of the trees they induce – this is the issue of indistinguishability of networks from a collection of trees that Nakhleh and colleagues discussed in [19]. In such a case, using the tree-based comparison, or equivalently the cluster-based comparison, is most appropriate. However, if the similarity/dissimilarity of two networks means something close to an isomorphism, then the tripartition-based measure is more appropriate. However, it is important to note that none of the three measures described here is a metric on the general space of all evolutionary networks labeled by a given set of taxa.

A practical distinction among the three measures can be derived based on the methods used to infer the evolutionary history of the set of species under study. Methods such as SplitsTree [23] and NeighborNet [24] represent the evolutionary history as a set of splits, or clusters, hence making it more natural to use cluster-based comparison to study their performance. Methods such as RIATA-HGT [16] and LatTrans [14] compute evolutionary networks that are rooted, directed, acylic graphs, where internal nodes have an evolutionary implication in terms of ancestry. For these two methods, all three measures are appropriate. When the evolutionary history of a set of species is represented as a collection of its constituent gene trees, the tree-based measure is most appropriate.

Finally, a clear distinction can be made among the methods in terms of computational requirements. In their current implementations, the tripartition-based measure is very fast in practice, taking time that is polynomial in the size of the two networks. On the other hand, the tree- and cluster-based measures are much slower, taking time that is exponential in the number of network nodes in the two networks (since these measures compute explicitly all trees inside each of the two networks). In light of recent complexity results that we obtained [31], it is very likely that no polynomial-time algorithms exist for computing the tree- and cluster-based measures in general.

Parsimony of evolutionary networks

Nakhleh and colleagues have recently formalized a maximum parsimony (MP) criterion for evolutionary networks [32] and demonstrated its utility in reconstructing evolutionary networks on both biological and synthetic data sets [33]. In this section, we describe a PhyloNet utility that allows for comparing two evolutionary networks in terms of their fitness to the evolution of a sequence data set, based on the MP criterion. We first begin with a brief review of the MP criterion, based on the exposition in [32].

The relationship between an evolutionary network and its constituent trees, as described in the background section, is the basis for the MP extension to evolutionary networks.

Definition 1 The Hamming distance between two equal-length sequences x and y, denoted by H(x, y), is the number of positions j such that x_j ≠ y_j.

Given a fully-labeled tree T, i.e., a tree in which each node v is labeled by a sequence s_v over some alphabet Σ, we define the Hamming distance of an edge e ∈ E(T), denoted by H(e), to be H(s_u, s_v), where u and v are the two endpoints of e. We now define the parsimony score of a tree T.

Definition 2 The parsimony score of a fully-labeled tree T, is Σ_{e ∈ E(T)} H(e). Given a set S of sequences, a maximum parsimony tree for S is a tree leaf-labeled by S and assigned labels for the internal nodes, of minimum parsimony score.

The parsimony definitions can be extended in a straightforward manner to incorporate different site substitution matrices, where different substitutions do not necessarily contribute equally to the parsimony score, by simply modifying the formula H(x, y) to reflect the weights. Let Σ be the set of states that the two sequences x and y can take (e.g., Σ = {A, C, T, G} for DNA sequences), and W the site substitution matrix such that W[σ₁, σ₂] is the weight of replacing σ₁ by σ₂, for every σ₁, σ₂ ∈ Σ. In particular, the identity site substitution matrix satisfies W[σ₁, σ₂] = 0 when σ₁ = σ₂, and W[σ₁, σ₂] = 1 otherwise. The weighted Hamming distance between two sequence is H(x, y) = Σ_{1 ≤ i ≤ k} W(x_i, y_i), where k is the length of the sequences x and y. The rest of the definitions are identical to the simple Hamming distance case. As described above, the evolutionary history of a single (non-recombining) gene is modeled by one of the trees contained inside the evolutionary network of the species containing that gene. Therefore the evolutionary history of a site s is also modeled by a tree contained inside the evolutionary network. A natural way to extend the tree-based parsimony score to fit a dataset that evolved on a network is to define the parsimony score for each site as the minimum parsimony score of that site over all trees contained inside the network.

Definition 3 ([32]) The parsimony score of a network N leaf-labeled by a set S of taxa, is

where TCost(T, s_i) is the parsimony score of site s_i on tree T.

Notice that as usually large segments of DNA, rather than single sites, evolve together, Definition 3 can be extended easily to reflect this fact, by partitioning the sequences S into non-overlapping blocks b_i of sites, rather than sites s_i, and replacing s_i by b_i in Definition 3. This extension may be very significant if, for example, the evolutionary history of a gene includes some recombination events, and hence that evolutionary history is not a single tree. In this case, the recombination breakpoint can be detected by experimenting with different block sizes.

The MP utility in PhyloNet allows the user to specify two evolutionary networks (either or both of which can be a tree) N₁ and N₂ and a sequence data set S, and computes the parsimony scores NCost(N₁, S) and NCost(N₂, S). The user can then compare the two scores and evaluate the fitness of the networks to the data set S based on the difference in the scores. Further, the utility allows the user, for example, to evaluate the significance of each network edge in a network N by comparing the parsimony scores of two different versions of N that contain different subsets of the network edges in N.

Computing multiple solutions and the graphical output

RIATA-HGT was designed originally to compute a single solution to the problem, and was mainly aimed at binary trees. Later, Than and Nakhleh [17] extended the method to compute multiple solutions and to handle non-binary trees. These two features are very significant: the former allows biologists to explore multiple potential HGT scenarios, whereas the latter allows for analyzing trees in which some edges were contracted due to inaccuracies (see [36] for example). We have conducted an experimental study to compare the performance of RIATA-HGT with LatTrans [18]. Although RIATA-HGT and LatTrans [14] have almost the same performance in terms of the number of HGT solutions and the solution size, the former runs much faster than the latter.

For a compact representation of multiple solutions, we introduce four terms:

• An event: this is a single HGT edge, written in the form of X → Y, where X and Y are two nodes in the species tree.

• A subsolution: this is an atomic set of events, which forms a part of an overall solution. In other words, either all or none of the events of a subsolution are taken in a solution.

• A component: a set of components and/or subsolutions. Two components at the same level of decomposition are independent, in that an element of each component is needed to form a solution.

• A solution: the union of a single element from each component at the highest level.

To illustrate these concepts, consider species tree (((a, b), c), (d, (e, f))) and the gene tree (((a, c), b), ((d, f), e)). Observe, that each HGT event required to reconcile the two trees has both endpoints in the subtree ((a, b), c) or both endpoints in the subtree (d, (e, f)), and no HGT event has endpoints in both subtrees. In this case, RIATA-HGT divides the pair of trees into two pairs:

• Pair 1: ((a, b), c) and ((a, c), b)

• Pair 2: (d, (e, f)) and ((d, f), e),

and solves the HGT Reconstruction Problem on each of the two pairs independently. The set of solutions of each pair is a component. Notice that for each pair there are three possible ways to reconcile them; each such way is a called a subsolution. Each subsolution is a set of events, which in this case is only one event. Figure 7(a) shows the screen captures of two graphical outputs that correspond to two solutions on this pair of trees. Notice that if a component can be further divided into independent components, RIATA-HGT would do so, which will result in components at different levels, with the largest components being at the highest level.

The compact representation of RIATA-HGT's output in terms of subsolutions and components is especially helpful when the number of solutions is large. RIATA-HGT also has an option to display all complete solutions. RIATA-HGT enumerates all complete solutions that are compactly represented as described in the preceding paragraphs. Each solution, which is a set of HGT events, along with the species tree defines an evolutionary network, which RIATA-HGT displays in the eNewick format. For example, for the trees (((a, b), c), (d, (e, f))) and (((a, c), b), ((d, f), e)), RIATA-HGT outputs 9 different networks in the eNewick format, if RIATA-HGT's option for displaying complete solutions is on. Figure 7(b) shows the corresponding eNewick representations.

From the multiple comparisons between a species and a set of trees, RIATA-HGT offers a (strict) consensus network. For each pair of species tree and gene tree, RIATA-HGT computes a set of HGT events for reconciling them. To obtain the consensus network, RIATA-HGT retains only HGT events that appear in every set of solutions for every pair of species tree and gene tree. Those events are then added to the species to build the consensus network.

We note here that while offering a simple summary of solutions, this way of computing consensus networks may not be appropriate in general; work is under way to address this issue more properly.

Finally, RIATA-HGT may report '[time violation?]' next to an inferred HGT X → Y. If this is the case, this indicates that node X lies on the path from Y to the root of the species tree. Theoretically, this indicates that two nodes that do not co-exist in time, X and Y in this case, shared genetic material, and hence the warning of 'time violation.' However, this may be the case simply due to incomplete taxon sampling, as discussed in [19]. Therefore, the warning is issued in this case so as to alert that user that this inferred HGT edge is worth further inspection.

Assessing the support of HGT edges

In [37] we proposed a method for assessing the support of HGT edges. Roughly speaking, the support value of HGT edge X → Y in the species tree, where Y'is the sibling of Y, is derived from the bootstrap values of the gene tree branches that separate the clade under Y from the clade under Y'. The rationale behind the idea is that if Y' and Y are well separated in the gene trees (i.e., some branches in the path from Y to Y' have high bootstrap values), HGT is necessary to move Y away from Y'). For example, the support of HGT edge X → Y in Figure 8(a) is calculated based on the bootstrap values of the branches separating B from A in the gene tree, and it is 80 (which is the maximum bootstrap value of all edges on the path separating A and B in the gene tree). More technical details can be found in [37].

Than et al. [37] have studied the reliability of this method for assessing the support of HGT edges on various data sets from [38]. In this paper, we illustrate the output of RIATA-HGT on a pair of species/gene trees from [38], as shown in Figure 9. The output of RIATA-HGT on this pair of trees is shown in Figure 10. RIATA-HGT computed four solutions, each of which has nine HGT edges. To allow for a compact representation of the solutions, they are divided into two components (which are computed automatically by RIATA-HGT), and each solution is formed by taking one subsolution from each component. HGT edges for the solutions are divided into two components, which means that a complete solution is formed by taking one solution from each component. Each component is labeled by the name of the internal node that is the root of the clade corresponding to that component. In the case of the solutions presented in Figure 10, each solution contains nine HGT edges, eight of which form a single subsolution in Component I18 and the ninth is the only edge in the only subsolution in Component I8. The value in parentheses next to each HGT edge is its support value computed from the bootstrap values of the gene tree branches (Figure 9(b)). Bergthorsson et. al. [38] reported three HGTs involving Amborella: one HGT with donor being a species in the Moss group (species Brachythecium, Hypnum, and Thuidium, under the internal node I14 in the species tree) and the other two with donors being species in the Eudicot group (species Arabidopsis, Beta, Brassica, Daucus, Petunia, and Oenothera, under the internal node I5 in the species tree). The HGT from Moss has high SH support value [39]. RIATA-HGT finds this event, I14 → Amborella_H_M, with bootstrap value 98.0%. The other two HGT events from Eudicot do not have significant SH support values. RIATA-HGT also finds these events, I5 → Amborella_H_E1 and I2 → Amborella_H_E2. However, their support values are 73.0%, much smaller than that of the event from Moss. In addition to these three HGT edges, RIATA-HGT identified six more edges, four out of which had support values smaller than 50.0% (RIATA-HGT does not display support values that are smaller than 50.0%). The HGT edge from Component I8, which is shared among all four solutions, has support value of 100.0%. This edge was not reported in [38]. A similar situation arises with the HGT edge I5 → I8, which is part of the three solutions that contain subsolutions 1, 3, and 4 from Component I18: the HGT edge has support of 100.0%, and was not reported in [38], which may be a reflection that the authors focused only on HGT events involving Amborella. The ninth HGT edge in Subsolution2 of Component I18 has support value smaller than 50.0%.