Stochastic errors vs. modeling errors in distance based phylogenetic reconstructions

Background 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. Results 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. Conclusions 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.


Introduction
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 [1][2][3]. Distance-based reconstruction algorithms tackle this task by first computing a set of n 2 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 *Correspondence: moran@cs.technion.ac.il 3 Computer Science Department, Technion -Israel Institute of Technology, Haifa, Israel Full list of author information is available at the end of the article 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. http://www.almob.org/content/7/1 /22 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 paper a , 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 [13][14][15]. 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, nonadditive 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.

Background
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. http://www.almob.org/content/7/1/22

Substitution Models
In this work, a DNA substitution model M is simply a set of stochastic 4 × 4 transition matrices closed under matrix product (i.e., P, Q ∈ M → PQ ∈ M). 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 M, or an M-tree, is an undirected tree T = (V , E) in which each edge e ∈ E is associated with a transition matrix P e ∈ M. An M-tree T implies an inter-leaf transition matrix P ij ∈ M for each pair of leaves {i, j} ⊂ L(T), namely P ij = e∈path T (i,j) P e . 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: 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: M R = {e tR : t > 0}. 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 M R 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 (A ↔ G, C ↔ T), and β, which is the rate of transversion-type (tv) substitutions ({A,G} ↔ {C,T}). 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 R = α 2β ≥ 1 2 . The Jukes-Cantor (JC) model [4] is a special homogeneous sub-model of K2P, in which R = 1 2 (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: (2) (3)

Substitution rate functions
A substitution rate (SR) function for a model M is a nonnegative continuous function : M → R + 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 M-tree T: Of particular interest in phylogenetic reconstruction are additive SR functions. 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. http://www.almob.org/content/7/1/22 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]. It is well known that additive SR functions imply additive metrics: if is an additive SR function for a model M, then for any M-tree T, D 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 Tadditive 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 algorithms c .
Formally, "sufficiently close" is defined by the following relation: where w min (D ) is the minimal weight assigned to an internal edge d 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 affineadditive mapping, defined below.

Definition 2.4 (Affine-additive mapping). A dissimilarity mapping D is said to be affine-additive w.r.t. a phylogenetic tree T, if there is a T-additive metric D, and scalars
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 w (·) corresponding to the affine additive mapping D = aD + b is given by: w (e) = aw(e) for all internal edges, and w (e) = aw(e) + 1 2 b for all external edges. When b is positive, D is actually an additive metric, but when b is negative, the weights of external edges implied by w (·) might be negative, and D 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 M, as follows: 1. If is additive for M, then for each M-tree T the mapping D T defined by D T (i, j) = (P ij ) for all i, j ∈ L(T), is a T -additive metric. http://www.almob.org/content/7/1/22

As the length of the input sequences grows, the estimated transition matrices
is sufficiently close to D T , and is thus near-additive. 4. The near-additivity of the estimated dissimilarity map D 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., [25][26][27]), and in all these cases an additive SR function is assumed. Notice, however, that this line of argument remains valid when D T is near additive w.r.t. T. For instance, consistent reconstruction of any M-tree is guaranteed by using an affine-additive SR function , which is an affine transformation of some additive SR function : = a + b (with a > 0). An SR function that is not affine-additive in a given substitution model M does not guarantee consistency across all M-trees, but it still can be consistent for specific M-trees.

Definition 2.5 (Consistent SR function). An SR function of a substitution model M is said to be consistent w.r.t. an M-tree T if D T 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 M, then it might be consistent with respect to many M-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 . 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 M R . The unit rate matrix R implies a 1-1 mapping between evolutionary time t and rate matrices in M R . It is thus useful to view an SR function for M R as a function : R + → R + 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 a ∈ R + , b ∈ R. We define the deviation of an SR function from a given affine-additive function at + b in an interval [ t 0 , t 1 ] as 1 a max{| (t)−at−b| : t ∈[ t 0 , t 1 ] } (the factor 1 a 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
: 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.
Proof. We need to show that D T is near-additive w.r.t.
For all i, j ∈ L(T), denote t ij = e∈path T (i,j) t e , 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 D = aD + b is an affine-additive mapping. The internal-edge-weights associated with D are given by w (e) = at(e) (see discussion following Definition 2.4), implying that w min (D ) = at min . 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 inter- and http://www.almob.org/content/7/1/22 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.

Lemma 2.8. Let
: R + → R + be an SR function in a homogeneous substitution model, and Proof. Let us start by introducing a couple of auxiliary notations: We are looking for a ∈ R + and b ∈ R which mini- 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 . Therefore, the bound on the error of polynomial interpolation (see, e.g., [28], p. 187) implies that Note. In Appendix 3 we prove that if does not intersect its linear interpolation int = At + B within the interval (t 0 , t 1 ), then the function At+b * mentioned in the proof above is, in fact, the affine-additive function which minimizes the deviation from additivity of in [ t 0 , t 1 ]. This means that, in such cases, the first inequality in (9) holds in equality. The last inequality in (9) also holds in equality in such cases, because we are guaranteed to have either [ t 2 , t 3 ] =[ t 0 , t 1 ] (when is bounded from above by its linear interpolation) or [ t 4 , t 5 ] =[ t 0 , t 1 ] (when is bounded from below by its linear interpolation). Thus, in such a case, the bound of Lemma 2.8 is reduced to the bound on interpolation error (middle inequality in (9)). Cases where does not intersect its linear interpolation are frequent among many SR functions of interest, as this condition holds when is either convex or concave.

Deviation of JC from Additivity in K2P
We now turn to study the deviation of JC from additivity in homogeneous sub-models of K2P with ti-tv ratio R > 1 2 . First, we express JC as a function of the ti-tv ratio R and the time t, using (5) and the relations α 2β = R and α + 2β = 1.
Note that the homogeneous K2P sub-model with R = 1 2 is the JC model; in this case the second term of (10) vanishes, leaving JC ( 1 2 , t) = t. For other homogeneous sub-models of K2P, where R > 1 2 , JC is not affineadditive (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 JC from additivity. Denoting ρ = 2R−1 R+1 , we get The affine-additive SR function minimizing its deviation from JC is * = int + 1 2 X. The two SR functions JC and * are shown with their stochastic error margins, assuming sequence length of 500 bp.
We get that for any given ti-tv ratio R > 1 2 , JC (R, t) is a concave monotone increasing function, and its second derivative attains a global minimum of − 3 16 ρ 2 at t = ln (2) ρ . 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 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 R > 1 2 , to the standard additive SR function K2P . The potential advantage of JC over 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 p(t, (3)). Figure 1 provides an illustrative comparison of JC and K2P under the homogeneous sub-model of K2P with titv ratio R = 10, and within the inter-leaf time interval of [ 0. 8,2]. Figure 1a shows the deviation of JC from additivity in that setting, using its linear interpolation int = At + B. Note that Lemma 2.8 and the subsequent note imply that dev( JC , [ 0. Figure 1b depicts JC in the same setting with its stochastic error margins ( JC ± σ ( JC )), alongside its closest affine-additive function * = int + 1 2 X 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 JC are actually more tightly concentrated around its affine-additive approximation * than the margins of * . This implies that, http://www.almob.org/content/7/1/22 despite its deviation from additivity in this setting, distances obtained using JC are 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 { d ij : {i, j} ⊂ {1, 2, 3, 4}}: it first partitions the six observed distances into three sums d 12 + d 34 , d 13 + d 24 , and d 14 + d 23 , and then determines the quartet split according to the minimal sum (the sum d ij + d kl corresponds to the split (ij|kl)). We will focus on the task of reconstructing homogeneous K2P quartets using FPM with distances { d ij } estimated using either JC or 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 12 and d 34 are 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.
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 function e .
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 ( JC and 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.

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 X ∼ N(μ 1 , σ 1 ) and Y ∼ N(μ 2 , σ 2 ) using the following measure f ( [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 http://www.almob.org/content/7/1/22  (3, 4). Therefore, the deviation from additivity of JC decreases its FPM separation, denoted SEP JC , compared to the FPM separation SEP int of int . Note that SEP int remains invariant in both types of quartets under fixed t i whereas SEP JC changes, depending on the type of quartet and the t s /t l ratio.
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 http://www.almob.org/content/7/1/22 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 titv 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).
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.
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 intertaxon 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 byp α andp β , were extracted for all 7 2 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. http://www.almob.org/content/7/1/22 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: tv (p α , p β ) = − 1 4 log 1 − 4p β (t) = βt, and the fourth SR function, R=2 , is based on a maximum likelihood (ML) estimator g of the time t from the estimated transition probabilitiesp α ,p β , 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 probabilitiesp α ,p β 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 affineadditive 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, http://www.almob.org/content/7/1/22  [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.
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 http://www.almob.org/content/7/1/22 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 10way subalignment was extracted from the original 163way 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 oversimplified distance methods are particularly beneficial when the sequence data conveys a stronger phylogenetic signal.

Conclusions
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 http://www.almob.org/content/7/1/22 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.
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 R = 1 2 . 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 affineadditive 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  (18) is illustrated for a < A on the right, and equation (19) is illustrated for a > A on the left. http://www.almob.org/content/7/1/22 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.
Endnotes a This is a WABI 2011 special issue invited paper. Extended abstract of this paper appeared in [47]. b Typically, the unit rate matrix is assumed to be the one corresponding to one substitution per site. c Many common distance-based algorithms, such as the Neighbor Joining (NJ) algorithm [31,32], are known to be robust in this sense. d In a tree, edges which touch leaves are external, and all other edges are internal. e Types A and B quartets represent the Farris zone and Felsenstein zone, resp. (see, e.g., [1], Chapter 9). f We 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]). http://www.almob.org/content/7/1/22 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.