A new phylogenetic reconstruction method based on invariants

A new phylogenetic reconstruction method based on invariants M. Casanellas J. Fern´andez-S´anchez Dpt. Matem`atica Aplicada I Universitat Polit`ecnica de Catalunya Av. Diagonal 647. 08028-Barcelona. Spain [email protected], [email protected] Abstract An attempt to use phylogenetic invariants for tree reconstruction was made at the end of the 80s and the beginning of the 90s by several authors (the initial idea due to Lake [Lake, 1987] and Cavender and Felsenstein [Cavender and Felsenstein, 1987]). However, the efficiency of methods based on invariants is still in doubt ([Huelsenbeck, 1995], [Jin and Nei, 1990]), probably because these methods only used few generators of the set of phylogenetic invariants. The method studied in this paper was first introduced in [Casanellas et al., 2005] and it is the first method based on invariants that uses the whole set of generators for DNA data. The simulation studies performed in this paper prove that it is a very competitive and highly efficient phylogenetic reconstruction method, especially for non-homogeneous phylogenetic trees. Introduction Since the introduction of phylogenetic invariants by Cavender and Felsenstein [Cavender and Felsenstein, 1987] and Lake [Lake, 1987], several attempts to give a generating set of polynomial phylogenetic invariants have been made (see for example [Steel et al., 1993], [Ferreti and Sankoff, 1995]) but it has not been until recently that algebraic geometers have managed to find them all ([Allman and Rhodes, 2004a], [Sturmfels and Sullivant, 2005], [Casanellas and Sullivant, 2005]). Methods based on invariants have already proved to be useful in comparative genomics [Sankoff and Blanchette, 1999]). However, they have not had much success in phylogenetic inference ([Jin and Nei, 1990], [Huelsenbeck, 1995]) because only a small set of invariants has been considered (for example, Lake’s method of invariants only used two phylogenetic invariants of degree one among the 795 generators of the 1 set of polynomial invariants of a quartet tree for the Kimura 2-parameter model [Garcia and Porter, ]). But as Felsenstein explains, invariants are worth more attention for what they might lead to in the future ([Felsenstein, 2003]). This future may be the present now since the studies of this paper validate the method based on invariants presented in this report as a promising method. Recently, other methods based on a large set of invariants have also been considered [Eriksson, 2005], [Kim et al., ]. Phylogenetic invariants are equations satisfied by the expected pattern frequencies of a given tree topology Tevolving under an evolutionary model. More precisely, if tis the set of model parameters on Tand pα(t) is the probability of observing the pattern αat the leaves of T, by letting tvary on an open subset of Rd, the probability vector p(t) = (pAA...A(t), . . . , pT T ...T (t)) defines a subset STof dimension ≤dof R4n. A phylogenetic invariant is a real-valued continuous f(x) on R4nsuch that f(p) = 0 for any p∈ST, but not for all the points on the subset ST0deter- mined by another tree topology T0. Essentially, the equations f(p) = 0 are satisfied whatever the parameters of the model are, so they might be used for recovering the tree topology. In practice, the vector of observed pattern frequencies ˆpobtained from an alignment of ntaxa with enough data, should approximate p(t) for some set of parameters ton a tree topology T. In other words, ˆpshould be a point close to the subset STso, if fis an invariant for the topology T, one should have that f(ˆp) is very close to 0. Thus, with probability one, there will be a unique tree topology Tfor which all its phylogenetic invariants are close to 0 when evaluated at ˆp. Therefore, using the phylogenetic invariants for tree reconstruction is a consistent method (see [Hagedorn and Landweber, 2000], [Felsenstein, 2003]). A practical introduction to the theory of invariants can be found in the book of J. Felsenstein [Felsenstein, 2003, chapter 22], whereas the book [Pachter and Sturmfels, 2005] provides a beautiful insight into the applications of algebraic statistics (and in particular polynomial phylogenetic invariants) to computational biology. There are two major motivations for using phylogenetic invariants in tree reconstruction: one of them is the prohibitive computational expense of a full maximum likelihood estimation of a tree and its edge lengths and the other is that the evolutionary models underlying the theory of invariants are non-homogeneous. Indeed, it is known that for biological species different rate matrices should be allowed in different lineages. Thus it is essential to have at our disposal phylogenetic methods for reconstructing trees admitting non-homogeneous models (see [Yang and Yoder, 1999]). In this paper, a phylogenetic reconstruction method that uses polynomial phylogenetic invariants (introduced in [Casanellas et al., 2005]) is studied and tested for quartet unrooted trees evolving under the Kimura 3-parameter model of nucleotide substitution [Kimura, 1981]. Actually, we consider an algebraic Kimura model: the parameters of the model are the entries of the transition matrices (not of the rate matrices). Hence the model is non-homogenous —because it allows different rate matrices along different lineages— but it is stationary and we always assume that all sites are independent and identically distributed (i.i.d. hypotheses). We performed 2 simulation studies to prove its efficiency. Our approach to test the efficiency of the method is taken from Huelsenbeck [Huelsenbeck, 1995] so that a large portion of the tree space is examined to get a general idea of how the algorithm performs. We present the results obtained for sequences of length 100 up to 10000. We also checked the performance of the method on simulated data from nonhomogeneous trees (different rate matrices in different branches) by carrying out a comparison to Neighbor-Joining algorithm [Saitou and Nei, 1987], maximum likelihood algorithm [Felsenstein, 1981], and a non-homogeneous algorithm from PAML [Yang, 1997] for sequences generated under a Kimura 2-parameter model [Kimura, 1980] and different rate matrices along different tree branches. Results Homogeneous data The performance of the invariants method studied here on homogeneous trees can be seen in figures 1 and 2. Using the approach of J.P. Huelsenbeck in [Huelsenbeck, 1995] for quartet trees, we considered two branch-length parameters a, b (see figure 6(a)) and simulated data for each couple of parameters. Parameters aand bwere varied from 0.01 to 0.75 in increments of 0.02 and 1000 alignments were simulated for each couple (a, b) (see figure 6(b)). The simulated trees evolve under the Kimura 3-parameter model [Kimura, 1981] with a fixed rate matrix of the form A C G T A C G T     ·γ α β γ·β α α β ·γ β α γ ·     along the tree (·=−α−β−γ). Figures 1 and 2 show the efficiency of the method considered in this paper for two different rate matrices and for nucleotide sequences of lengths 100, 500, 1000 and 10000. The performance of the method is similar in both cases, although the results are slightly better in the first case (notice that in this case the transition:transversion bias is 1:1 against 2.93:1 in figure 2). This difference in the performance is shown, for example, considering the 95% isocline in the graphic corresponding to sequences of length 100. These figures are to be compared with those shown in Figure A2 of [Huelsenbeck, 1995] (corresponding to phylogenetic inference for sequences generated under a Kimura 2-parameter model of substitution [Kimura, 1980] with 5:1 transition:transversion bias). Though this is of course a biased comparison because our method admits non-homogeneous data and Kimura 3-parameter model, it is worth noticing that our method outperforms many of the methods considered there. In particular, it is clearly better than Lake’s invariant method (which is not surprising because Lake only used linear invariants). At least for sequences of length 3 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Length=100 parameter a parameter b 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Length=500 parameter a parameter b 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Length=1000 parameter a parameter b 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Length=10000 parameter a parameter b Figure 1: The graphics above represent the probability of reconstructing the correct tree in the parameter space (see Figure 6). Black areas correspond to couples (a, b) for which the tree was correctly estimated, while white regions correspond to couples (a, b) for which the tree is never estimated; grey tones indicate areas of intermediate probability. The 95% isocline is drawn in white, while the 33% is drawn in black. The four graphics above show the results obtained under the Kimura 3-parameter model when using a rate matrix with parameters γ= 0.2, α = 0.5, β = 0.3 (hence a 1:1 transition:transversion bias). 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Length=100 parameter a parameter b 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Length=500 parameter a parameter b 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Length=1000 parameter a parameter b 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Length=10000 parameter a parameter b Figure 2: The four graphics above show the results obtained in the four-taxon space when the sequences are generated under the Kimura 3-parameter model of substitution using a rate matrix with parameters γ= 0.12, α = 0.73, β = 0.13 (hence a 2.93:1 transition:transversion bias). 5 500 or larger, our method performs better than neighbor-joining —referred as Minimum Evolution (Kimura1) in figure A2 of Huelsenbeck’s paper. Notice also that the shape of the 95%-isocline for lengths around 500 or larger is quite different from the corresponding shapes in the methods tested by Huelsenbeck: for large values of a(near 0.7), the performance of our method of invariants does not drop drastically (as it happens in most methods considered there). For length 1000, the efficiency of the invariants method considered here is similar to that obtained for lengths ≥10000 in many of the methods tested in [Huelsenbeck, 1995]. From this, it can be inferred that in order to reconstruct the correct tree, much less data is needed in the invariants method presented here than in many other methods (contrary to what was thought until now [Hagedorn and Landweber, 2000]). Non-homogeneous data We tested the invariants method studied in this paper on non-homogeneous data by comparing it with other methods. 1. Comparison with neigbor-joining First of all we compared the performance of the invariants method presented here with neighbor-joining (the algorithm of [Saitou and Nei, 1987]) using Kimura 3-parameter distance [Kimura, 1981]. As it can be seen in figure 3, considering certain non-homogeneous sets of data, we found that the invariants method is more efficient than neighbor-joining. Indeed, we simulated data on an unrooted quartet tree evolving under the Kimura 3-parameter model where different rate matrices at each edge were choosen (see figure 4) and we studied the efficiency of the method when varying the length of nucleotide sequences. In this case, the mean of correctly reconstructed trees for the invariants method is 90.2% whereas the mean for neighbor-joining is 84%. For this particular tree, the maximum likelihood algorithm for the Kimura 3-parameter model reconstructed the tree correctly almost all times, so we do not include the results here. 2. Comparison with maximum likelihood We compared the invariants method with two versions of maximum-likelihood: the usual maximum likelihood for Kimura 2-parameter [Kimura, 1980] and nonhomogeneous maximum likelihood method developped in the package PAML [Yang, 1997] for Kimura 2-parameter model (see details in the methods section). This last method allows different transition/transversion ratio in different tree branches. To perform this comparison we simulated data according to the tree in figure 4 evolving under the Kimura 2-parameter model. The entries in the rate matrices are functions of a parameter εthat is varied from 1 to 9. When ε= 1 the tree is homogeneous and we studied the efficiency of the three methods as ε increases up to 9. 6 0 500 1000 1500 2000 2500 3000 0 20 40 60 80 100 Comparison with Neighbour joining Sequence length number of correctly estimated trees Invariants Neighbor−joining Figure 3: Performance of neighbor-joining and the invariants method on nonhomogeneous data for sequence length varying from 100 to 3000. The non-homogeneous tree used for simulations is described in figure 4. Figure 5 summarizes the results relative to the comparison between our method (invariants), the non-homogeneous method in PAML for Kimura 2-parameter model (PAML non-homogeneous) and the maximum likelihood homogeneous in PAML (PAML homogeneous) for Kimura 2-parameter. It shows the percentage of correctly reconstructed trees for each value of parameter ε. As expected, the homogeneous maximum-likelihood algorithm decreases its efficiency as εincreases. Thus, it is not recommendable to use homogeneous methods on data that can be non-homogeneous. Already for ε= 5, the invariants method overtakes the homogenous maximum likelihood algorithm. As it is deduced from figure 5, our method performs worse than the non-homogeneous algorithm in PAML. However, it is worth noticing that in this test we generated data according to Kimura 2-parameter model and the invariants method presented here was developed under Kimura 3-parameter model. Methods The phylogenetic reconstruction method used in this paper is based on phylogenetic invariants and was first introduced by the first author and L.D. Garcia and S. Sullivant in [Casanellas et al., 2005]. We consider the Kimura 3-parameter model [Kimura, 1981] on an unrooted binary (trivalent) tree of nspecies. The taxa are given by an alignment of nDNA 7 Figure 4: A non-homogeneous tree used for simulations. The numbers labelling each edge correspond to branch lenghts (i.e. expected percent change between the two taxa on the edge). 1. In comparing our method with neighbor-joining, the parameters αi, βi, γiin the Kimura 3-parameter rate matrices Qiwere chosen as: γ1= 1, α1= 4, β1= 1, γ2= 5, α2= 14, β2= 3 ,γ3= 4, α3= 15, β3= 3, γ4= 2, α4= 6, β4= 2, γ5= 2, α5= 3, β5= 1. 2. In comparing our method with maximum likelihood algorithm based on the Kimura 2-parameter model, we chose γ=β= 1 in all rate matrices whereas parameter αwas set as follows: in Q1,α= 4, in Q2,Q4and Q5,α= 3 + ε2, in Q3,α= 3 + ε. When ε= 1 the tree is homogeneous and we studied the efficiency of three methods as εincreases up to 9. sequences of length N. A continous time Markov process along the tree is assumed and we consider that all sites are independently and identically distributed (i.i.d. hypotheses). At each branch ithe substitution matrix has the form A C G T Si= A C G T 0 B B @ aibicidi biaidici cidiaibi dicibiai 1 C C A for some ai, bi, ci, dipositive parameters representing the substitution probabilities along the branch iand satisfying ai+bi+ci+di= 1. Usually a rate matrix Qis fixed and common to the whole tree and Siis the exponential eQti(for some parameter ti representing time). However, in our method we do not make use of rate matrices Q —we only use the substitution matrices— so, as we will see later, the rate matrices might vary along different lineages. We want to stress that the parameters of the model we are considering are the entries of the substitution matrices (we should rather speak of an algebraic Kimura 3-parameter model, according to the book [Pachter and Sturmfels, 2005, chapters 1 and 4]). Phylogenetic invariants Sturmfels and Sullivant [Sturmfels and Sullivant, 2005] gave an explicit description of the generators of the set of polynomial phylogenetic invariants I(T) for an arbitrary tree evolving under a group-based model. For the Kimura 3-parameter model 8 2 4 6 8 0 20 40 60 80 100 Comparison with Maximum Likelihood parameter ε number of correctly estimated trees 123456789 Invariants PAML non homogeneous PAML homogeneous Figure 5: The effect of rate heterogeneity among lineages on method performance. The simulated nucleotide sequences have length 1000 and evolve under the non-homogeneous Kimura 2-parameter model is described in figure 4 (homogeneous for ε= 1). The trees were reconstructed with our method (invariants), with PAML under a non-homogeneous method (PAML non-homogeneous) and the usual PAML maximum likelihood homogeneous (PAML homogeneous). on an unrooted 4-taxa tree the ideal of phylogenetic invariants has 8002 minimal generators (see the webpage [Garcia and Porter, ], [Casanellas et al., 2005] and the discussion in [Sturmfels and Sullivant, 2005], section 7 to see why a smaller subset of invariants does not suffice). According to the results in [Sturmfels and Sullivant, 2005], we produced this generating set for an unrooted tree with 4 leaves under the Kimura 3-parameter model. This requires doing a Fourier transform (or Hadamard conjugation) on the vector of probabilities (pA...A, . . . , pT ...T ) and the phylogenetic invariants are described in terms of the Fourier coordinates. It turns out that the phylogenetic invariants in this case are generated by 144 binomials of degree 2, 1984 binomials of degree 3 and 5874 binomials of degree 4. They can be found in the small trees website [Garcia and Porter, ] (see [Casanellas et al., 2005] for an explanation on the website), and more precisely at http://www.math.tamu.edu/∼lgp/smalltrees/small-trees 30.html. Algorithm Our tree reconstruction algorithm performs the following tasks. Given 4 aligned DNA sequences s1,s2,s3,s4, it first computes the observed relative frequencies of 9 [Pachter and Sturmfels, 2005] Pachter, L. and Sturmfels, B., editors (2005). Algebraic Statistics for computational biology. Cambride University Press. ISBN 0-521-85700-7. [Paradis et al., 2006] Paradis, E., Strimmer, K., Claude, J., Jobb, G., Opgen-Rhein, R., Dutheil, J., Noel, Y., Bolker, B., and Lemon, J. (2006). ape: Analyses of Phylogenetics and Evolution. R package version 1.8-3. [Rambaut and Grassly, 1997] Rambaut, A. and Grassly, N. (1997). Seq-Gen: An application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees. Comput. Appl. Biosci., 13:235–238. [Saitou and Nei, 1987] Saitou, N. and Nei, M. (1987). The neighbor joining method: a new method for reconstructing phylogenetic trees. Mol. Biol. Evol., 4(4):406– 425. [Sankoff and Blanchette, 1999] Sankoff, D. and Blanchette, M. (1999). Phylogenetic invariants for genome rearrangements. J. Comput. Biol., 6:431–445. [Steel et al., 1993] Steel, M., Sz´ekely, L., Erd˝os, P., and Waddell, P. (1993). A complete family of phylogenetic invariants for any number of taxa under kimura’s 3st model. NZ J. Bot., 31:289–296. [Sturmfels and Sullivant, 2005] Sturmfels, B. and Sullivant, S. (2005). Toric ideals of phylogenetic invariants. J. Comput. Biol., 12:204–228. [Team, 2005] Team, R. D. C. (2005). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3- 900051-07-0. [Yang, 1997] Yang, Z. (1997). PAML: A program package for phylogenetic analysis by maximum likelihood. CABIOS, 15:555–556. [Yang and Roberts, 1995] Yang, Z. and Roberts, D. (1995). On the use of nucleic acid sequences to infer early branchings in the tree of life. Mol. Biol. Evol., 12:451–458. [Yang and Yoder, 1999] Yang, Z. and Yoder, A. D. (1999). Estimation of the transition/transversion rate bias and species sampling. J. Mor. Evol., 48:274–283. 16