How Fitch-Margoliash Algorithm can Benefit from Multi Dimensional Scaling.

Whatever the phylogenetic method, genetic sequences are often described as strings of characters, thus molecular sequences can be viewed as elements of a multi-dimensional space. As a consequence, studying motion in this space (ie, the evolutionary process) must deal with the amazing features of high-dimensional spaces like concentration of measured phenomenon.TO STUDY HOW THESE FEATURES MIGHT INFLUENCE PHYLOGENY RECONSTRUCTIONS, WE EXAMINED A PARTICULAR POPULAR
METHOD: the Fitch-Margoliash algorithm, which belongs to the Least Squares methods. We show that the Least Squares methods are closely related to Multi Dimensional Scaling. Indeed, criteria for Fitch-Margoliash and Sammon's mapping are somewhat similar. However, the prolific research in Multi Dimensional Scaling has definitely allowed outclassing Sammon's mapping.Least Square methods for tree reconstruction can now take advantage of these improvements. However, "false neighborhood" and "tears" are the two main risks in dimensionality reduction field: "false neighborhood" corresponds to a widely separated data in the original space that are found close in representation space, and neighbor data that are displayed in remote positions constitute a "tear". To address this problem, we took advantage of the concepts of "continuity" and "trustworthiness" in the tree reconstruction field, which limit the risk of "false neighborhood" and "tears". We also point out the concentration of measured phenomenon as a source of error and introduce here new criteria to build phylogenies with improved preservation of distances and robustness.The authors and the Evolutionary Bioinformatics Journal dedicate this article to the memory of Professor W.M. Fitch (1929-2011).

Introduction

Phylogenetic trees building methods

Since the early works of Hitchcock and Hitchcock (1840),1 Darwin (1859, the only one illustration of the book)2 and those of Ernst Haeckel (1866),3 species are classified within trees. This idea has never been challenged and is supported by the evolutionary theory, which states that species genotypes/phenotypes drift from a common ancestor according to phylogenic proximity of the considered species. Trees are thereby convenient tools to express phylogenic relationships.

However, evolution has not been observed per se, and we must infer it from contemporary information (see Brocchieri4 for a review on phylogenetic methods). Phylogenetic trees are often achieved from aligned sequences using character-based methods. Maximum parsimony methods minimize the number of changes.5 Maximum likelihood methods are based on models of substitution process.6

Molecular character based methods can be launched from genomic or proteomic sequences but also from other molecular features like structural or thermodynamic features.7,8 Conversely, distance based methods are often less powerful,9 but they can be applied to any kind of matrices (of possibly large size) collecting distances between items (whatever the distances used—from aligned sequences to any measurable features). For example, UPGMA (Unweighted Pair Group Method with Arithmetic mean)10 iteratively aggregates the most similar items. UPGMA has been said to suffer of long-branch attraction due to the (frequently false) hypothesis of the molecular clock (every sequence is supposed to evolve at the similar rate). In order to solve this problem, Neighbor Joining (NJ)11 generates trees with the minimum cumulative branch lengths. We highly recommend the use of character based methods (such as maximum likelihood, minimum evolution, maximum parsimony, Bayesian approaches; they are more powerful9 when the related hypotheses are valid). However, in some difficult cases, distance based methods are the only option4,12 in large scale comparisons13 and analyses on oligonucleotides frequencies (the so-called “genomic signatures”),14,15 etc.

Note that every distance method assumes that a tree that preserves distances reflects the evolution between data. Fitch and Margoliash16 followed this concept and described a method designed to generate a tree that preserves distances at best. The tree that minimizes squared difference between original and resulting distances is assumed to express the evolution of considered genes (or proteins, if it is the case). Here, we hold on this widely shared hypothesis.

Biological sequences and amino acids are clearly multi-dimensional objects.17–20 It has been recently demonstrated that the high-dimensionality of biological sequences leads to emergent computational features like similarity measure and particular probability distribution of this similarity.21–23 Distance methods are the first approach for phylogeny reconstructions and consist merely in considering the data as a point in an η-dimensional phase space (where η is in first approximation the length of the sequence). If one knows the relative position of sequences (distance matrices) and the velocity modulus corresponding to the sequence evolution, it should be possible to reconstruct the divergence history of the sequences (which corresponds to a trajectory in the η-dimensional space). However, starting from a point of this phase space, another point cannot be reached by a simple straight line because only a few likely paths between the two points in this phase space are allowed (Biological features of sequences must be preserved all along the path).24 The topology of all a priori possible evolution process space has an E topology, where E is the space of the states of one residue position. Nevertheless, the topology of actual evolution process is clearly unknown, since the E subspace where the process occurs is also unknown. This is one of the reasons why maximum likelihood and parsimony methods are usually more accurate than distance methods. Indeed, ML and parsimony do take the intrinsic high-dimensional nature of the original data (ie, the sequences) into account, by modeling resp. optimizing character state changes. Conversely, distance-based methods such as Neighbor Joining or Fitch-Margoliash algorithm (which do not account for the data dimensionality) are commonly performed from genomic signature comparisons, very large data sets, such as DNA barcode libraries, or indeed largescale genomic comparisons, and so on. However, these often are high dimensional data and we will show in the following why Fitch-Margoliash algorithm could be misled by features of high dimensional space, and how modifying it to avoid these traps.

Multi dimensional scaling

The general aim of the Multi Dimensional Scaling methods is the dimensionality reduction. The most used technique is by far the Principal Components Analysis (or PCA).25,26 To obtain this representation, data are projected (by orthogonal projection) on a selected vector subspace according to a criterion equivalent to the maximization of distances preservation.27 The “Multi Dimensional Scaling” (MDS) term was extracted from a method that allowed a mapping from a distance matrix27 (this method is known thereafter as “Classical Multi Dimensional Scaling” or Classical MDS). When dealing with Euclidean distances, Classical MDS can be seen as a linear projection method. However, non-linear projection can be achieved as well. For that purpose, the preservation of small distances is favored through a weighting system. The Sammon’s mapping28 is one of the first non-linear MDS methods. Many MDS techniques have been developed since then, including Curvilinear Components Analysis (or CCA),29 Kernel Principal Component Analysis (or KPCA),30 Isomap,31 Curvilinear Distances Analysis (CDA)32 and Generative Topographic Mapping (GTM).33

We focus here on a new efficient method called Data-Driven High Dimensional Scaling (DD-HDS)34 which specificities are i) the penalization of “false-neighborhoods” and “tears” and also ii) the consideration of the concentration of measured phenomenon. As many MDS methods (including Sammon’s mapping, CCA and Local-MDS) DD-HDS is driven by a given criterion in order to iteratively converge toward a satisfying map (several examples of criteria are given in the beginning of section 3).

The present paper shows how close are Multi Dimensional Scaling and tree building methods fields and extracts from Multi Dimensional Scaling new ways to improve actual tree building methods from distance matrices. The article is organized as follows: Section 2 compares the Fitch-Margoliash algorithm with Multi Dimensional Scaling methods (with a special attention for Sammon’s mapping)28 Section 3 describes the possible contribution of the recently introduced Multi Dimensional Scaling methods to tree building methods through a new criterion. Section 4 sums up the resulting algorithm, used to generate results presented and analyzed in Section 5.

Existing links between Fitch-Margoliash and Sammon’s Mapping criteria

Least Square (LS) methods are designed to generate a classification (through a tree) from distances between data. This distance matrix is noted d (d is the element in the iest row and the jest column, which corresponds to the distance between item i and item j; d is then a N × N matrix where N is the number of items). In the field of molecular evolution, data are DNA (sometimes RNA) or protein sequences, represented as strings of characters, and distances are calculated between aligned sequences. Least Square methods purpose is to generate a tree that preserves distances between data “as much as possible” in the sense that it should minimize a selected criterion (eg, Eq. 1).

The general aim of Multi Dimensional Scaling is also to provide a configuration of points that preserves “as much as possible” distances between data (and a criterion is also selected to that purpose). Multi Dimensional Scaling can be then considered as a function that associates N items in a metric space to the original distance matrix. The output space is frequently a Euclidian two-dimensional space (which is a well-known metric space). However, any metric space could be obviously used here instead. Least Square methods (such as the Fitch-Margoliash method) can be then reformulated as Multi Dimensional Scaling where the distance in the output space is an additive distance (ie, a distance calculated on a tree).

In fact, striking parallel features can be drawn between Multi Dimensional Scaling and Least Square methods: inputs (distances matrices), purposes (preservation of distances “as much as possible”) and outputs (configuration of items in an intuitive metric space) are similar. The only difference between Multi Dimensional Scaling methods (including Sammon’s mapping) and tree building Least square methods (including the Fitch-Margoliash method) can be found in the “nature” of the output (and consequently, in the optimization methods). Moreover, we will show thereafter that Sammon and Fitch-margoliash’s criteria are similar. Multi Dimensional Scaling methods might therefore provide a framework to better understanding criteria for tree building methods.

Note that we consider unrooted trees in the present paper. Moreover, branches with negative length, a common feature, are considered as 0-length branches,35,36 distances in trees will have all the three mathematical properties of distances (∀x,y,z, 1) d(x,y) = d(y,x); 2) d(x,y) = 0⇔x = y; 3) d(x,z)≤d’x,y) + d(y,z)).37 Note also that no hypothesis is made on the input distance; moreover all mathematical properties of distances are not required here (triangle inequality can possibly be transgressed). As a consequence, any dissimilarity could be considered (in tree reconstruction as well as in Multi Dimensional Scaling).34

Fitch-Margoliash criterion

Least Square methods are achieved in order to minimize difference between original distances and distances in the resulting tree. In particular, Fitch-Margoliash proposes a criterion (Eq. 1) to quantify the misfit between the original distance d and the distance matrix between items positioned in a given tree (we will designate this tree-distance matrix as ∂). where w = f(d). Several functions have been proposed for f (f is thereafter designed as the weighting function). Within the Least Square framework, w is often considered to be proportional to the reverse of an estimation of the variance of d, seen as an estimator of the “real” evolutionary distance between i and j.

Generally, f is chosen as f(x) =1/x, where p is positive. For example, w = f0(d) = 1 (ie, P = 0)38 corresponds to the Ordinary Least Square methods (OLS). It assumes that every distance is measured with the same error rate. Conversely, by choosing a strictly positive value for p leads to Weighted Least Square methods (WLS), we can naturally consider that the larger the evolutionary distance, the poorer its estimation and then the weaker its weight. In that case, w = f1(d) = 1/d could be considered appropriate (ie, P = 1). The most commonly used weight is (ie, P = 2).16 It is based on the observation that the variance of the error for the estimation of the evolutionary distance is often approximately proportional to in simulated datasets (these datasets are aligned sequences following classical evolutionary models).35 This estimation as been refined recently, to propose the following weight: (ie, P = 1.823).39

Similarity between Fitch-Margoliash and Sammon’s mapping criteria

The Sammon’s mapping28 has been designed to embed data in a low-dimensional space while preserving distances “as much as possible” with a special consideration for small distance. The method is known as one of the earliest non-linear projection proposed by John W. Sammon in 1969.

Striking parallel features can be drawn between Fitch-Margoliash and Sammon’s mapping criteria (for P = 1 in Eq. 1 that corresponds to the case where w = 1/d). Indeed, Eq. 1 exactly expresses criterion that drives the Sammon’s mapping if ∂ assigns distances in a Euclidean output space.

Multi Dimensional Scaling has been improved since it was first introduced (for an overview of the field, please report to Lee et al)40 and many non-linear Multi Dimensional Scalings have outclassed Sammon’s mapping. Progresses obtained in the field of non-linear mapping can be applied to Fitch-Margoliash method. In particular, we will focus on Data-Driven High Dimensional Scaling (DD-HDS), 34 an effective mapping method that follows the Sammon’s mapping concept.

Finding the best configuration

While w is chosen, (which lead to define the criterion ς, eq. 1), the best fitting tree has to be discovered. However, the number of possible tree topologies increases explosively with the number of items (number of possible topologies for N items = (2 N-5)!/[2 × (N-2)!] = (2 N-5) × (2 N-3) × … × 5 × 3).41 Of course, the exhaustive exploration is not an option, and finding the tree that best preserves distances is a NP-difficult problem.42,43 Numerous algorithms have been proposed to approach this goal including genetic algorithm,44–46 ant colony optimization47 and TABU search.48 In particular, the popular FITCH software from the PHYLIP package49 iteratively introduces items (each possible place is tested for the new item, the best place from the criterion point of view is validated). After the introduction of each item, a “tree swapping” optimization process is run to reconsider the global organization of the tree: the Nearest Neighbor Interchanges (NNI).50–52 Some alternative methods can be used in such framework rather than NNI. In particular, the “Tree Bisection-Reconnection” (TBR) and the “Subtree Pruning Regrafting” (SPR) can be used and allow a wider tree reorganization.51 However, such advantage is expensive in term of complexity: complexities are O(N) for NNI, O(N2) for TBR and O(N3) for SPR (N is the number of leaves).53 For that reason, and to be mimetic with the original Fitch-Margoliash algorithm, we used NNI in the present work.

A corresponding problem occurred in Multi Dimensional Scaling. However, even when methods agree about the criterion that should be optimized, various ways can be followed to generate the best configuration of points in the output space. Among optimization methods used in this framework, we can cite Generalized Newton-Raphson algorithm,28,54 TABU search,55,56 genetic algorithms,57,58 simulated annealing,59,60 neural networks29 and Force Directed Placement.34,60–62

Improvement of the Fitch-Margoliash Criterion

We have shown the proximity between Sammon’s mapping and Fitch-Margoliash criteria. However, some drawbacks have been identified for the Sammon’s mapping criterion. They obviously impact Fitch-Margoliash and have then to be avoided.

Penalization of “false neighborhoods” and “tears”

It has been found that Sammon’s mapping poorly penalizes “false neighborhoods”.34,63–65 Indeed, let suppose that d is high; f(d) is small; (d – ∂)2 × f(d) does not contribute much to ς, even if ∂ is small. A case where d is high and ∂ is small (which corresponds to a “false neighborhood”) is then poorly penalized.

This situation can be avoided using the Curvilinear Components Analysis criterion29 (Eq. 2) (the weight depends now on the distance in the representation space).

Unfortunately, “tears” are now poorly penalized: if d is small and ∂ is high (which corresponds to a “tear”), (d – ∂)2 × f(∂) remains small and does not contribute much to ς. Note that, consideration of “tears” and “false neighborhoods” are sometimes addressed under concepts of “continuity” and “trustworthiness”.65,66

However, “tears” and “false neighborhoods” can now be simultaneously penalized through a new criterion like the following:34 f1(min(d,∂)) is high if and only if d or ∂ is small: each small distance is considered as important (whatever if the distance is small in original or output space). This criterion proposed for Data-Driven High Dimensional Scaling has been found effective for avoidance of “tears” and “false neighborhoods” on several simulated and real datasets.

When considering trees, this analysis is still valid. In this context, two items that are connected whereas they are widely different is a “false neighborhood”, when close items that are placed in different parts of the tree corresponds to a “tear”. Naturally, within phylogenetic inference, both “false neighborhoods” and “tears” should be avoided. However, we demonstrated that Fitch-Margoliash (as well as Sammon’s mapping) does not penalize efficiently “false neighborhoods”. As a consequence, the consideration of criterion proposed in Eq. 3 could provide benefits in tree building like in Multi Dimensional Scaling.

Consideration on concentration of measure phenomenon (CMP)

The concentration of measure phenomenon

Dealing with high-dimensional data is complicated. Indeed, high-dimensional spaces have several surprising properties that discomfit human spirit (the famous “curse of dimensionality”).67,68 Among them, the concentration of measured phenomenon (illustrated in Fig. 1) is critical; it makes the relative difference between small and large distances tends to zero when the dimension increases.68–71 Such phenomenon can be illustrated by the distribution of distances between data points randomly drawn in a hypercube according to the dimension of the space.72

Figure 1

200 random data are uniformly distributed in a unit cube (of a given dimensionality). Histograms of distances between every pairs of items (19 500 distances) are displayed according to space dimension. Distributions of distances for dimensions larger than 200 would have the same Gaussian-like shape, but their centers would be shifted to the right proportionally to the square root of the dimension.

The concentration of measured phenomenon has a main impact on Multi Dimensional Scaling: even if small distances are believed to be emphasized in the mapping (through the weighting function f defined in section 2.1), the weights for small and large distances could be very close.15,34

As an illustration, let us considered two distances noted a and a +ɛ, where a is small and a + ɛ is large. When the dimension of space increases, ɛ/a may tend to 0 ( ). Then,

Weights for a and a + ɛ (noted f(a) and f(a + ɛ)) are then very close for high-dimensional data. Here, p (the parameter of the weighting function f; in Fitch-Margoliash, p often equals 2) is supposed to be a strictly positive integer for the calculus simplicity (this could be generalized to P ∈ [1, + ∞[using the gamma function).

Impact of CMP on phylogenic data

Proteins and DNA sequences are high dimensional objects (including when the only available material is the distance matrix between items).17 The question of their dimensionality has been theoretically addressed through the question of intrinsic dimension (or degree of freedom) of data.18,19 However, determining the dimensionality is not needed in the present framework. Still, the distance matrix may show properties close to ones of distances between high-dimensional data, ie, a low “relative contrast” (contrast can be defined as ([distance of the farthest to the origin of the vector space]—[distance of the nearest to the origin])/[ distance of the nearest to the origin]) in the distance matrix.69,73 In such case, the phenomenon described by equation 4 may occur: every distance has almost the same weight in contradiction with the spirit of Fitch-Margoliash method. In case of phylogenetic data, the origin of an eventual concentration of measured phenomenon could be discussed. However, such discussion can be avoided here, indeed we do not consider the dimensionality here, but the concentration of measure itself (whatever its origin).

The Multi Dimensional Scaling method DD-HDS34 has treated this problem: a sigmoid function adapted on the distances distribution allows us to assign a high weight to smaller distances and a low weight to larger ones whatever the dimension. We propose to apply the same solution to Weighted Least Square methods: the weighting function f is modified according to equation 5. where g(u,μ(λ),σ (λ)) is the probability density function of a Gaussian variable with mean μ(λ) and standard deviation σ (λ). μ(λ) and σ (λ) are suggested to be: and where λ is a user-defined positive parameter, usually equal to 0.1, that reflects the trade-off between preservation of local versus global distances. λ can be initially set at a high value (0.9 for example) and be decreased during the optimization, initially to take into account the global organization of data and progressively focus on local relationships.

Noticeably, the use of a Gaussian density function does not result from a hypothesis of Gaussian distribution of measured but only allows a continuous decrease of function f with high values for smaller distances and low values for larger distances. Equations 6 and 7 allow adjusting parameters μ and σ for equation 5; in case of an atypical distances distribution, the user can therefore adjust parameters.

Within the phylogenetic analysis framework, sequences are somehow related to each other and distances are usually expected to follow patterns. A similar argument applies to dimensionality reduction methods, regarding distances between data points (that are supposed to lie onto a low-dimensional manifold). However, we can observe that patterns are sometimes blotted out by the concentration of measure phenomenon when dealing with real data. Obviously, we cannot rely on this hypothesis in the context of tree reconstruction (of course, in the favorable case where patterns can be found, the goal will be reached more easily).

A new criterion

DD-HDS introduces a new criterion which is optimized (the absolute difference between distances allows to be consistent with an optimization process called Force Directed Placement that simulates a spring system, not used in the present tree building method). In order to develop a Weighted Least Square method (which is more common in tree building framework), we can use a corresponding criterion that rather considers square difference between distances (Eq. 8). For sure, close minimums should be found for these criteria. The optimized criterion is now: where λ is the user-defined parameter defined in the previous section. λ must be set between 0 (for maximally emphasizing of small distances) and 1 (where large distance still weight).

In Least Square methods, w is often considered as the reverse of an estimation of the variance of d, which is not the case for our criteria. Classically, w is selected based on this perspective from observations on simulated data. These simulated data follow a definite form, since the alignment of artificial sequences is strictly based on assumption regarding the evolutionary model. When the user’s data follow a similar model, such a criterion must be chosen. Conversely, when data do not follow the usual standard evolutionary models or when distances are not measured in the same way or when data are not aligned sequences etc, the classical Fitch-Margoliash criterion is not necessarily more suitable than any other one. In such cases, we gave important arguments to prefer using weights as presented in section 3, subsection “Consideration on concentration of measure phenomenon (CMP)”. As a consequence, eq. 8 provides a model-free criterion (in the sense that it is based on no explicit model but it is designed to adapt to the distance distribution).

Algorithm

We exposed a new criterion and gave theoretical arguments to prefer its use instead of the classical Fitch-Margoliash’s one. To test our assumptions, an algorithm using the new criterion is set up.

Optimization

Our algorithm mimics FITCH software49 (except, naturally, for the criterion). Three randomly selected data are initially chosen as leaves of a 3-branch tree (one single possible topology). Remaining leaves are then randomly selected to be introduced one by one according to an iterative procedure: each position is tested for the new leaf; when the best position is chosen, the whole topology is challenged by testing branch permutations (the so called “tree swapping” procedure) according to the NNI (Nearest Neighbor Interchange) algorithm. When all leaves have been introduced, the procedure is stopped. The only difference between algorithms lies in the criterion (c.f. eq. 1 for example).

A modification is nevertheless necessary for the calculus of branch lengths. Originally, branch lengths depend on original distances and could then be analytically deduced from a given topology and chosen weighting simply using matrix algebra.74 However, weights can no more be fixed a priori, because they now also depend on distances in output tree. A slight adaptation using an incremental approach is then necessary. Firstly, weights are fixed as w = f(d,λ)2 and branch lengths are calculated. We then get estimations for ∂, which allows the modification of weights as w = f(min(d,∂),λ)2 that generates new estimations for ∂ (each of these iterations is a linear operation). This process is iterated until stability of ∂ is reached. This modification does not strongly affect the calculus time. Indeed, few quasi-instantaneous steps (each one is a simple matrix algebra operation) are generally needed to converge.

The use of progressive reduction of λ (c.f. section 3.2), at worst multiplies the calculus time by the number of steps (here 2 steps). Overall these procedures increase the calculus time by less than a 2-fold factor (and, as a consequence, has no impact in term of complexity). It is worth to note that the overall complexity of the algorithm is O(N4), just as the classical Fitch-Margoliash (FITCH implementation).75

Expected benefits

If distances are additive there is a perfect solution to Least Square methods. Obviously, Fitch-Margoliash (just as most of classical methods like Neighbor Joining, etc) will solve this: no benefit can be expected in this ideal case. However, distances are often far from being additive, and it may be hard to reach the best tree in these cases.

Since additivity is a very restrictive property for distances, there are few chances to be in that case when dealing with real data. Fitch-Margoliash has shown efficiency in such a situation. However we highlight some drawback in this paper that may be a matter of concern. If, in addition, the distribution of distances is tight (which corresponds for example to distances in a high-dimensional space), we have shown that Fitch-Margoliash can be overtaken for theoretical reasons. As a consequence, better results can be expected in such case (and/or when the user’s data do not follow classical evolution models, cf. section 3, subsection “A new criteria”) because of a criterion that quantifies more appropriately the correspondence between the distances and the tree.

Moreover, during inference processes, the optimization to find the best tree is the major difficulty. We believe that the penalization of both “false neighborhoods” and “tears” will allow us to avoid the exploration of many poor solutions that were lightly penalized by classical criteria. As a consequence, the risk of reaching a local minimum is reduced: our method may therefore drive us more closely to the optimal tree.

Eventually, we propose the choice of a parameter λ that balances the matter of small versus large distances in order to allow a control by user.

Results

Methods and criteria

We evaluated trees within six criteria to test the benefits obtained by each one of our modifications (Fig. 2). The benefit supplied by every proposed modification can thus be sized up.

Figure 2

Summary of the eight tested criteria and links between them. Every combination of criterion components is tested to evaluate each improvement. Components allows: penalizing either “tears” (as in original criteria) or “tears” and “false neighborhoods” together (as it is suggested here); emphasizing small distances either thanks to the traditional reverse function or while considering the concentration of measure phenomenon with more (cases ζ2 and ζ5) or less (cases ζ1 and ζ4) focusing on smaller distances. ζFM corresponds to the original Fitch-Margoliash’s criterion and ζSa to the Sanjuán and Wróbel’s one. ζ1 and ζ2 avoid both risks of “false neighborhoods” and risks related to the concentration of measure phenomenon. ζ3 and ζ6 avoid risks related to the concentration of measure phenomenon and risks of “tears” (but do not considered the risk of “false neighborhoods”). ζ4 and ζ5 avoid risks related to the concentration of measure phenomenon (while highly focusing on smaller distances in case of ζ5) but lightly penalize “false neighborhoods”.

These criteria allow us to build trees. Ti corresponds to trees generated by minimizing ζi. TFM corresponds to the tree based on the original Fitch and Margoliash’s criterion ζFM and TSa corresponds to the tree based on the Sanjuán and Wróbel39 criterion ζSa. Several variations around the criterion that have been previously discussed are tested (criteria ζ1, ζ2, ζ3, ζ4, ζ5 and ζ6, described in Table 1). Note that T2 (respectively T5) is calculated from T1 (respectively T4) thanks to NNI transformations, within the process of progressive reduction of λ. TNJ corresponds to trees generated by Neighbor Joining,11 which is an agglomerative method for tree building from the distance matrix. Neighbor Joining (often noted NJ) is designed to generate the tree topology that gives the least total branch length. It is known to be a very fast method for phylogeny, which often provides a good preservation of distances even in maximum likehood methods, which are slower methods, were shown to be much more accurate than NJ. As a consequence, Neighbor Joining will be also compared to our method. We will then test here eight criteria and nine ways to create trees.

Table 1

Kappa coefficients (c.f. section 5.3) for each method (in column) and for various k (in row). Note that the whole number of distances between 15 species is 105.

	T1	T2	T3	T4	T5	TFM	T6	TSa	TNJ
1	0.85	0.98	0.91	0.85	0.98	0.9	0.88	0.89	0.86
2	0.83	0.96	0.9	0.84	0.95	0.89	0.90	0.89	0.87
3	0.8	0.96	0.87	0.81	0.95	0.87	0.86	0.86	0.88
4	0.81	0.96	0.88	0.82	0.94	0.88	0.87	0.87	0.85
5	0.82	0.94	0.87	0.83	0.92	0.87	0.86	0.86	0.86
10	0.85	0.88	0.87	0.85	0.87	0.86	0.86	0.86	0.87
15	0.85	0.85	0.86	0.85	0.85	0.86	0.86	0.86	0.87

Evaluation plots: “continuity” and “trustworthiness”

Within Multi Dimensional Scaling framework, the distance preservation is often analyzed thanks to the Shepard’s dy-dx diagram.29,76 However, in tree-building framework, an equivalent diagram is not as easily interpretable (data not shown). To explore the distance preservation in a more appropriate way, we split the distance preservation through two curves that display “continuity” and “trustworthiness” (concepts introduced by Venna and Kaski).65

On the one hand, “continuity” quantifies the preservation of small original distances. As a consequence, a “tear” (close data in original space that are displayed widely separated) is a violation of the continuity. On the other hand, “trustworthiness” corresponds to the ability of the tree to express small distances. As a consequence, a “false neighborhood” (far distanced data in original space that are displayed as neighbors) is a violation of the trustworthiness.

Trustworthiness and continuity are displayed here as two curves: 1) the curve that draws the square differences between original and output distances as a function of the original distance allows evaluating the continuity, and 2) the curve that draws the square differences between original and output distances as a function of the output distance allows evaluating the trustworthiness. Of course, these curves are often not slick; according to the user’s purpose, these curves can be smoothened or not (a smooth curves allows the visualization of the general behaviors while an unsmooth curve reveals default). In the present paper, most figures show smoothened curves of “continuity” and “trustworthiness”.

Some data randomly embedded (uniform distribution) in a two-dimensional Euclidian space were analyzed here (Fig. 3, upper insert). Trees are first generated from the matrix distance according to the Fitch-Margoliash method and then by the new method. Some differences can be noticed between the topology of the two resulting trees (lateral inserts). The comparison between original and resulting tree distances are presented in the lower insert in terms of “continuity” and “trustworthiness” (lower insert). Figure 3 illustrates the algorithm progress but does not allow conclusion about the characteristics of each criterion.

Figure 3

Example of comparison between two tree-building methods: Tree built according to the classical Fitch-Margoliash method versus the tree built thanks to the new criterion. Original data (and the associate distance matrix) are displayed in the upper insert. Left and right inserts express the two trees. “Continuity” and “trustworthiness” can then be compared on the lower insert (the grey curve corresponds to the Fitch-Margoliash method and the black curve is related to the new method).

Kappa coefficients: preservation of k smallest distances

In order to compare the preservation of short distances in each method, the concordance between the k smallest original distances and the k smallest distances in trees is analyzed by mean of Kappa coefficients.77 Kappa measures the degree of concordance on a scale from minus infinity to 1. A Kappa of one indicates full concordance, a Kappa of zero indicates that there is no more concordance than expected by chance and negative values are observed if concordance is weaker than expected by chance (a very rare situation). The measure could be seen as another mean to quantify the quality of resulting trees.

Kappa values are generated for each method and for various size of neighborhood (ie, various k).

Simulated data

Randomly positioned data

To test the capability of each criterion to generate fine trees from matrix of no additive distances, we created distance matrices from data randomly embedded in a Euclidean space. We tested sets of 15 data (ie, species) in two-dimensional and 100-dimensional spaces. Using 15 data points only could seem low, but it is clearly large enough to compare algorithms. Indeed, there is 7.9 × 1014 possible unrooted trees (which provide many chances for differences between results). Moreover, using 15 data points fits with the order of magnitude of real phylogenic problems. Lastly, 15-leaves trees are computable in a reasonable time; a large number of datasets can be simulated, which ensure to limit biases due to atypical situations.

Simulated sequences

Biological sequence evolution was simulated using a branching process. So, they are always related to a tree topology. For each time unit, all sequences present in the process can both be duplicated (with a fixed probability of duplication) and be mutated. The process starting from i) a sequence a0 ii) a probability distribution with the a0 amino acids position as support for this distribution reflecting the profile from which is derived the a0 protein family and iii) a duplication rate which can be allowed to vary with time or not. This program named EvolSeq was written in language C and can be requested to the authors.

In the present paper we mainly test the tree building methods with randomly positioned data; but several results reached with simulated sequences are also shown as supplementary material. Why was randomly positioned data preferred? Because randomly positioned data can maintain a perfect control on the topology and the dimensionality of the dataset. Conversely, as far as we know, whereas degrees of freedom are often controlled in case of simulated sequences, the sequences can however lie onto a subspace that dimensionality and topology are unknown. Indeed, let us consider biological sequences under evolution process as data that follow paths in a sequence space (the so called Configuration Space of Homologous Proteins, CSHP, after Bastien et al).18 The dimensionality of such a space is the number of varying sites, but sequences more often live in a subspace with a much lower dimension (due to the co-varying sites and so on). In such a case, the possible configuration of sequences allowed by the process (which corresponds to the topology of the subspace) is unknown. However, highlighting the critical impact of the subspace dimensionality is one of the main claims of the present paper.

Two-dimensional space

A distance matrix has been created from 15 points embedded in a two-dimensional Euclidean space (random points, uniform distribution). Nine trees have been produced according to the nine tested criteria (with the same order for the insertion of data). We then observe various topologies according to the considered criterion. The average of Robinson and Foulds distances (Robinson and Foulds)50 allows observing the proximity of tree topology obtained by each method (up and left insert in Figure 4 display these distances on a two-dimensional map thanks to the Data-Driven High Dimensional Scaling algorithm). For each method, “continuity” and “trustworthiness” are displayed as curves. To avoid considering the peculiarities of a given dataset, 200 sets of 15 data have been successively generated. For each dataset, trees have been calculated according to the various tested methods. This procedure led to 200 results that are merged in order to compare methods: the presented Robinson and Foulds distances corresponds to the average of the 200 Robinson and Foulds distances and trustworthiness and continuity curves are generated from every resulting couple of distances (that is 200 × 15 × 14/2 = 21000 distances).

Figure 4

Evaluation of distance preservation by nine tree building methods (analysis on 200 sets of 15 random data in a two-dimensional space). The up and left insert shows the Robinson and Foulds distance between trees generated by the various methods: the distance that separates two methods in the graph accounts for the average Robinson and Foulds distance between methods’ trees (see Hillis et al)78 for a somewhat similar display of Robinson and Foulds distances). Other inserts express “continuity” (lefts inserts) and “trustworthiness” (right inserts). Every curve is reported on each graph to allow an easy comparison. On each insert, the related method corresponds to the black curve.

Methods lead to different results, even from the topology point of view. Criteria that use the weighting function with high consideration of small distance (criteria ζ2 and ζ5 ) provide the best “continuity”. Note that they both lead to close results (c.f. Robinson and Foulds distance). However, ζ5 does not provide a high “trustworthiness”: many small distances in the tree correspond to large original distances which relies on “false neighborhoods” predicted by theoretical arguments (indeed ζ5 has been said to few penalize “false neighborhoods”). Robinson and Foulds distance also reveals that other criteria (ζ1 ζ3 ζ4 ζ6 ζ and ζ) lead to close results except for the Neighbor Joining (NJ). Such a simple dataset highlights the reality of the problem mentioned in section 3, and shows that some improvement can be done to circumvent these problems.

Kappa values are globally close to 1, which indicates that the distances sorting is generally well preserved (Table 1). Here again, criteria that focus on small distances (ζ2 and ζ5) clearly stand out. However, ζ2 (the proposed criterion) seems to lightly outclass ζ5. These results are clearly statistically significant. Indeed, after reproducing the estimation of Kappa values procedure 30 times, we observed a standard deviation around the values equal to 0.01 (a difference higher than 0.02–2 × σ–can be considered as significant).

High-dimensional space

The same procedure has been used on datasets randomly generated in a 100-dimensional space (the quantity of data and the number of iterations remain stable) and the result is illustrated in Figure 5.

Figure 5

Evaluation of distance preservation (data in a 100-dimensional space). The up and left insert shows the Robinson and Foulds distance between trees generated by the various methods. Other inserts express “continuity” (lefts inserts) and “trustworthiness” (right inserts). Every curve is reported on each graph to allow an easy comparison. On each insert, the related method corresponds to the black curve.

As previously, ζ2 and ζ5 show the best preservation of short original distance (“continuity”). ζ5 provides a poor “trustworthiness” due to the presence of “false neighborhoods” predicted by theoretical arguments (same in the previous test). However, in case of high-dimensional data, results from ζ2 and ζ5 are somewhat different (high Robinson and Foulds distance). Moreover, results obtained with ζ1 and ζ4 can be discriminated from those obtained with ζ3,ζ, ζ6,ζ and NJ. This difference is expected because ζ1 and ζ4 use the sigmoid weighting function: in case of high-dimensional data, only the sigmoid differentiates large and small distances (c.f. section 3).

Obviously, preserving the distances sorting is much harder for high dimensional data (values in Table 2 are weaker that the ones in Table 1). This situation can surely be related with the concentration of measured phenomenon that ensures that all distances are more or less similar (c.f. section 3).

Table 2

Kappa coefficients (c.f. section 5.3) for each method (in column) and for various k (in row). Note that the whole number of distances between 15 species is 105.

	T1	T2	T3	T4	T5	TFM	T6	TSa	TNJ
1	0.77	0.87	0.75	0.78	0.89	0.77	0.74	0.77	0.77
2	0.71	0.87	0.68	0.71	0.89	0.69	0.67	0.69	0.69
3	0.68	0.86	0.66	0.71	0.84	0.66	0.65	0.67	0.66
4	0.69	0.85	0.65	0.73	0.83	0.65	0.65	0.65	0.67
5	0.69	0.84	0.65	0.72	0.81	0.64	0.64	0.64	0.66
10	0.67	0.72	0.64	0.7	0.71	0.64	0.64	0.64	0.65
15	0.68	0.64	0.65	0.69	0.67	0.65	0.66	0.65	0.65

Contrast between Kappa values is much weaker here, which indicates that most methods are more or less equivalent from the distance sorting point of view. However, methods using ζ2 and ζ5 are clearly ahead of the curve in quality. The standard deviation of the Kappa values estimation equals to 0.015: a difference higher than 0.03 can then be considered as statistically significant.

Test for the robustness of optimization processes

The NNI procedure can be seen as a gradient descent algorithm used to optimize the tree topology (monitored by various criteria). It is well known that a risk of gradient methods is to fall into local minimums. To quantify this risk according to each criterion, we used a randomly embedded dataset (8 data; two-dimensional at first, and then in 100-dimensional spaces). 20 random trees are then used as initializations of a NNI process (same data and same initializations are used for every criterion). The tree with the minimal value for the criterion is considered as the global minimum: the percentage of resulting trees corresponding to this solution quantifies then the robustness of optimization using a given criterion. To be free from the feature of a particular dataset, this test is iterated 1000 times, and results are averaged (Table 3).

Table 3

Robustness of NNI according to the criterion that drives optimization. Each method is evaluated by the percentage of generated tree that corresponds to the optimal tree according to its criterion.

	ζ1	ζ2	ζ3	ζ4	ζ5	ζFM	ζ6	ζSa
dim = 2	40.9%	61.7%	44.8%	30.9%	56%	33.4%	39.1%	33.9%
dim = 100	45.2%	48.8%	45%	41.3%	48.6%	42.1%	45.2%	42.2%

Because NNI is run from a tree containing all items, we were able to treat here an unfavorable case: indeed, the progressive introduction of items in the tree should often help to reach the optimal minima. This can explain the somewhat low robustness for most methods (Table 3). However, obviously, an iterative introduction of items (as implemented by most methods) cannot guarantee that the optimal minima will be reached and the robustness is then an important property of tree building methods.

It can be seen that the use of the sigmoid weighting function with a low value for λ (that allows focusing on shortest distances) strongly increases the efficiency of the algorithm, as well as its robustness. Indeed, In Tables 1 and 2, we can observe the benefice in terms of efficiency though a gap between results reached from criterion 2 and the ones from criterion 3 as well as results reached from criterion 5 and the ones from the FM criterion (in both cases original Fitch-Margoliash weighting function is replaced by the sigmoid weighing function, please report to Fig. 2). Similarly in Table 3, the robustness of criterion 2 (respectively 5) is clearly higher than the robustness of criterion 3 (respectively, the FM criterion).

The penalization of both “false neighborhoods” and “tears” also clearly increases the robustness (by comparing criterion 3 to FM criterion or criterion 2 to criterion 5 in Table 2), but less strongly than the use of a sigmoid weighing function. However, it is not clear if it improves or damages the efficiency: trees from low-dimensional data seems to be very lightly improved (Table 1) but trees from high-dimensional data seems to be very lightly decreased (Table 2) when both defaults are similarly penalized. In order to evaluate the statistical significance of results presented in Table 3, we reproduced the estimation of robustness evaluation 30 times: the standard deviation around the values equal to 0.5% (a difference higher than 1% can thus be considered as significant).

Conclusion

Several properties derived from Multi Dimensional Scaling can be used within tree-building framework. In particular, we highlight the importance of taking “false neighborhoods” and “tears” into account while computing phylogeny. A new criterion is then proposed to overcome these risks and new evaluation processes using “continuity” and “trustworthiness” concepts allow us to check in it.

We also described the concentration of measured phenomenon and its possible influence when a tree is calculated according to classical methods. We then rectified the criterion in order to take account of the curse of dimensionality.

The analysis of simulated datasets highlighted that topology of trees may be different according to the chosen criterion. The new criterion outperforms other ones from distances preservation point of view. Moreover, a better robustness is observed.

As a consequence, the proposed criterion should be considered as a worthy alternative to the Fitch-Margoliash’s one. Future works should i) include simulations conducted in a more realistic way, for example by adding noise to the evolution process, ii) be evaluated on several real datasets and iii) be extended to maximum-likehood approaches.

A	B	C	D	E	F	G	H	I	J	K	L	M	N	O
0	17998	17998	17998	17998	17998	17998	17998	17998	17998	17998	17998	17798	17998	17998
17998	0	13858	13858	13858	13858	13858	13858	13858	13858	13858	13858	13858	16750	16750
17998	13858	0	10530	10800	10800	10800	10800	11344	11344	11344	11344	13678	16750	16750
17998	13858	10530	0	10800	10800	10800	10800	11344	11344	11344	11344	13678	16750	16750
17998	13858	10800	10800	0	2108	2472	2472	11344	11344	11344	11344	13678	16750	16750
17998	13858	10800	10800	2108	0	2472	2472	11344	11344	11344	11344	13678	16750	16750
17998	13858	10800	10800	2472	2472	0	804	11344	11344	11344	11344	13678	16750	16750
17998	13858	10800	10800	2472	2472	804	0	11344	11344	11344	11344	13678	16750	16750
17998	13858	11344	11344	11344	11344	11344	11344	0	4134	4134	8352	13678	16750	16750
17998	13858	11344	11344	11344	11344	11344	11344	4134	236	0	8352	13678	16750	16750
17998	13858	11344	11344	11344	11344	11344	11344	8352	8352	8352	0	13678	16750	16750
17798	13858	13678	13678	13678	13678	13678	13678	13678	13678	13678	13678	0	16750	16750
17798	16750	16750	16750	16750	16750	16750	16750	16750	16750	16750	16750	16750	0	722
17798	16750	16750	16750	16750	16750	16750	16750	16750	16750	16750	16750	16750	722	0

	A	B	C	D	E	F	G	H	I	J	K	L	M	N	O
A	0	3000	3000	3000	5296	5296	18406	18406	18406	18406	18406	18406	18406	18406	18406
B	3000	0	2088	2088	5296	5296	18406	18406	18406	18406	18406	18406	18406	18406	18406
C	3000	2088	0	444	5296	5296	18406	18406	18406	18406	18406	18406	18406	18406	18406
D	3000	2088	444	0	5296	5296	18406	18406	18406	18406	18406	18406	18406	18406	18406
E	5296	5296	5296	5296	0	5014	18406	18406	18406	18406	18406	18406	18406	18406	18406
F	5296	5296	5296	5296	5014	0	18406	18406	18406	18406	18406	18406	18406	18406	18406
G	18406	18406	18406	18406	18406	18406	0	2716	15744	15744	15744	15744	15744	15744	15744
H	18406	18406	18406	18406	18406	18406	2716	0	15744	15744	16248	16248	16248	16248	16248
I	18406	18406	18406	18406	18406	18406	15744	15744	0	7154	16248	16248	16248	16248	16248
J	18406	18406	18406	18406	18406	18406	15744	15744	7154	0	16248	16248	16248	16248	16248
K	18406	18406	18406	18406	18406	18406	16248	16248	16248	16248	0	3690	13248	13248	13248
L	18406	18406	18406	18406	18406	18406	16248	16248	16248	16248	3690	0	13248	13248	13248
M	18406	18406	18406	18406	18406	18406	16248	16248	16248	16248	13248	13248	0	3178	9164
N	18406	18406	18406	18406	18406	18406	16248	16248	16248	16248	13248	13248	3178	0	9164
O	18406	18406	18406	18406	18406	18406	16248	16248	16248	16248	13248	13248	9164	9164	0