Stochastic errors vs. modeling errors in distance based phylogenetic reconstructions

  • Daniel Doerr1,

    Affiliated with

    • Ilan Gronau2,

      Affiliated with

      • Shlomo Moran3Email author and

        Affiliated with

        • Irad Yavneh3

          Affiliated with

          Algorithms for Molecular Biology20127:22

          DOI: 10.1186/1748-7188-7-22

          Received: 23 December 2011

          Accepted: 28 June 2012

          Published: 31 August 2012



          Distance-based phylogenetic reconstruction methods use evolutionary distances between species in order to reconstruct the phylogenetic tree spanning them. There are many different methods for estimating distances from sequence data. These methods assume different substitution models and have different statistical properties. Since the true substitution model is typically unknown, it is important to consider the effect of model misspecification on the performance of a distance estimation method.


          This paper continues the line of research which attempts to adjust to each given set of input sequences a distance function which maximizes the expected topological accuracy of the reconstructed tree. We focus here on the effect of systematic error caused by assuming an inadequate model, but consider also the stochastic error caused by using short sequences. We introduce a theoretical framework for analyzing both sources of error based on the notion of deviation from additivity, which quantifies the contribution of model misspecification to the estimation error. We demonstrate this framework by studying the behavior of the Jukes-Cantor distance function when applied to data generated according to Kimura’s two-parameter model with a transition-transversion bias. We provide both a theoretical derivation for this case, and a detailed simulation study on quartet trees.


          We demonstrate both analytically and experimentally that by deliberately assuming an oversimplified evolutionary model, it is possible to increase the topological accuracy of reconstruction. Our theoretical framework provides new insights into the mechanisms that enables statistically inconsistent reconstruction methods to outperform consistent methods.


          Phylogenetic reconstructions Substitution models Additive substitution rate functions


          Phylogenetic reconstruction is the task of determining the topology of an evolutionary tree underlying a given set of samples (species) using sequence data extracted from them. This is typically done by assuming some simplified model for DNA sequence evolution, in most cases modeling it as a homogeneous continuous-time Markov process [13]. Distance-based reconstruction algorithms tackle this task by first computing a set of pairwise distances between the n input samples and then finding a tree which fits these distances. The distance measures used for this purpose typically reflect the rates of certain substitution events along the evolutionary paths in question. We thus refer to these distance measures as substitution rate (SR) functions. The distance-based approach is based on the fact that if the SR function used is additive for the underlying substitution model, and the input sequences are sufficiently long, then the topology of the true tree can be efficiently recovered with high probability. However, since the underlying evolutionary model is usually unknown, this assumption is rarely satisfied in practice.

          Substitution models used for phylogenetic reconstruction range from the simplest Jukes-Cantor (JC) model [4], through slightly more complex and flexible models, such as Kimura’s two-parameter (K2P) model [5] and the Hasegawa-Kishino-Yano model (HKY) [6], to the General Time-Reversible (GTR) model [7, 8]. In previous works [9, 10] we observed that substitution models which are not too restrictive or too general have many inherently different additive SR functions. We used this basic observation to demonstrate that it is possible to adjust for each given set of DNA sequences a “good” additive SR function, which leads to significantly increased phylogenetic reconstruction accuracy, compared to other additive SR functions. This exploits our ability to predict the stochastic noise associated with each SR function. When the SR function used for distance estimation is additive for the underlying substitution model, this stochastic noise is the only cause for inaccurate reconstruction. However, in the scenario, which is very common in practice, where the SR function in use is not additive for the model, an additional systematic bias is introduced in the distance estimates. This systematic bias in distance estimation results in a phylogenetic reconstruction method that might be statistically inconsistent in some cases. In this papera, we extend our previous line of research to this scenario, by removing the constraint of additivity. We do this by considering both the stochastic noise and systematic error.

          Several previous studies have demonstrated the utility of phylogenetic reconstruction methods that are not generally statistically consistent. The maximum parsimony method has been long known to be inconsistent in some cases [11, 12]. However, in other cases it was shown to be more likely to produce accurate reconstructions, compared with the maximum likelihood method [1315]. More recently, it has been demonstrated that reconstruction accuracy can be improved by deliberately assuming an oversimplified substitution model, when reconstructing a tree using maximum likelihood [16, 17]. In the context of distance-based reconstruction, non-additive distance measures have been shown in several cases to lead to improved accuracy when compared with additive measures [18, 19]. Overall, these studies provide convincing evidence for the need to consider inconsistent phylogenetic reconstruction methods. However, none of them provide a rigorous framework for characterizing the cases in which inconsistent methods outperform consistent ones.

          In this paper we develop a theoretical framework which provides a practical and systematic way to quantify the effect of distance-estimation-bias on the accuracy of distance-based reconstruction. This framework is based on a novel method for measuring the deviation from additivity of SR functions. Coupled with the results in [9], this method enables evaluation of both the systematic bias and stochastic noise of SR functions. Such evaluation is important, because there is often a tradeoff between these two sources of error, stemming from the fact that simpler models with fewer parameters (such as JC) have smaller stochastic noise at the expense of greater estimation bias. Our framework allows us to consider this tradeoff when deciding which SR function to use for a given data set. This allows us to characterize a wide range of cases in which an SR function associated with an oversimplified evolutionary model results in increased reconstruction accuracy.

          This finding falls in line with previous studies demonstrating the usefulness of phylogenetic reconstruction methods that are not generally consistent. Previous studies have attributed the increased accuracy of inconsistent methods mainly to the fact that these methods have a bias toward reconstructing certain topologies, leading to increased accuracy in cases where the phylogeny being reconstructed has the “favored topology”. We notice a similar behavior using our theoretical characterization of non-additive SR functions. However, somewhat surprisingly, we find that non-additive SR functions often have an advantage even when the phylogeny being reconstructed has an “unfavorable topology”. This is due to the reduced stochastic noise of the non-additive SR function (compared with it additive alternatives), which compensates for its topological bias.

          Our paper is organized as follows. Section “Background” outlines some of the required background and introduces several new concepts that are central in our analysis. Section “Deviation from additivity in homogeneous substitution models” provides the main analytic results in the paper, and introduces deviation from additivity as a measure of distance estimation bias. In that section we prove a general upper bound for this deviation and establish a connection with reconstruction accuracy. We then study deviation from additivity and stochastic error of the JC distance formula when applied to data generated under the K2P model. In Section “Performance of Non affine-additive SR functions in quartet resolution” we study the effect of deviation from additivity and stochastic error on the accuracy of quartet reconstruction. In the case of quartets we can draw a tight connection between the different sources of error in distance estimation and inaccuracy of reconstruction. We present a useful heuristic, based on the so-called Fisher criterion ( [20, 21]), for comparing the expected accuracy of two SR functions in this context. In Section “Simulations on Hasegawa’s Tree” we extend our study to larger trees using experiments on simulated data based on the tree obtained by Hasegawa in [6]. Finally, In Section “Inferring trees from genomic sequences” we demonstrate our approach through a series of experiments reconstructing trees from bacterial gene sequences.


          In this section we provide a brief exposition of DNA substitution models and substitution rate functions used for distance estimation. We concentrate on details essential to this study and refer the reader to a previous paper [9] and standard textbooks [1, 2] for a more complete survey.

          Substitution Models

          In this work, a DNA substitution model is simply a set of stochastic 4×4 transition matrices closed under matrix product (i.e., P Q P Q ). These matrices serve to describe the substitution process along evolutionary paths in a phylogenetic tree. All substitution models addressed in this paper are time-reversible [7]. A model tree in a time reversible substitution model , or an tree, is an undirected tree t =(V E ) in which each edge eE is associated with a transition matrix . An -tree t implies an inter-leaf transition matrix for each pair of leaves , namely . Most common models are defined using rate matrices, which are 4×4 matrices whose off-diagonal elements are non-negative substitution rates, and whose rows sum to 0. A stochastic transition matrix P is obtained from a rate matrix R through matrix exponentiation: P= e R .

          A common assumption made on the substitution process is that it is homogeneous throughout time. This means that all rate matrices in the model are proportional to each other. Such a substitution model is thus termed homogeneous, and it is defined by a unit rate matrix R as follows: . Note that the definition of the unit rate matrix associated with a given homogeneous model is somewhat arbitrary b, but once the unit R is defined, it implies a bijection (or equivalence) between rate matrices in and the parameter t, which corresponds to evolutionary time. We will make use of this equivalence extensively throughout this paper.

          We use the Kimura’s two-parameter (K2P) model [5] as a concrete example for demonstrating our approach. A rate matrix in this model is defined by two rate parameters: α, which is the rate of transition -type (ti) substitutions ( , ), and β, which is the rate of transversion -type (tv) substitutions ( ). Each K2P rate matrix can be represented as a product of a unit rate matrix, in which α + 2β =1, and a scalar t corresponding to evolutionary time .


          Each unit rate matrix of the K2P model defines a homogeneous sub-model, which is identified by its unique transition-transversion (ti-tv) ratio . The Jukes-Cantor (JC) model [4] is a special homogeneous sub-model of K2P, in which (i.e., α =β ). Although the K2P model is defined in (1) as a union of its homogeneous sub-models, it is important to note that this union is closed under matrix product, implying that K2P adheres to our definition of a proper substitution model. Conversely, some commonly used substitution models, such as GTR and HKY, are defined as a union of homogeneous models, but are not themselves closed under matrix product [22].

          Transition matrices in the K2P model have the same symmetric structure as the underlying rate matrices, with two distinct transition parameters: p α – the probability of a transition-type substitution; p β – the probability of a transversion-type substitution. The transformations between (α,β,t ) and (p α ,p β ) are given by the following equations:



          Substitution rate functions

          A substitution rate (SR) function for a model is a non-negative continuous function that maps each transition matrix onto a numerical value of “substitution rate”. An SR function Δ induces the following dissimilarity mapping over the leaves of an -tree t : , for all . Of particular interest in phylogenetic reconstruction are additive SR functions.

          Definition 2.1 (Additive SR function)

          An SR function Δ is said to be additive for a substitution model if for all , Δ (P Q)=Δ (P) + Δ (Q).

          It is often explicitly required that an SR function be additive for the assumed model (see [9]). The evolutionary time, t, typically serves as the standard additive measure in most common substitution models. Throughout this study we follow the special case of K2P, focusing on the two SR functions defined below.



          The first SR function, Δ K2P, is the common SR function suggested for the K2P model in [5], and it is clearly additive, as it maps the transition probabilities onto evolutionary time t . The second SR function, Δ JC, maps the transition probabilities onto evolutionary time only in the special case of the JC model where α =β . Under other homogeneous sub-models of K2P, it is non-additive. This non-additivity is analyzed in details in section Deviation from additivity in homogeneous substitution models.

          Additive metrics, Affine-additive mappings, and Near-additivity

          The core idea behind distance-based phylogenetic reconstruction is that a phylogenetic tree t can be accurately and efficiently reconstructed from pairwise distances which are additive with respect to t[23, 24].

          Definition 2.2 (Additive metric)

          A metric d defined over the leaf-set L (t ) of a tree t is t -additive (or additive w.r.t t ), if there exists a positive edge-weighting function , such that for each i,jL (t ), . d is additive for a set S if it is t -additive for some tree t where L (t )=S .

          It is well known that additive SR functions imply additive metrics: if Δ is an additive SR function for a model , then for any -tree t, (the dissimilarity mapping induced by Δ on t ) is a t -additive metric. The inherent difficulty in reconstructing phylogenies using additive SR functions is that computing the implied t -additive metric requires the exact values of the inter-taxon transition matrices {P ij }, and getting these exact values from alignments of finite length is practically impossible. Therefore, a distance-based reconstruction algorithm is useful in a realistic setting only if it has some robustness to error in distance estimation. In [25], Atteson observed that the topology of a phylogenetic tree t can be accurately (and efficiently) reconstructed from any dissimilarity mapping d which is sufficiently close to a t -additive metric, using certain “robust” distance-based algorithmsc. Formally, “sufficiently close” is defined by the following relation:

          Definition 2.3 (Near-additive mapping)

          A dissimilarity mapping d on L (t ) is said to be near-additive w.r.t. t iff there exists a t -additive mapping d s.t.


          where is the minimal weight assigned to an internal edged by the edge weighting function corresponding to the additive metric d.

          For our results we will be using a generalization of this criterion, in which the mapping d can be any affine-additive mapping, defined below.

          Definition 2.4 (Affine-additive mapping)

          A dissimilarity mapping is said to be affine-additive w.r.t. a phylogenetic tree t, if there is a t -additive metric d, and scalars a >0,B s.t. (i.e., for all ).

          As with additive metrics, affine-additive mapping are also associated with edge weights. Let d be a t -additive mapping corresponding to the edge-weighting function w (·). Then the edge weighting function corresponding to the affine additive mapping is given by: for all internal edges, and for all external edges. When B is positive, is actually an additive metric, but when B is negative, the weights of external edges implied by might be negative, and might even yield negative dissimilarities. The generalization of Atteson’s theorem to cases where d is affine-additive follows from the observation that the robust distance-based reconstruct algorithms considered by Atteson are invariant to affine transformations of their input distances. From this point on, when we say a dissimilarity mapping d is near additive, we mean it satisfies (6) with respect to some affine-additive mapping d .

          Local consistency

          Atteson’s result plays a central role in arguing the statistical consistency of distance-based phylogenetic reconstruction. Typically, this is done by assuming that the inter-leaf distances are computed using an SR function Δ which is additive for the underlying substitution model , as follows:
          1. 1.

            If Δ is additive for , then for each -tree t the mapping defined by for all i,jL (t ), is a t -additive metric.

          2. 2.

            As the length of the input sequences grows, the estimated transition matrices converge (w.h.p.) to the true matrices {P ij }.

          3. 3.

            When are sufficiently close to {P ij }, the estimated dissimilarity map defined by is sufficiently close to , and is thus near-additive.

          4. 4.

            The near-additivity of the estimated dissimilarity map implies accurate topological reconstruction, assuming a robust distance-based algorithm is used.


          This line of argument has been used in numerous works studying statistical consistency of distance-based algorithms (e.g., [2527]), and in all these cases an additive SR function is assumed. Notice, however, that this line of argument remains valid when is near additive w.r.t. t . For instance, consistent reconstruction of any -tree is guaranteed by using an affine-additive SR function , which is an affine transformation of some additive SR function Δ : (with a >0). An SR function that is not affine-additive in a given substitution model does not guarantee consistency across all -trees, but it still can be consistent for specific -trees.

          Definition 2.5 (Consistent SR function)

          An SR function Δ of a substitution model is said to be consistent w.r.t. an -tree t if is near-additive w.r.t t .

          The main idea endorsed in this paper is that if an SR function only deviates slightly from some SR function which is affine-additive for , then it might be consistent with respect to many -trees of interest, and as such should be considered for use in distance based reconstructions.

          Deviation from additivity in homogeneous substitution models

          In order to assess whether a given SR function Δ is consistent w.r.t. a given model tree t, one has to find an affine-additive mapping d which minimizes the ratio (see Definition 2.3). This task seems hard in a general setting, but in the special case of homogeneous substitution models it is tractable. Consider a homogeneous substitution model . The unit rate matrix Rimplies a 1-1 mapping between evolutionary time t and rate matrices in . It is thus useful to view an SR function for as a function which maps the evolutionary time t to a dissimilarity measure Δ (t ).

          It can be shown that such Δ is affine-additive in the model if and only if Δ (t )=at + B for some . We define the deviation of an SR function Δ from a given affine-additive function at + B in an interval [t 0,t 1] as (the factor normalizes the deviation to units of evolutionary time). The deviation from additivity of Δ within [t 0,t 1] is defined as the minimum deviation of Δ from any affine-additive function in that interval.

          Definition 2.6 (Deviation from additivity)

          Let be an SR function in a homogeneous substitution model. The deviation from additivity of Δ in an interval [t 0,t 1] is defined by:


          Lemma 2.7 below presents the basic relation between deviation from additivity and consistency. In Section Performance of Non affine-additive SR functions in quartet resolution we demonstrate the tightness of this relation.

          Lemma 2.7

          Let be a homogeneous model, and let t be an -tree with edge lengths (measured in time units) denoted by {t e }. Let tmin=min{t e :et }, and assume that all inter-leaf distances in t fall within the interval [t 0,t 1]. Then any SR function Δ in for which is consistent w.r.t. t .


          We need to show that is near-additive w.r.t. t . Since , there are which satisfy

          For all i,jL (t ), denote , and let d be the dissimilarity map associated with evolutionary time: d (i,j )=t ij . Clearly, d is an additive metric, and the dissimilarity mapping is an affine-additive mapping. The internal-edge-weights associated with are given by (see discussion following Definition 2.4), implying that . We thus have:

          An upper bound on the deviation of an SR function Δ from additivity in a given interval [t 0,t 1] is implied from the error associated with its linear interpolation At + B within that interval ( and ). Figure 1a demonstrates this for Δ JC under a homogeneous sub-model of K2P, and Lemma 2.8 below presents a general upper bound on the deviation from additivity. For this purpose, we assume that the SR function Δ is a monotone increasing continuous function of t with continuous first and second derivatives.

          Figure 1

          Deviation from additivity and stochastic error. (a) Δ JC is portrayed (green) in the homogeneous sub-model of K2P with R =10 in the interval t ∈[0.8,2]. Its linear interpolation in that interval, Δ int=At + B, is plotted in blue, and the maximum difference between the two functions is designated by X. The deviation of Δ JC from additivity in this setting is . (b) The affine-additive SR function minimizing its deviation from Δ JC is . The two SR functions Δ JC and Δ are shown with their stochastic error margins, assuming sequence length of 500 bp.

          Lemma 2.8

          Let be an SR function in a homogeneous substitution model, and let [t 0,t 1] be an interval. Let Δ int(t )=At + B be the linear interpolation of Δ in [t 0,t 1] defined above, and let . Then



          Let us start by introducing a couple of auxiliary notations:

          We are looking for and which minimize . Let ψmin=mint ∈[t 0,t 1]{ψ (A,B,t )}, ψ max=maxt ∈[t 0,t 1]{ψ (A,B,t )}, and let . Then . A bound for dev (Δ,[t 0,t 1]) will thus follow by showing that .

          Since Δ int(t )=At + B is a linear interpolation of Δ in t 0 t 1, we have ψ (A B t 0)=ψ (A B t 1)=0. Let t min be an arbitrary point in the interval t 0 t 1s.t. ψ (A B t min)=ψ min≤0 and let (t 2 t 3) be the maximal open interval in t 0 t 1 containing t min in which ψ (A B t )<0 (this interval can be empty if ψ min=0). We define a similar interval (t 4 t 5) in which ψ (A B t )>0 around some arbitrary t maxs.t. ψ (A B t max )=ψ max. Note that the intervals (t 2 t 3) and (t 4 t 5) are disjoint, and that Δ int is the linear interpolation of Δ in both these intervals (since ψ (A B t 2)=ψ (A B t 3)=ψ (A B t 4)=ψ (A B t 5)=0). Therefore, the bound on the error of polynomial interpolation (see, e.g., [28], p. 187) implies that

          Combining these, we get


          Deviation of Δ JCfrom Additivity in K2P

          We now turn to study the deviation of Δ JC from additivity in homogeneous sub-models of K2P with ti-tv ratio . First, we express Δ JCas a function of the ti-tv ratio R and the time t, using (5) and the relations and α + 2β =1.

          Note that the homogeneous K2P sub-model with is the JC model; in this case the second term of (10) vanishes, leaving . For other homogeneous sub-models of K2P, where , Δ JC is not affine-additive (i.e., not of the form at + B for a >0), and we can use the result in Lemma 2.8 to bound the deviation of Δ JCfrom additivity. Denoting , we get




          We get that for any given ti-tv ratio , Δ JC(R,t ) is a concave monotone increasing function, and its second derivative attains a global minimum of at . By the note following Lemma 2.8, the deviation of Δ JC from additivity in an interval [t 0,t 1] can be evaluated by computing the linear interpolation Δ int=At + B of Δ JC in [t 0,t 1], and finding t ∈[t 0,t 1] which maximizes Δ JC(t )−Δ int(t ) (see Figure 1a). A bound on this deviation from additivity can be obtained through Lemma 2.8 by plugging in the slope of the linear interpolation, A, and the maximum value, F, attained by the second derivative of Δ JC in [t 0,t 1]. Using Lemma 2.7 and an expression for dev (Δ JC(R,t ),[t 0,t 1]), it is possible to map out coherent collections of homogeneous K2P-trees for which Δ JC is guaranteed to be consistent. Each collection is defined by a range of ti-tv ratios [0.5,R max], a range of inter-leaf distances [t 0,t 1], and a lower bound on the weights of internal edges in the tree, given by t min=2dev (Δ JC(R max,t ),[t 0,t 1]).

          After determining a collection of trees for which a given non-affine-additive SR function, Δ, is consistent, one can compare the performance of Δ with additive alternatives. In our case, we compare Δ =Δ JC, which is not affine additive when , to the standard additive SR function Δ K2P. The potential advantage of Δ JCover Δ K2P lies in its reduced stochastic noise . Informally, this occurs because JC relies on the accuracy of estimating a single parameter - the sum p =p α + 2p β , while Δ K2P relies on the accuracy of estimating each of the two parameters p α and p β separately. The stochastic noise of an SR function is measured by the standard deviation of the statistical estimator associated with it, denoted σ (Δ JC) and σ (Δ K2P), respectively. We use the result in [9] to get a first order approximation (based on the delta method [29]) of σ (Δ K2P) for sequences of length k and model parameters R t :

          By a similar application of the delta method to Δ JC, we obtain:


          where k is the sequence length and (see (3)).

          Figure 1 provides an illustrative comparison of Δ JCand Δ K2P under the homogeneous sub-model of K2P with ti-tv ratio R =10, and within the inter-leaf time interval of [0.8,2]. Figure 1a shows the deviation of Δ JCfrom additivity in that setting, using its linear interpolation Δ int=At + B . Note that Lemma 2.8 and the subsequent note imply that , where X =maxt ∈[0.8,2]{Δ JC(t )−Δ int(t )}. Figure 1b depicts Δ JC in the same setting with its stochastic error margins (Δ JC±σ (Δ JC)), alongside its closest affine-additive function and its stochastic error margins (Δ ±σ (Δ )). These stochastic error margins are determined by assuming a sequence length of 500 bp in the first-order approximations given in (14) and (15), where σ (Δ ) is given by scaling σ (Δ K2P) by the slope A of the linear interpolation. Note how the margins of Δ JCare actually more tightly concentrated around its affine-additive approximation Δ than the margins of Δ . This implies that, despite its deviation from additivity in this setting, distances obtained using Δ JCare actually more likely to be near-additive than distances obtained using Δ K2P.

          Performance of Non affine-additive SR functions in quartet resolution

          The quartet tree is the smallest phylogenetic tree with non-trivial topology. Focusing on quartets enables a close study of the effects of deviation from additivity and stochastic noise on reconstruction accuracy. The topology of a quartet spanning four taxa {1,2,3,4} can be represented by the split notation (ij |kl ) (where {i j k l }={1,2,3,4}), indicating that the internal edge of the quartet separates i j from k l . All distance based quartet resolution algorithms essentially reduce to the four-point method (FPM) [26, 30], which resolves this split using the six observed pairwise distances : it first partitions the six observed distances into three sums , , and , and then determines the quartet split according to the minimal sum (the sum corresponds to the split (ij |kl )). We will focus on the task of reconstructing homogeneous K2P quartets using FPM with distances estimated using either Δ JCor Δ K2P. We note that most of our findings easily generalize to more sophisticated homogeneous substitution models, replacing Δ JC by any concave distance function and Δ K2P by some SR function corresponding to the evolutionary time t .

          For concreteness, we assume henceforth that the quartet split is (12|34), meaning that the sum of the exact evolutionary times t 12 + t 34 is minimal. We start by analyzing the impact of the deviation from additivity of Δ JC on the consistency of quartet resolutions. First, observe that any monotone distance function is consistent for quartets in which t 12 and t 34 are the smallest interleaf distances - as is the case with symmetric quartets, in which all external edges are of the same length. Therefore, we study two prototypes of asymmetric quartets. The length of the internal edge in both types is t i , and each type has two long external edges of length t l , and two short external edges of length t s . In type A quartets (Figure 2a), the short edges are on one side of the split and the long edges are on the other side. In this case d 12and d 34 re the smallest and largest interleaf distances (resp.). Hence, the concavity of Δ JC increases the separation between the sum d 12 + d 34 and the other two competing sums, leading to an expected improvement in reconstruction accuracy. The other quartet configuration (type B; Figure 2b) has a short edge and a long edge on both sides of the split. In this case, the interval of interpolation is [d 13,d 24], and the distance d 12=d 34 is near the center of this interval. Thus the concavity of Δ JC decreases the separation between the sums d 13 + d 24 and d 12 + d 34 by approximately twice the deviation from additivity of Δ JC in that range.

          Figure 2

          Performance of the Four Point Method using Δ JC on K2P quartets with ti-tv ratio R = 2. The concave non affine-additive SR function Δ JC is shown (dashed green line) in the interval [t 0,t 1], where t 0 and t 1 are the smallest and largest of the six pairwise distances (resp.). The dashed blue line shows the linear interpolation Δ int=At + B of Δ JC in the interval [t 0,t 1]. Horizontal dotted lines correspond to half of the two competing sums computed by FPM under the two SR functions (see legend). (a) In quartets of type A, t 0=t 12 and t 1=t 34, and so Δ int(1,2) + Δ int(3,4)=Δ JC(1,2) + Δ JC(3,4). However, for i ∈{1,2} and j ∈{3,4}, Δ int(i,j )<Δ JC(i,j ). Therefore, the deviation from additivity of Δ JC increases its FPM separation, denoted SEPJC, compared to the FPM separation SEPint of Δ int. (b) In quartets of type B, t 0=t 13 and t 1=t 24, and so Δ int(1,3) + Δ int(2,4)=Δ JC(1,3) + Δ JC(2,4). However, Δ int(1,2)=Δ int(3,4)<Δ JC(1,2)=Δ JC(3,4), and so Δ int(1,2) + Δ int(3,4)<Δ JC(1,2) + Δ JC(3,4). Therefore, the deviation from additivity of Δ JC decreases its FPM separation, denoted SEPJC, compared to the FPM separation SEPint of Δ int. Note that SEPint remains invariant in both types of quartets under fixed t i whereas SEPJC changes, depending on the type of quartet and the t s /t l ratio.

          When the deviation from additivity exceeds half the length of the internal edge, the sum d 13 + d 24 becomes the minimal sum, and Δ JC becomes inconsistent. Note that this demonstrates the tightness of the condition stated in Lemma 2.7, and in this sense, type B quartets provide a worst case scenario for quartet resolution by a concave SR functione.

          Next we turn to compare the accuracy of Δ JC with that of Δ K2P when used to reconstruct its “worst case scenario” quartets of type B. Interestingly, Δ JC ends up outperforming Δ K2P on many of these quartets, due to its reduced stochastic noise (as predicted in our discussion revolving around Figure 1b). For example, consider a series of homogeneous K2P quartets of type B with ti-tv ratio R =5, whose edge lengths were set as follows: t i =0.2, t l =1.0, and t s ∈[0.2,1.0]. We assessed reconstruction accuracy for both SR functions (Δ JCand Δ K2P) across this series of quartets, by generating 100,000 simulations of the substitution process using 1,000 bp long sequences for each quartet (Figure 3a). Despite its deviation from additivity, Δ JC outperforms the additive SR function Δ K2P on many of these quartets (as long as t l /t s <3.6) . Note that as t s shrinks, the deviation of Δ JC from additivity increases, since the interval [t 0,t 1] expands. This experiment appears to indicate that the deviation of Δ JC from additivity has to be quite large for Δ K2P to outperform it.

          Figure 3

          Performance of Δ JC and Δ K2P on a series of quartets of type B. A series of homogeneous K2P quartets is considered (left illustration), with ti-tv ratio of R =5, and edge lengths t i =0.2, t l =1, and t s ∈[0.2,1]. (a) Reconstruction accuracy using FPM and either Δ JC (dashed green) or Δ K2P (solid red) plotted against t l /t s . Accuracy ratio is estimated using 100,000 independent replicates for each value of t s in the interval [0.2,1] (in steps of 0.01), with sequence length 1,000 bp. (b) Fisher’s Criterion (FC) for the sums corresponding to splits (12|34) and (13|24) under either Δ JC (dashed green) or Δ K2P (solid red) plotted against t l /t s .

          Fisher’s criterion for separability

          We now present a simple and general framework based on the so-called Fisher Criterion (FC) for predicting the relative accuracy of two SR functions in resolving quartets. FC measures the effective separation between normal random variables Xn (μ 1 σ 1) and Yn (μ 2 σ 2) using the following measuref ( [20, 21]):


          We use FC to measure the separability of the distance sum corresponding to the true split (which should be the minimal sum for consistent SR functions) from the two remaining sums. For the expectation μ of each sum we use the true distances as computed by the SR function on the actual model parameters. For the variance σ 2, we use the sum of the approximate variances of the two distances involved in the sum. We expect that an SR function which provides a larger separation of the smallest sum from the two other sums will imply a better reconstruction probability.

          We note that FC is not an exact indicator of the separability in our case, because the necessary criteria for this are not satisfied in our model. Namely, the two distance sums are not normally distributed, and they are correlated through the substitution process along the external edges of the quartet. Nevertheless, as Figure 3b suggests, FC turns out to provide a quite reliable comparison of the expected performance of Δ JC and Δ K2P for the quartet series considered in the aforementioned experiment. Figure 3b exhibits for each quartet the FC of Δ JC alongside that of Δ K2P, both associated with the comparison of the true split (12|34) and the “Δ JC favored split” (13|24). As shown, the trends observed in both FC plots closely resemble the trends observed in the reconstruction accuracy plot (Figure 3a), and the the equilibrium point of the FC values of Δ JC and Δ K2P is very close to the equilibrium point of the accuracy of reconstructions of these two functions (near t l /t s =3.6).

          A useful feature of this framework is the natural way in which it teases apart the stochastic noise from the deviation from additivity. If we denote the numerator of FC by SEP (for “separation”) and its denominator by NOISE, then a comparison of FC estimates between two SR function Δ 1,Δ 2 can be represented as a ratio of ratios:

          Figure 4 illustrates how a comparison between the expected performance of Δ JC and that of Δ K2P can be carried out by tracing the SEP and NOISE ratios along four series of homogeneous K2P quartet: the bottom-left plot corresponds to the quartet series considered in Figure 3; the plot above it corresponds to the same series with ti-tv ratio R =2; the two plots on the right describe two quartet series in which the weight of the short edges is constant t s =0.2, and the weight of the long edges ranges in [0.2,1]. These four series demonstrate several typical trends in the behavior of the SEP and NOISE ratios. First, we observe that the NOISE ratio decreases (favoring Δ JC) as the diameter of the quartet (t 24) increases (it is almost constant in the two series on the left, and monotone decreasing in the series on the right). This is because the diameter provides the major contribution to the stochastic noise (for both Δ JC and Δ K2P), and as it increases, the ratio between the stochastic noise of Δ K2P and Δ JC increases as well. We also observe a natural decrease in the NOISE ratio with an increase in the ti-tv ratio (the NOISE ratio for R =5 is consistently smaller than for R =2). Concerning the SEP ratio, we see it becomes smaller (favoring Δ K2P) as the quartet becomes more unbalanced (the SEP ratio decreases along the X axis in each of the four plots). This is because the deviation of Δ JC from additivity increases as the inter-leaf distance interval [t 0,t 1]=[t 13,t 24] expands. Deviation of Δ JC from additivity also increases with the ti-tv ratio, as the substitution model further departs from the assumptions of JC (the SEP ratio for R =5 is consistently smaller than for R =2).

          Figure 4

          SEP and NOISE ratios. SEP (Δ JC)/SEP (Δ K2P) (dashed) and NOISE (Δ JC)/NOISE (Δ K2P) (dotted) plotted against t l /t s for four series of homogeneous K2P quartets of type B. Top two series have ti-tv ratio of R =2, and bottom two series have ti-tv ratio of R =5. Left two series have external edge lengths t l =1 and t s ∈[0.2,1], and right two series have external edge lengths t l ∈[0.2,1] and t s =0.2. The length of the internal edge is constant t i =0.2 in all four series.

          The two series on the right side of Figure 4 demonstrate well the tradeoff between the effects of stochastic noise and deviation from additivity. In both series, the SEP and NOISE ratios decrease as the quartets become more unbalanced (due to the trends listed above). However, the rates of decrease of these two ratios are different due to the different ti-tv ratios, and this determines the expected relative performance of the two SR functions across the series. When R = 2, the SEP ratio decreases at a slower rate than the NOISE ratio, and Δ JC is expected to outperform Δ K2P across the entire series. When R = 5, the SEP ratio decreases at a faster rate than the NOISE ratio, and when the quartets are sufficiently unbalanced (t l /t s >4) Δ K2P is expected to outperform Δ JC.

          Simulations on Hasegawa’s Tree

          In this section we describe experiments done on simulated data sets generated along the seven-taxon tree assembled by Hasegawa, Kishino, and Yano in 1985 [1, 6]. This tree, spanning seven eutherian mammals (Figure 5a), was reconstructed originally using mitochondrial DNA sequences. It has a caterpillar topology (meaning that every internal node is incident to an external edge), and it has long external edges and short internal edges, making it a suitable representative of small phylogenetic trees spanning moderately distant species. These features also make it particularly challenging for distance-based reconstruction.

          Figure 5

          Simulations on Hasegawa’s Tree. (a) Reconstruction accuracy of four different SR functions on different scaled versions of Hasegawa’s tree [6]. The tree with scaled edge weights is depicted (left) next to the graph (right) plotting reconstruction accuracy of four SR functions. Different scales of the tree are considered, indicated by the diameter of the tree (X axis). Reconstruction accuracy (Y axis) is measured for each scaled tree by the average normalized RF distance between the reconstructed tree and the true tree across 10,000 simulated data sets. Simulations were carried out assuming a ti=tv ratio of R = 2 and sequence length of 500 bp. (b) A similar plot is shown for a semi-symmetric caterpillar tree.

          In our study we used the tree structure and edge lengths to generate simulated data sets. We considered the tree in various scales, by setting the tree diameter (largest inter-taxon path length) to values in the interval [0.1,2.0]. For each scale considered, 10,000 simulations were carried out, where in each simulation 500 bp sequences were evolved along the tree according to a homogeneous K2P substitution model with ti-tv ratio of R =2. For each simulated data set, estimated values of the K2P statistics p α and p β , denoted by and , were extracted for all pairs of taxa. Subsequently, several distance matrices were computed for each data set by applying different SR functions to these estimated statistics. Reconstruction accuracy was evaluated by applying the Neighbor Joining (NJ) algorithm [31, 32] to these distance matrices and recording the Robinson-Foulds topological distance (RF) [33] between the reconstructed tree and the Hasegawa tree. Sequence simulation was performed using SeqGen [34] (by choosing the HKY model with uniform base frequencies), and tree reconstruction was performed using the version of NJ implemented in the PHYLIP package [35].

          We studied the reconstruction accuracy associated with four different SR functions: Δ JC, Δ K2P, Δ tv, and Δ R=2. The first two are as described in Equations (5) and (4), respectively. The third SR function, Δ tv, considers only tv-type substitutions: , and the fourth SR function, Δ R=2, is based on a maximum likelihood (ML) estimatorg of the time t from the estimated transition probabilities , given that R =2. Informally, this function, which uses knowledge of the true value of R (which is typically unknown to the user), is optimal in our setting, because it has similar stochastic noise as Δ JC, and it is additive since it coincides with Δ K2P when applied to transition probabilities that are consistent with a ti-tv ratio of R =2.

          The performance of these four SR functions is traced across the different tree scales in Figure 5a. For each SR function Δ and scale s, we recorded the average normalized RF distance from the true tree to each of the 10,000 trees reconstructed using Δ . The RF distance was normalized by its maximum value which is twice the number of internal edges in the tree (in our case 2×4=8). As observed previously in [9], Δ K2P performed well in shorter scales, and Δ tv performed well in longer scales. However, both additive SR functions were significantly outperformed in nearly all cases by Δ JC. Surprisingly, Δ JC even slightly outperformed Δ R=2. We speculate that this happened due to a bias similar to the one observed in type A quartets in Section Performance of Non affine-additive SR functions in quartet resolution, improving the performance of concave SR functions such as Δ JC on certain K2P-trees.

          To test this hypothesis, we went through a similar experiment with a more symmetric seven-taxon caterpillar tree, with internal edges of uniform length t int , and external edges of uniform length t ext =5t int (Figure 5b). The symmetry of this tree was expected to reduce the effect of the reconstruction bias observed in Hasegawa’s tree, and indeed, Δ JC performed much more poorly on this tree. Despite this fact, Δ JC still outperformed Δ K2P in all scales and Δ tv in the smaller scales (s <1.1).

          Inferring trees from genomic sequences

          In this section we describe our study comparing various SR functions on genomic DNA sequences. Next to Δ JC and Δ K2P we also considered the well known LogDet SR function [36, 37], denoted here as Δ LogDet. Extending our study to this setting is challenging in two respects. First of all, unlike the simulated case, the true tree is not known with complete confidence, and accuracy of reconstruction can only be determined by using a well-accepted reference tree that may contain some errors. Secondly, the true substitution model is also unknown and is likely to violate the assumptions of both JC and K2P models and even the relaxed assumptions of the general time-reversible model (in which Δ LogDet is additive). Hence, we have to assume in this case that Δ JC, Δ K2P, and Δ LogDet are all non affine-additive, where Δ JC and Δ K2P are still likely to exhibit higher deviation from additivity than Δ LogDet, since they make stronger assumptions on the substitution model.

          The genomic data set

          In building the genomic data set, we made use of a set of 31 clusters of orthologous groups (COGs) which was compiled by Ciccarelli et al. and used for inferring phylogenetic relationships amongst a large number of species in [38, 39]. These 31 gene families were selected to capture the evolutionary history of the species containing them. This was done in [38] by making sure that the genes in these families have the following properties: (1) they are highly conserved across species, (2) they have a small number of paralogs, and (3) they are weakly affected by horizontal gene transfer. We scanned the NCBI genome database and found 199 bacterial genomes that contained all annotated COGs. For each of the 31 COGs, we extracted the appropriate protein sequence in each of the 199 bacterial species, choosing an arbitrary paralog in cases of multiple hits. We followed a procedure similar to the one described in [38, 39] to obtain reliable multiple-sequence alignments for each COG: we computed a 199-way multiple alignment of the protein sequences of each COG using HMMalign [40] and then mapped each protein sequence back to its coding DNA sequence. The conserved parts of each of the 31 DNA alignments were extracted using GBLOCKS [41] to filter out alignment columns with 50% or more gap symbols. The alignments were manually scanned, and 36 species which contributed a large number of gaps to the alignments were removed from the subsequent analysis. The 31 different alignments were concatenated to form one long 163-way multiple sequence DNA alignment.

          For the reference tree we used the phylogenetic tree of microbial species provided by the ARB-SILVA Living Tree Project [42]. This tree, spanning 8,029 species at the time of writing, is based on a widely accepted analysis of the small subunit (SSU) 16S RNA. A subtree spanning our 163 bacterial species was extracted from this tree and treated as the true phylogenetic tree in our analysis.

          Reconstruction accuracy for ten-species subsets

          We used the base set of 163 species to generate 40,000 random 10-species sub-alignments. The random selection process was guided to generate species subsets corresponding to a wide range of diameter scales (a blind random selection process is biased toward subsets with large diameters). For each of the 40,000 subsets, a 10-way subalignment was extracted from the original 163-way alignment, and in this alignment we extracted only columns corresponding to four-fold degenerate sites that do not have any gap symbol. This is done to make sure the sites used for distance estimation have undergone a substitution process that is as uniform as possible along the different lineages and across the different sites. Each sub-alignment was used to compute three distance matrices – one under Δ JC, one under Δ K2P, and one under Δ LogDet. The latter was calculated by the version that is implemented in the PHYLIP package. The NJ algorithm was then applied to the three matrices and the resulting trees were compared to the true tree (as depicted by the appropriate LTP subtree) according to the RF distance.

          As an additional comparison, we used a fourth reconstruction technique. This method (termed BIONJ-GTR) used the BIONJ reconstruction algorithm [43] on distances obtained under the general time-reversible model with invariant sites and Gamma distribution of rates across variant sites (GTR+Γ +I) [8, 44].

          The PhyML package [45] was used to infer this tree for each of the 40,000 subsets. We selected the GTR+Γ +I model since it was found by the MEGA5 software [46] to provide the best fit to the sequence data. The 40,000 sampled instances were partitioned into eight bins according to the RF distance observed between the BIONJ-GTR tree and the true (LTP) tree, and average RF distances were recorded for each of the three SR functions in each bin. This allowed us to observe trends throughout these 40,000 samples (Figure 6). Of the 40,000 trees inferred under Δ JC, 83.1% showed an equal or lower RF distance than those reconstructed by the BIONJ-GTR method. Moreover, Δ JC outperformed Δ K2P and Δ LogDet on average in all partitions, and Δ LogDet showed by far the worst performance with 48.7% of all reconstructed trees achieving higher RF distances to the reference tree than those inferred by BIONJ-GTR. As with our results on simulated data sets, we see that the SR functions with lower stochastic error but inferior model fit performed best. Unsurprisingly, the GTR+G+I model itself, which was predicted to have the best fit to the sequence data, was often outperformed by the simpler JC and K2P models. Note that the difference in performance between Δ JC and the two other SR functions is greater for subsets that are more accurately reconstructed by the BIONJ-GTR approach (the lower bins). This appears to indicate that over-simplified distance methods are particularly beneficial when the sequence data conveys a stronger phylogenetic signal.

          Figure 6

          Evaluation against BIONJ-GTR tree. The 40,000 subsets of size 10 were partitioned according to the the RF-distance between the reference LTP tree and the tree reconstructed using BIONJ-GTR (X axis). The (left) Y axis describes the mean difference between the RF-distance associated with a tree reconstructed using a particular SR function (Δ K2P, Δ JC, or Δ LogDet) and the RF-distance associated with the BIONJ-GTR tree. The bar plot in the background depicts the number of subsets in each bin.


          In this paper we explored the basic properties of methods for estimating evolutionary distances, and studied how these properties affect the accuracy of distance-based phylogenetic reconstruction. We considered both the systematic bias and the stochastic noise (variance) of the distance estimators, and examined the tradeoff between these two factors. We focused on the common task of phylogenetic reconstruction under homogeneous substitution models. Assuming homogeneous models simplifies the analytical framework, since in such models each SR function is reduced to a univariate function of the evolutionary time t . However, obtaining accurate estimates of t is still a hard task in this setting, since the unit rate matrix is unknown. An SR function Δ is guaranteed to yield consistent reconstruction across all trees in a homogeneous model only if it is additive, meaning that it is a linear function of t . When Δ is not additive, it introduces a systematic bias in distance estimates, which we denoted here as deviation from additivity . Some SR functions are only additive in one homogeneous model, whereas others are additive across a wider collection of homogeneous models. This less constrained additivity is typically achieved at a price of increased estimation noise. We studied the tradeoff between “deviation from additivity” and “estimation noise” via a case study where the model tree is a homogeneous K2P-tree with an unknown ti-tv ratio R . In this case, Kimura’s distance formula Δ K2P is always additive, while the less noisy Jukes Cantor’s formula, Δ JC, is additive only when .

          A study of this type requires a way to measure the deviation from additivity of a non-additive SR function Δ in a given range of distances [t 0,t 1]. To this end, we introduced the concept of affine-additive distance functions, and defined the deviation from additivity of Δ in [t 0,t 1] as the distance of Δ from its closest affine-additive function in [t 0,t 1]. We established a tight connection between this measure and statistical consistency of reconstruction (Lemma 2.7) and derived an upper bound for deviation from additivity in homogeneous models (Lemma 2.8). We applied these results in analyzing the deviation from additivity of Δ JC, and its effect on the accuracy of reconstructing homogeneous K2P-trees. We then showed, both analytically (in Section Deviation from additivity in homogeneous substitution models) and through experiments on simulated data sets (in Sections Performance of Non affine-additive SR functions in quartet resolution and Simulations on Hasegawa’s Tree), that, compared to Δ K2P, it is often better to use the non-additive but less noisy estimates of Δ JC, even when R is quite high. Somewhat surprisingly, we found this to be the case even when the tree being reconstructed has an “unfavorable” topology. Our experiments on bacterial gene sequences (Section Inferring trees from genomic sequences) also indicate that the simple and less noisy SR functions perform better on average than ones that are expected to better fit the true substitution process.

          The framework presented in this paper implies a practical way for selecting SR functions which are likely to increase the accuracy of distance estimation. The practicality of the method is drawn from the fact that the criteria by which we select an SR function depend only a relatively crude information about the tree being reconstructed. For instance, in the case of a homogeneous K2P-tree, one can easily obtain from the input sequences rough estimates of both the ti-tv ratio R and the range of inter-leaf times [t 0,t 1]. These estimates can then be used to compare the expected accuracies of Δ JC and Δ K2P on the given input, and determine which of them is more likely to yield an accurate phylogeny. For quartets, a tight comparison can be made using the FC-based approach suggested in Section Fisher’s Criterion for Separability, and for larger trees, a cruder comparison can be made using a plot like the one presented in Figure 1b. A promising avenue of further research is to extend the FC-based approach to allow tighter prediction of reconstruction accuracy of trees spanning more than four taxa.


          aThis is a WABI 2011 special issue invited paper. Extended abstract of this paper appeared in [47]. bTypically, the unit rate matrix is assumed to be the one corresponding to one substitution per site. cMany common distance-based algorithms, such as the Neighbor Joining (NJ) algorithm [31, 32], are known to be robust in this sense. dIn a tree, edges which touch leaves are external, and all other edges are internal. eTypes A and B quartets represent the Farris zone and Felsenstein zone, resp. (see, e.g., [1], Chapter 9). fWe use here the square root of the criterion commonly used in the literature, because we prefer to think in terms of distances rather than squares of distances. This has no practical influence, since we use FC only for comparing between different choices, not for assessing the quality of a give choice.g This ML estimate is obtained by a simple numerical method for maximizing the likelihood function (see, e.g., [1]).


          Tightness of Lemma 2.8

          Let f (t ) be a (continuous) function on some interval [t 0,t 1]. We prove below that if f does not intersect its linear interpolation At + B in that interval, then . We use the following notations, conforming to the notations in the proof of Lemma 2.8:

          Lemma 2.9

          Let f (t ) be a monotone increasing function in the interval [t 0,t 1] and let At + B be its linear interpolation in [t 0,t 1]. If either f (t )≥At + B for all t ∈[t 0,t 1] or f (t )≤At + B for all t ∈[t 0,t 1], then for all a >0, we have .


          We prove the minimality of in the case where f (t )≥At + B for all t ∈[t 0,t 1]. The other case (where f (t )≤At + B for all t ∈[t 0,t 1]) can be proven in an identical fashion.

          For a >0, let B a be the maximum value of B s.t. for all t ∈[t 0,t 1]. Evidently, . If the linear interpolation of f (t ) in [t 0,t 1] is given by At + B, then B A =B . We need to show that for every a >0, it holds that (a,B a )> (A,B A ). Let t A be a point in [t 0,t 1] s.t. ψ (A,B A ,t A )=ψ (A,B A ). Note that if a <A, then the two linear functions At + B A and at + B a intersect at (t 0,f (t 0)), and if a >A, then they intersect at (t 1,f (t 1)) (see Figure 7).

          Figure 7

          Proof of Lemma 2.9. A function f(t) is portrayed (bold) with its linear interpolation At + B =At + B A (green) in the interval [t 0 , t 1 ], s.t.f (t)≥At + B for allt ∈[t 0 ,t 1 ]. Equation (18) is illustrated for a <A on the right, and equation (19) is illustrated for a >A on the left.

          For a <A, we get the following equality (Figure 7; right):

          Hence, since ψ (a,B a )≥ψ (a,B a ,t ) for every t ∈[t 0,t 1], and since a <A, we get

          Similarly, if a >A, we get the following equality (Figure 7; left)

          and a >A implies that



          This research was supported by the Israel Science Foundation (ISF) grant No. 509/11. We also acknowledge the support for the publication fee by the Deutsche Forschungsgemeinschaft and the Open Access Publication Funds of Bielefeld University. The third author would also like to thanks the Max Planck Institute for Informatics for supporting a visit under which part of this research was carried out.

          Authors’ Affiliations

          Center for Biotechnology, Bielefeld University
          Department of Biological Statistics and Computational Biology, Cornell University
          Computer Science Department, Technion - Israel Institute of Technology


          1. Felsenstein J: Inferring Phylogenies. Sunderland: MA Sinauer Associated Inc; 2004.
          2. Semple C, Steel M: Phylogenetics. Oxford University Press; 2003.
          3. Papoulis A, Pillali SU: Probability, Random Variables and Stochastic Processes. New York: McGraw Hill Higher Education; 2002.
          4. Jukes T, Cantor C: Evolution of Protein Molecules. In Mammalian Protein Metab. Edited by: Munro H. New York: Academic Press; 1969:21–132.
          5. Kimura M: A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J Mol Evol 1980,16(2):111–120.PubMedView Article
          6. Hasegawa M, Kishino H, Yano T: Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J Mol Evol 1985,22(2):160–174.PubMedView Article
          7. Tavaré S: Some Probabilistic and Statistical Problems in the Analysis of DNA Sequences. Lectures on Mathematics in the Life Sci 1986, 17:57–86.
          8. Lanave C, Preparata G, Saccone C, Serio G: A new method for calculating evolutionary substitution rates. J Mol Evol 1984, 20:86–93.PubMedView Article
          9. Gronau I, Moran S, Yavneh I: Towards Optimal Distance Functions for Stochastic Substitution Models. J Theor Biol 2009,260(2):294–307.PubMedView Article
          10. Gronau I, Moran S, Yavneh I: Adaptive Distance Measures for Resolving K2P Quartets: Metric Separation versus Stochastic Noise. J Comp Biol 2010,17(11):1391–1400.View Article
          11. Felsenstein J: Cases in which parsimony or compatability methods will be positively misleading. Syst Zool 1978, 27:401–410.View Article
          12. Cavender J: Taxonomy with confidence. Math Biosci 1978, 40:271–280.View Article
          13. Steel M, Penny D: Parsimony, likelihood, and the role of models in molecular phylogenetics. Mol Biol Evol 2000, 17:839–850.PubMedView Article
          14. Sober E: A likelihood justification of parsimony. Cladistics 1985, 1:209–233.View Article
          15. Felstenstein J, Sober E: Parsimony and likelihood: an exchange. Syst Zool 1986, 35:617–626.View Article
          16. Yang Z: How often do wrong models produce better phylogenies? Mol Biol Evol 1997, 14:105–108.PubMedView Article
          17. Bruno WJ, Halpern AL: Topological bias and inconsistency of maximum likelihood using wrong models. [http://​www-t10.​lanl.​gov/​billb/​BrunoHalpern99.​pdf] Mol Biol Evol 1999,16(4):564–566.PubMedView Article
          18. Zharkikh A: Estimation of evolutionary distances between nucleotide sequences. J Mol Evol 1994,39(3):315–329.PubMedView Article
          19. Gascuel O, Guindon S: Efficient Biased Estimation of Evolutionary Distances When Substitution Rates Vary Across Sites. Mol Biol Evol 2002,19(4):534–543.PubMedView Article
          20. Fisher R: The use of multiple measurements in taxonomic problems. Ann of Eugenics 1936, 7:177–188.
          21. Duda R, Hart P: Pattern Classification and Scene Analysis. Hoboken: John Wiley and Sons; 1973.
          22. Sumner J, Fernandez-Sanchez J, Jarvis P: Lie Markov Models. J Theor Biol 2012, 298:16–31.PubMedView Article
          23. Buneman P: The recovery of trees from measures of dissimilarity. In Mathematics in the Archeological and Historical Sciences. Edited by: Hodson F, Kendall D, Tautu P. Edinburgh University Press; 1971:387–395.
          24. Sattath S, Tversky A: Additive similarity trees. Psychometrica 1977,42(3):319–345.View Article
          25. Atteson K: The Performance of Neighbor-Joining Methods of Phylogenetic Reconstruction. Algorithmica 1999, 25:251–278.View Article
          26. Erdos P, Steel M, Szekely L, Warnow T: A few logs suffice to build (almost) all trees (I). Random Struct Algorithms 1999, 14:153–184.View Article
          27. Erdos P, Steel M, Szekely L, Warnow T: A few logs suffice to build (almost) all trees (II). Theoret Comput Sci 1999, 221:77–118.View Article
          28. Johnson L, Riess R: Numerical Analysis. Boston: Addison Wesley; 1977.
          29. Oehlert G: A note on the delta method. Am Statistician 1992, 46:27–29.
          30. Zaretskii K: Constructing a tree on the basis of a set of distances between the hanging vertices. Uspekhi Mat Nauk 1965,20(6):90–92. [In Russian]
          31. Saitou N, Nei M: The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol 1987, 4:406–425.PubMed
          32. Studier J, Keppler K: A note on the neighbor-joining algorithm of Saitou and Nei. Mol Biol Evol 1988,5(6):729–731.PubMed
          33. Robinson F, Foulds R: Comparison of phylogenetic trees. Math Biosci 1981, 53:131–147.View Article
          34. Rambaut A, Grass NC: Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees. Comput Appl Biosci 1997,13(3):235–238.PubMed
          35. Felsenstein J: PHYLIP - Phylogeny Inference Package (Version 3.2). Cladistics 1989, 5:164–166.
          36. Steel M: Recovering a tree from the leaf colourations it generates under a Markov model. Appl Math Lett 1994,7(2):19–24.View Article
          37. Lockhart P, Steel M, Hendy M, Penny D: Recovering evolutionary trees under a more realistic model of sequence evolution. Mol Biol Evol 1994,11(4):605–612.PubMed
          38. Ciccarelli FD, Doerks T, von Mering C, Creevey CJ, Snel B, Bork P: Toward Automatic Reconstruction of a Highly Resolved Tree of Life. Science 2006,311(5765):1283–1287.PubMedView Article
          39. von Mering C, Hugenholtz P, Raes J, Tringe SG, Doerks T, Jensen LJ, Ward N, Bork P: Quantitative Phylogenetic Assessment of Microbial Communities in Diverse Environments. Science 2007,315(5815):1126–1130.PubMedView Article
          40. Durbin R, Eddy SR, Krogh A, Mitchison G: Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge University Press; 1999.
          41. Talavera G, Castresana J: Improvement of phylogenies after removing divergent and ambiguously aligned blocks from protein sequence alignments. Syst Bio 2007, 56:564–577.View Article
          42. Yarza P, Ludwig W, Euzeby J, Amann R, Schleifer KH, Glockner FO, Rossello-Mora R: Update of the All-Species Living Tree Project based on 16S and 23S rRNA sequence analyses. Syst Appl Microbiol 2010, 33:291–299.PubMedView Article
          43. Gascuel O: BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data. Mol Biol Evol 1997,14(7):685–695.PubMedView Article
          44. Rodriguez F, Oliver JL, Marin A, Medina JR: The general stochastic model of nucleotide substitution. J Theor Biol 1990, 142:485–501.PubMedView Article
          45. Guindon S, Gascuel O: A simple, fast and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol 2003, 52:696–704.PubMedView Article
          46. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S: MEGA5: Molecular Evolutionary Genetics Analysis using Maximum Likelihood, Evolutionary Distance, and Maximum Parsimony Methods. Mol Biol Evol 2011, 28:2731–2739.PubMedView Article
          47. Doerr D, Gronau I, Moran S, Yavneh I: Stochastic Errors vs. Modeling Errors in Distance Based Phylogenetic Reconstructions. In Algorithms in Bioinformatics, Volume 6833 of Lecture Notes in Computer Science. Edited by: Przytycka T, Sagot MF. Berlin /Heidelberg: Springer; 2011:49–60.

          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.