Peter Richmond
email: peterrichmond@ymail.com

Bertrand M Roehner
email: roehner@lpthe.jussieu.fr

A joint explanation of infant and old age mortality

Infant deaths and old age deaths are very different. The former are mostly due to severe congenital malformations of one or a small number of specific organs. On the contrary, old age deaths are largely the outcome of a long process of deterioration which starts in the 20s and affects almost all organs. In terms of age-specific death rates, there is also a clear distinction: The infant death rate falls off with age, whereas the adult and old age death rate increases exponentially with age in conformity with Gompertz's law. Nevertheless it would be satisfactory to explain the two phenomena as two variants within the same explanatory framework. In other words a mechanism providing a combined explanation for the two forms of mortality would be welcome. This is the purpose of the present paper. We show here that the same biological effects can account for the two cases provided there is a difference in their severity: death triggered by isolated lethal anomalies in one case and widespread wear-out anomalies in the second. We show that quite generally this mechanism leads indeed, respectively, to a declining and an upgoing death rate. Moreover, this theoretical framework leads to the conjecture that the severity of the death effects, whether in infancy or old age, are higher for organisms which are comprised of a larger number of organs.

Introduction: infant versus old age mortality

In this paper we consider the shape of the curves of death rates1 as a function of age. Deaths in infancy versus old age death can be characterized in two ways [START_REF] Berrut | Deciphering infant mortality[END_REF][START_REF] Bois | Infant mortality across species. A global probe of congenital abnormalities[END_REF]):

Shape of the age-specific deathb rate

In infancy the death rate decreases with age whereas in old age it increases, see Fig. 1. Between birth and the age of 10 (note the log-log scale) the infant mortality rate falls off as a power law: µ b = A/x γ where the exponent γ is usually of the order of 1. After the infant phase comes the aging phase (note the lin-log scale) during which the death rate increases exponentially µ(x) = µ 0 exp(αx) in agreement with Gompertz's law. Source: Wonder-CDC data base for detailed mortality.

In medical terminology infancy refers to new born under one year of age. However, as the decrease continues until the age of 10, it seems appropriate to extend the meaning of the term to the whole age interval over which the death rate is decreasing. This is what was done in the two papers cited above and we use the same terminology here.

For humans, the increase of the death rate is described by the well-known law of [START_REF] Gompertz | On the nature of the function expressive of the law of human mortality, and on the mode of determining the value of life contengencies[END_REF]. This law can be summarized by saying that the death rate doubles approximately every 10 years of age.

Wear-in versus wear-out

In the terminology of reliability studies, infant mortality is described by a wear-in process, that is to say a phase during which the organs of the new-born start to work which results in the elimination of the organisms which are beset with an organ which does not work appropriately.

On the contrary, old age death is described as a wear-out process in which all organs experience damages due to continuous use. The lungs catch less oxygen, the bones become more fragile, the arteries become less elastic, and so on. Death occurs eventually due to the failure of a crucial organ but actually this failure is favored by the degradation of the whole organism. For instance, when the arteries become less elastic, when the lungs become less effective, it becomes more difficult for the heart to ensure blood circulation. This means that a heart failure does not come about in isolation but rather in relation with the wear of other organs.

The purpose of this article is to show how the feature (1) results from the feature (2). However, before coming to that we wish to explain how our study fits into the broad framework of aging and senescence studies.

Why do we think that infant mortality is an essential component in the understanding of aging? There are several reasons that are developed below.

2

The key-role of congenital malformations

The real challenges of aging models

Not surprisingly, the modeling of aging and senescence has received great attention. A comprehensive review can be found in two papers by Leonid Gavrilov andandNatalia Gavrilova (2001, 2006).

Readers will find in these papers a comprehensive and very readable account of the literature of aging models2 . In the following subsections we explain why infant mortality is a simpler and more fundamental effect than old age mortality.

Great diversity in the shape of the death rate in old age

Across species the hyperbolic decrease of the death rate in the infancy phase appears to be a phenomenon that is more widespread than its exponential (i.e. Gompertz-like) increase in adulthood. Indeed, there is much more diversity in old age death rates than in infancy death rates; see in [START_REF] Berrut | Deciphering infant mortality[END_REF] the graph based on zoo species.

In addition for some spłcies documented in Gavrilov et al. (2006, p. 18 and 33) there is a marked effect of leveling-off in old age. For instance, house flies have a maximum life span of 40 days but around the age of 15 days the exponential growth of the death rate is replaced by a section where it is practically flat.

In humans the dominant diseases in old age are not the same nowadays as one century ago. Presently, there is a prevalence of heart, cancer and Alzheimer's disease whereas around 1900 infectious diseases were still common. Thus, with organisms being confronted to different challenges, one should not be surprised to see changes in the shape of the death rate in old age.

Finally, the individuals who reach old age were "filtered" and selected by the diseases to which they were confronted. If one could observe the signature of the immune system one would see that the immune system of persons of old age is not the same in 2020 than in 1900, and also not the same in developing countries than in developed countries.

Common characteristic of embryonic and infancy death rates

It has been shown recently [START_REF] Chen | Magnitude and significance of the peak of early embryonic mortality[END_REF]) that for zebrafish3 the embryonic death rate is by far highest at the beginning of the embryogenesis, an observation which suggests that most of these deaths are due to mistakes in the manufacturing processes of the oocyte (femelle egg) and sperm cell.

If instead the deaths would be due to mutations during the embryogenesis they would be uniformly distributed or even (through a cumulative effect) concentrated in the late phase of embryogenesis.

Similarly, the infant death rate is by far highest immediately after birth. By the same argument, it appears that most of these deaths are due to faults in the manufacuring of the embryo. For instance, in mamals lung malformations are without consequence as long as the fetus receives its blood from the mother but they will lead to death as soon as this connection is interrupted. Embryonic and infant deaths along with the malformations which are not immediately lethal give us global information about the underlying manufacturing processes.

The effects of congenital defects and of aging occur jointly

At first sight it may seem that the infant death rate can be easily described and explained through the process of elimination of individuals with malformations. Clinical data show that in the first weeks after birth most of the deaths are due to congenital anomalies (percentage data are given in [START_REF] Berrut | Deciphering infant mortality[END_REF]. When the most serious malformations have been eliminated the rest of the cohort is less likely to die.

However, the previous explanation is not really satisfactory for the following reason. In fact, deaths due to congenital anomalies are not limited to young age but continue during the whole life. For instance, a congenital defect of heart valves may be of no consequence until the age of 60 or 70 when the defect becomes more serious because the valve's leaflets become stiffer (see [START_REF] Bois | Physical models of infant mortality: implications for defects in biological systems[END_REF].

In other words, the wear-in and wear-out processes should not be seen as occurring successively but rather simultaneously; it is their strength, not their existence, which changes in the course of time.

Immediately after birth, wear-in is completely dominant, whereas in old age it is wear-out which is predominant. I short, taken alone the elimination of congenital malformations cannot explain the decrease of the death rate. In order to make it work we need to define both wear-in and wear-out more precisely.

In the next section we will use the feature already mentioned above, namely that the infant mortality is usually due to a congenital defect in one important organ (e.g. heart, lung, brain, liver, and so on) whereas the wear-out is due to parallel degradation of various important organs.

Modeling the wear-in and wear-out processes

The first step is to decompose any organism into its vital organs. For instance, Fig. 2 shows a decomposition into 4 organs, that could be heart, lung, brain and temperature regulation.

Secondly, we must find a way to describe mathematically whether each organ (as well as the whole organism) is alive. We do this by defining for each part a random variable X i which is its age at failure. In the case of humans we could make this description fairly realistic by giving to the X i The second line shows a wear-out death as the consequence of uniform deterioration of all vital organs. The graphs on the right-hand side show the implications of these mechanisms in terms of age-specific death rates: decreasing for wear-in as observed in infant death, increasing for wear-out as seen in old-age death.

values from interval (0, 125) for it should be remembered that 125 years is an upper bound of human life4 . However, as we do not wish to restrict ourselves to only the human species, we will normalize the interval of the X i to (0, 1) where 1 represents the maximum life span of the species. Now comes the most important step which is to describe the wear-in and wear-out mechanisms. Let us begin with the simplest case which is the wear-out mechanism.

Wear-out

The fact that the death of an individual occurs when the last surviving organ fails is expressed by saying that if X 1 = 0.5, X 2 = 0.3, X 3 = 0.7, X 4 = 0.1 (the X i denote the age at death of vital organs as shown in Fig. 2), then the age of death represented by the random variable Z will be Z = 0.7, in other words:

Z = Max(X 1 , X 2 , X 3 , X 4 )
For the sake of simplicity we assume that the X i are independent and identically distributed random variables. This assumption has the merit of making the analytical derivation possible. However, in specific applications one can take realistic distributions based on clinical data.

• If f (x) and F (x) respectively represent the density function and the cumulative distribution function of the X i , what will be the density function, f Z (x), of Z? Let us first consider the case of only two organs.

F Z (x) = P {Z ≤ x} = P {Max(X 1 , X 2 ) ≤ x} = P {X 1 ≤ x and X 2 ≤ x}
Now, the fact that X 1 and X 2 are independent means that:

P {X 1 ∈ A and X 2 ∈ B} = P {X 1 ∈ A}P {X 2 ∈ B}
where A and B are two subsets of the set of real numbers. Therefore:

F Z (x) = P {X 1 ≤ x}P {X2 ≤ x} = [F (x)] 2
which, by differentiation leads to the density function of Z:

f Z (x) = 2F ′ (x)F (x) = 2f (x)F (x)
For p parts instead of only two, one gets similarly:

F Z (x) = [F (x)] p ⇒ f Z (x) = pf (x) [F (x)] p-1
(1)

In order to see what is the shape of this function we consider the simple case of a random variable with a uniform density over the interval (0, 1); in this case:

for x ∈ (0, 1) : f (x) = 1, F (x) = x Thus, for x ∈ (0, 1) : f Z (x) = px p-1
This function is shown in Fig. 3b for p = 2, 4, 8, 15 We see that it is a fast increasing function of age.

According to the analytical expression f Z (x) is a power law function but when p becomes large it has the shape of an exponential (as shown in Fig. 2b for p = 15), a result which is qualitatively consistent with Gompertz's law according to which the probability of death increases exponentially with age.

Wear-in

For the example considered above wear-in death would mean that the age of death is: W = 0.1, i.e.:

W = Min(X 1 , X 2 , X 3 , X 4 )
For the distribution fonction F W (x) we can write:

F W (x) = P {W ≤ x} = 1-P {W > x} = 1-P {Min(X 1 , X 2 ) > x} = 1-P {X 1 > x and X 2 > x}
Again using the independence property, one gets:

F W (x) = 1 -P {X 1 > x}P {X 2 > x} = 1 -[1 -F (x)] 2
Then, as above, this result generalizes to:

F W (x) = 1 -[1 -F (x)] p ⇒ f W (x) = pf (x) [1 -F (x)] p-1 (2)
For the case of uniform random variables, one gets:

for x ∈ (0, 1) : f W (x) = p(1 -x) p-1
which means that the probability of death is a decreasing function of age, consistent with what is expected for infant mortality. The decrease is ilustrated in Fig. 3a for p = 2, 4, 8, 15.

Special cases

One can gain an intuitive understanding of the theoretical framework by considering a number of special cases.

Firstly, we can consider the extreme case of an organism with a high number of components, say one million or if you prefer a number p which tends to infinity. Then, intuitively, the wear-in assumption gives a probability of death equal to 1 because for such a large number of components there will always be one which will fails almost immediately.

For the same reason the wear-out assumption gives a probability of death equal to zero for almost all ages because it will take a very long time to eliminate all and any components.

The other extreme case is an organism with only one vital component. Then, obviously, the two assumptions should give the same result. Indeed, the formulas ( 1) and ( 2) give:

f W (x) = f Z (x) = f (x).
Here the shape of f W and f Z is completely determined by f (x) which can have any shape, whether increasing or decreasing.

As the number p of components increases the factor (1 -F ) p , which is a decreasing function will become more and more predominant. Similarly, for old age, the factor F p , which is an increasing function, will become predominant when p increases.

In summary, this discussion makes clear that equations ( 1) and ( 2) do not describe only one model but, by playing with p and f (x), they can describe a whole spectrum of cases. It is in this sense that the model is really predictive.

Graphs of infancy and old age death rates

The density functions of the variables W and Z give the infant and old age death rates respectively (Fig. 3a,b). The graphs suggest the conjecture that the larger the number of components, the steeper the death rates with respect to age. When p → ∞ it can be seen directly on the formulas that f W falls vertically from pf (0) to zero, whereas f Z jumps vertically from zero to pf (1).

As technical systems have often many well-defined, separate components it should be possible to test our conjecture. For instance microprocessors comprising a large number of chips should have steeper death rate curves than those with only a small number.of chips. Unfortunately, for such technical systems there are almost no life-time data publicly available, probably for reasons of commercial confidentiality.

Conclusion

We have proposed two paradigms of death: (i) single-organ death which occurs through the failure (e.g. due to congenital malformation) of one crucial organ, (ii) multi-organ death which comes about through the deterioration of almost all organs.

At this point a distinction should be made between the underlying cause of death (e.g. cancer or an infectious disease) and the immediate cause of death. It is the former which is of interest and is reported in the death statistics by cause of death. On the contrary, the immediate cause of death is almost always the same, namely a heart failure. For instance, a liver cancer will lead to blood poisening which makes the heart unable to perform its function properly. It is because other organs (e.g. the kidneys) are also in poor shape that the partial defect of the liver eventually proves fatal. It is in this sense that death occurs as as a kind of overall collapse.

The model predicts that the slope of the death rate is higher (whether for fall or increase) when the number of organs is larger. Can this conjecture be tested on biological systems? The answer is "yes" and "no". "Yes", because a comparison of various species shows that there are indeed great differences in the number of organs. For instance, rotifers (a small swimming animal about 200 micrometers in length) have no heart, no blood, no lungs, no kidneys. In short, they have much less organs than fish. "No" because there is a serious obstacle, namely the great difference in lifetimes. Whereas rotifers live about 5 days, zebrafish live about 5 years. Unless one knows how to normalize the respective times there can be no meaningfull comparison for, needless to say, age normalization affects the measurement of the slope.

A Appendix A. Statistical and probabilistic descriptions of the death process

In this appendix we discuss some particular aspects of the theoretical framework. The first subsection clarifies the connection between age and time.in aging processes. The second subsection establishes the important connection between the density function of the ages of death and its frequency counterpart commonly refered to as the death rate.

A.1 Age versus time

Should we use age (noted x) or time (noted t)? At first sight, the question may seem irrelevant for if time is counted from the moment of birth the two variables are identical. It is not so simple however, as shown by the fact that in our min-max argument we had to introduce as many age of death variables X 1 , X 2 , . . . as there are vital organs. Moreover, for any organ, its real age is the time elapsed since it was created. For instance, in zebrafish the heart appears and starts to beat some 20 hours after fertilization. Most vital organs are created during embryogenesis which means that their age is not identical with time (measured after birth which for fish means hatching of the eggs). That is why X was a more appropriate.variable than time.

However, in this appendix our perspective is different for we wish to consider the evolution of a whole population or more precisely of a cohort of individuals born at the same moment. Taking this moment as origin of the time axis makes the age of each individual numerically identical with the time given by an external clock.

Actually, to describe the evolution of a cohort, time seems a better variable than its age. Why? Although the population exists at any time, for external observers, it becomes real only when we can know its size and that occurs only when a census (or a survey) takes place. Censuses are conducted at specific time intervals (e.g. in the US every decade) and concern simultaneously all cohorts. That is why in this part calendar time seems to be the natural variable. Age will play a role only if we wish to consider different age groups.

A.2 Death rate versus probability density

As always, the tricky point is the relationship between the probabilistic notions and their statistical counterpart. The goal of this appendix is to recall the main notions and how they are related.

Let t denote the age of individuals.in a cohort Let y(t) denote the size of the cohort at time t We wish to describe the decrease of the population in the course of time.

The probability that an individual would die in the time interval (t 1 , t 2 ) (which is also an age interval) is:

(Number of those who die)/(number of those alive initially) = [y(t 1 ) -y(t 2 )]/y(t 1 )

to get the probability of dying per unit of time we must divide by the length of the time interval ∆t = t 2 -t 1 .

µ(t) = (1/∆t) [(y(t 1 ) -y(t 2 )) /y(t)]

Note that µ(t) represents what is usually called the death rate, sometimes also called the hazard rate or the force of death. Note also that: y(t 1 ) -y(t 2 ) = -∆y.

Now, let us consider the case of a constant probability of dying. One is led to: (1/y)(dy/dt) = -a which gives: y(t) = y 0 exp(-at). In other words, the survival function is a decreasing exponential. Now let us consider a random variable T which represents the age of death of an individual. Its density function, defined by: f (t)dt = P {t < T < t + dt}, is the derivative of the distribution function: F (t) = P {T ≤ t}. f (t)dt is the probability that the death of the individual occurs in the age interval (t, t + dt) ; f (t) is the probability per unit of time.

Note that: f (t) ∼ µ(t). In words, µ(t) is the statistical counterpart of the probability density function of T .

If we consider again the case of a constant probability of dying (for ages in a bounded interval and zero elsewhere), it means: f (t) = µ(t) = a. Then, the distribution function is: F (t) = at, at least

Fig. 1 :

 1 Fig.1: The two phases of human mortality. The data are for the United States over the period 1999-2016.Between birth and the age of 10 (note the log-log scale) the infant mortality rate falls off as a power law: µ b = A/x γ where the exponent γ is usually of the order of 1. After the infant phase comes the aging phase (note the lin-log scale) during which the death rate increases exponentially µ(x) = µ 0 exp(αx) in agreement with Gompertz's law. Source: Wonder-CDC data base for detailed mortality.

Fig. 2

 2 Fig.2 Decomposition of an organism into vital organs and difference between wear-in and wear-out mechanisms. The first line shows a wear-in death. It is the consequence of the failure of a single vital organ.The second line shows a wear-out death as the consequence of uniform deterioration of all vital organs. The graphs on the right-hand side show the implications of these mechanisms in terms of age-specific death rates: decreasing for wear-in as observed in infant death, increasing for wear-out as seen in old-age death.

  Fig.3a,b Graph of the density functions of the Min and Max of a set of random variables.which represent the age-specific death rate in young and old age respectively. When the number of organs becomes large, the death rates become exponential.

We use the standard definitions of death rates, namely: µ(x) = ∆y/(∆x × y) where ∆y is the number of deaths in a given age interval of size ∆x and y is the size of the population at the beginning of the age interval under consideration. With this definition µ(x) is the probability (per unit of time) that a person who has reached age x, will die in the subsequent age interval, see appendix A.

Whereas these authors share our approach based on reliability science, and whereas infant mortality is a standard notion in reliability, they devote only one page (in a total of 58 for the two papers) to the question of infant mortality. This disproportion reflects the overwhelming predominance of aging and senescence studies.

For a study of embryonic death rates, zebrafish have two great advantages. (i) As for most species of fish, fertilization of the eggs occurs outside of the body of the female (ii) The shell of the eggs is transparent. Taken together, these two features imply that one can observe the embryos immediately after fertilization, something that is impossible either for humans, birds or rotifers.

This was shown to be a consequence of Gompertz's law in[START_REF] Richmond | Predictive implications of Gompertz's law[END_REF].

until at is equal to 1. The decreasing distribution function is:

Note that G(t) is different from the survival function.


Le Si Quang

Olivier Gascuel

Nicolas Lartillot
email: nicolas.lartillot@lirmm.fr

Empirical profile mixture models for phylogenetic reconstruction

Motivation: Previous studies have shown that accounting for sitespecific amino acid replacement patterns using mixtures of stationary probability profiles offers a promising approach for improving the robustness of phylogenetic reconstructions in the presence of saturation. However, such profile mixture models were introduced only in a Bayesian context, and are not yet available in a Maximum Likelihood framework. In addition, these mixture models only perform well on large alignments, from which they can reliably learn the shapes of profiles, and their associated weights. Results: In this work, we introduce an expectation-maximization algorithm for estimating amino-acid profile mixtures from alignment databases. We apply it, learning on the HSSP database, and observe that a set of 20 profiles is enough to provide a better statistical fit than currently available empirical matrices (WAG, JTT), in particular on saturated data.

INTRODUCTION

Capturing the evolutionary properties of the amino-acid replacement process in protein sequences has traditionally been done using empirical matrices. Such matrices are meant to account for the biochemically conservative pattern of amino-acid replacement in natural protein sequences. In practice, they are empirically derived from databases of pairwise or multiple alignments, either using counting methods [START_REF] Dayhoff | A model of evolutionary change in proteins[END_REF][START_REF] Jones | The rapid generation of mutation data matrices from protein sequences[END_REF], or by direct likelihood maximization [START_REF] Adachi | Model of amino acid substitution in proteins encoded by mitochondrial DNA[END_REF][START_REF] Adachi | Plastid genome phylogeny and a model of amino acid substitution for proteins encoded by chloroplast DNA[END_REF]Whelan and Goldman, 2001;[START_REF] Le | An improved general amino-acid replacement matrix[END_REF].

Empirical matrices appear to lead to simple and fairly accurate phylogenetic models (Whelan et al., 2001). On the other hand, they do not explicitely encode the fact that biochemical constraints are essentially site-specific features, related to context-dependent purifying selection. More precisely, depending on their position and role in the protein's overall shape and structure, sites will in general * to whom correspondence should be addressed accept only a very specific subset of the 20 amino-acids, all other possibilities being strongly selected against. For instance, buried sites will preferentially accept hydrophobic amino-acids, whereas residues at an active site are more likely to be electrostatically charged. This suggests, as an alternative way of describing protein evolution, probabilistic models explicitely formulated in terms of variations of the amino acid propensities across sites [START_REF] Bruno | Modeling residue usage in aligned protein sequences via maximum likelihood[END_REF][START_REF] Halpern | Evolutionary distances for protein-coding sequences: modeling site-specific residue frequencies[END_REF][START_REF] Crooks | An alternative model of amino-acid replacement[END_REF]. Such models need to be devised with caution, however, so as to control the overall parameterization. In this respect, two essential simplifications can be proposed.

First, assuming that biochemical constraints can be explained through equilibrium frequencies, one can rely on simple processes, such as F81 [START_REF] Felsenstein | Evolutionary trees from DNA sequences: a maximum likelihood approach[END_REF], that are entirely characterized by their profile of stationary probabilities over the 20 amino acids; the replacement rate between two amino acids does not depend on the initial state, and is simply proportional to the stationary probability (or equilibrium frequency) of the target state. Second, to avoid relying on site-specific stationary probability profiles, which may create an over-parameterization problem, we can assume that the variation of amino acid constraints is explainable as a finite mixture.

Mixture models have been used several times in phylogenetics, to model variations across sites of various aspects of the substitution or replacement process, such as the dN/dS ratio [START_REF] Yang | Codon-substitution models for heterogeneous selection pressure at amino-acid sites[END_REF], or even the complete substitution matrix [START_REF] Koshi | Models of natural mutations including site heterogeneity[END_REF][START_REF] Pagel | A phylogenetic mixture model for detecting patternheterogeneity in gene sequence or character-state data[END_REF]. More sophisticated mathematical devices, such as Hidden Markov Models, of which mixture models can be seen as a degenerate case, have also been proposed [START_REF] Felsenstein | A Hidden Markov Model approach to variation among sites in rate of evolution[END_REF][START_REF] Goldman | Using evolutionary trees in protein secondary structure prediction and other comparative sequence analyses[END_REF][START_REF] Thorne | Combining protein evolution and secondary structure[END_REF][START_REF] Goldman | Assessing the impact of secondary structure and solvent accessibility on protein evolution[END_REF]. In another direction, Markov Modulated Models have been introduced, as a way of allowing sites to switch their substitution behavior at any time along the lineages [START_REF] Holmes | An expectation maximization algorithm for training hidden substitution models[END_REF][START_REF] Gascuel | Modelling the variability of evolutionary processes[END_REF]. Mixture models can be seen as a degenerate case of Markov Modulated Models, in which the switching rate is equal to 0.

A mixture of profiles has been introduced in a previous work [START_REF] Lartillot | A Bayesian mixture model for across-site heterogeneities in the amino-acid replacement process[END_REF]. It is based on a Dirichlet process mixture [START_REF] Ferguson | A bayesian analysis of some nonparametric problems[END_REF][START_REF] Antoniak | Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems[END_REF], which is a flexible non-parametric device often used in Bayesian studies [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF][START_REF] Huelsenbeck | A dirichlet process model for detecting positive selection in protein-coding DNA sequences[END_REF][START_REF] Huelsenbeck | A nonparametric method for accomodating and testing across-site rate variation[END_REF]. The resulting model was shown to provide a better fit than empirical matrices in several cases [START_REF] Lartillot | Computing Bayes factors using thermodynamic integration[END_REF]. It also appears to better accomodate saturation (i.e. multiple substitutions), making it more robust to phylogenetic artifacts [START_REF] Lartillot | Suppressing long branch attraction artefacts in the animal phylogeny using a site-heterogeneous model[END_REF].

Non-parametric modelling based on Dirichlet process mixtures is a powerful method. One of its essential advantages is its flexibility, adapting to each particular dataset. However, a large amount of data is necessary for this non-parametric device to converge to a reasonably faithful description of the true distribution of profiles across sites. For this reason, it performs well on large alignments (more than 1,000 aligned positions), such as multiple gene concatenations. In contrast, it appears to be less efficient on smaller datasets, in particular, on single gene alignments. Another limiting feature is that Dirichlet process mixtures can only be practically implemented in a Monte Carlo framework, which makes them a method of choice for Bayesian studies, but a much less attractive tool in a Maximum Likelihood (ML) paradigm. Yet, the availability of a profile mixture model in a ML context could be useful in many situations.

Thus, we still lack an efficient, versatile (ML and Bayes) version of the profile mixture model, whose parameters would be pre-determined on empirical grounds, so that it could be used on small single-gene alignments as well as on large phylogenomic datasets. In the present paper, we develop such a model. Starting from a variant of the Expectation-Maximization method [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] especially devised for phylogenetic models [START_REF] Holmes | An expectation maximization algorithm for training hidden substitution models[END_REF], we specialize it to profile mixture models, and apply it on a subset of HSSP database [START_REF] Sander | Database of homology-derived protein structures and the structural meaning of sequence alignment[END_REF], to estimate a series of profile mixture models with an increasing number of components. We find that a good compromise is obtained in the form of a mixture of 20 profiles. When tested on phylogenetic alignments, this model gives, on average, a better fit than classical replacement matrices, such as WAG and JTT. Interestingly, profile mixtures are particularly more fit on saturated data, confirming that matrices specifically fail at correctly describing multiple conservative substitutions [START_REF] Lartillot | Suppressing long branch attraction artefacts in the animal phylogeny using a site-heterogeneous model[END_REF]. On the other hand, matrices are often better than profile mixtures when a low level of saturation is detected in the data. Finally, we show that 20 profiles are enough to reproduce some of the essential phylogenetic robustness properties previously observed using the non-parametric Bayesian version.

METHODS

Notations

The training database consists of a series of alignments (Dr) r=1..R . For simplicity, we explain the overall method for a single alignment D. The generalization to several alignments is straightforward. The alignment D is made of P sequences, of length N , which are assumed to be related by an unknown, rooted phylogenetic tree τ . In all cases investigated in this work, the topology is considered fixed. The sequences are taken from an alphabet of size S. In practice, we will work on amino-acid alignments (S = 20), although the method introduced here is also valid for nucleic acid data.

Let i index the columns C i , or sites, and j the nodes (0 ≤ j ≤ 2P -3), with root node having index 0 and leaf node indexed in [1, P ]. Branches are also indexed by j, 1 ≤ j ≤ 2P -3, with the convention that a branch has the same index as the node which is at its tip. Branch lengths are denoted by l = (l j ) 1≤j≤2P -3 , the relative rate of substitution at each site by r = (r i ) 1≤i≤N .

Each cell of the data matrix is referred to as x ij for i and j running over sites (columns) and taxa (rows), respectively. Thus, x ij is simply the state of the process operating at site i, at the leaf indexed by j (1 ≤ j ≤ P ). In addition, one will need to specify the state of the substitution process at all interior nodes as well. Let us simply expand the range of j so that now, x ij is defined for 0 ≤ j ≤ 2P -3.

Profile mixture models

Substitutions occur independently at each site, according to Poisson (or F81, Felsenstein, 1981) Markov processes running along the branches of the tree. One such Markov process is entirely characterized by a vector of stationary probabilities, or equilibrium frequencies π = (π(a)) a=1..S , such that P a π(a) = 1. In the following, we will call π the profile of the process.

The finite-time probability transition for a Poisson process takes a simple form: the probability of observing state b, after time t = l, and given that the process was in state a at t = 0 is

p(b | a, l) = e -l δ ab + (1 -e -l )π b , (1)
where δ ab is the Kronecker symbol (i.e. is equal to 1 iff a = b).

We will consider a mixture of K Poisson processes. π = (π k ) k=1..K , where each π k is a vector of frequencies over the S states of the alphabet. To each component of the mixture is associated a weight w k , such that P k w k = 1. We denote by w = (w k ) k=1..K the weight vector.

In addition, we assume variations of the overall rate of substitutions across sites. As in most standard phylogenetic models, we do this by introducing a discretized gamma distribution [START_REF] Yang | Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods[END_REF], of mean 1 and (inverse variance) parameter α. In effect, this amounts to introducing a series of Q rates r(α) = (rq (α)) q=1..Q , each of which is the mean rate within the qth. bin of the discretization. In general, bins are chosen so that all the weights are equal to 1/Q. The values of the discretized rates depend on the α parameter. We use Q = 4 throughout the present article.

The full parameter vector is thus Θ = (l, π, w, α). The likelihood at site i is a weighted average over all K × Q combinations of rates and profiles:

p(C i | Θ) = p(C i | l, π, w, α) = 1 Q X k w k X q p(C i | l, π k , rq(α)),
(2) and the total likelihood is the product over all sites:

p(D | Θ) = Y i p(C i | Θ).
(3)

Expectation-Maximization algorithm

An Expectation-Maximization (EM) method [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] for estimating amino-acid replacement models has been introduced previously [START_REF] Holmes | An expectation maximization algorithm for training hidden substitution models[END_REF]. It assumes a general (possibly Markovmodulated) time-reversible process, applied uniformly across sites. In the present context, we adapt the method so as to deal with mixture models, and to account for rate variations across sites. Concomitantly, we specialize it to profile models. The EM algorithm is particularly convenient for Markov processes, compared to other numerical likelihood maximization methods, essentially because the conditional expectations involved can often be analytically maximized with respect to the parameters of the model [START_REF] Holmes | An expectation maximization algorithm for training hidden substitution models[END_REF]. In the present case, the use of profiles leads to further simplifications in the analytical part of the computations.

Specifically, our EM algorithm relies on a combination of data augmentation and parameter expansion. The parameter expansion part consists in specifying the allocation of each site to one of the K components of the mixture, and to one of the Q rates of the discretized gamma distribution.

To do this, we introduce two allocation vectors y = (y i ) i=1..N , and z = (z i ) i=1..N , such that y i ∈ [1..K] and z i ∈ [1..Q] are two integer indicators specifying the allocation status of site i. The likelihood can be rewritten as a sum over all possible joint realizations of these two allocation vectors:

p(D | Θ) = X y,z p(D | y, z, Θ)p(y, z | Θ) (4) where p(D | y, z, Θ) = Y i p(C i | l, πz i , ry i (α)) (5) and p(y, z | Θ) = p(y)p(z | w) = 1 Q N Y i wz i . (6)
The data augmentation part consists in specifying the state of the process at all internal nodes of the tree, and for each site (x = (x ij )), as well as the number of substitutions along each branch and for each site, which we denote by n = (n ij ), n ij being the number of substitutions on branch j and for site i. We call Ξ = (x, n) a substitution mapping. Note that a more complete account of the substitution history would be possible [START_REF] Nielsen | Mapping mutations on phylogenies[END_REF][START_REF] Holmes | An expectation maximization algorithm for training hidden substitution models[END_REF], by including the states of the process after each of the n ij substitutions, and the exact position of each successive substitution event along the branch. However, we do not need it in the present case. Conditional on a particular allocation (y, z), the likelihood is a sum over all possible mappings compatible with the data at the leaves D:

p(D | y, z, Θ) = X Ξ|D p(Ξ | y, z, Θ). (7)
The function to be maximized in our EM is then the expectation of ln p(Ξ | y, z, Θ) over all possible realizations of mappings Ξ and allocations (y, z), conditional on the current value Θ * of the parameter:

Q(Θ, Θ * ) = E Ξ,y,z [ln p(Ξ | y, z, Θ) | Θ * , D] (8)
This expectation takes on a very simple form (see Appendix for details), that makes maximization with respect to Θ straightforward.

In practice, the optimization proceeds by alternating between four modules, performing an expectation followed by a maximization with respect to either l, α, π or w. In addition, as the alignments of the training database are independent given the parameters of the mixture (w and π), we have implemented two separate programs: one performing the maximizations with respect to l and α conditional on the current value of w and π for each alignment taken in turn, and returning the relevant expectations to another program, whose task is then to sum up the expectations obtained from all alignments, and perform the maximization with respect to w and π.

DATA, LEARNING AND TESTING PROCEDURE

Models were estimated on a subset of the HSSP database [START_REF] Sander | Database of homology-derived protein structures and the structural meaning of sequence alignment[END_REF]. This database contains more than 32,012 alignments, corresponding to globular proteins for which the three dimensional structure is known. Many alignments contain a large number of sequences, which we assume might help inferring more accurate amino-acid replacement models. On the other hand, most alignments contain a high proportion of missing data. Furthermore, there is a high level of redundancy, many sequences being represented more than one time in the database. We therefore proceeded to a cleaning process (Supplementary Information), so as to retrieve a subset of 1,030 ungapped and non-redundant alignments (with an average of 40 sequences and 253 aligned positions per alignment).

We also tried an alternative database, HOGENOM release 3 [START_REF] Dufayard | Tree pattern matching in phylogenetic trees: automatic search for orthologs or paralogs in homologous gene sequence databases[END_REF]. Each alignment of HOGENOM was scanned with Gblocks [START_REF] Castresana | Selection of conserved blocks from multiple alignment for their use in phylogenetic analysis[END_REF], using the default options, and 1,200 alignments, with number of taxa ranging from 15 to 50, were selected at random.

For each alignment, the topology of the phylogenetic tree was inferred under the WAG model (Whelan and Goldman, 2001), using PhyML (Guindon and [START_REF] Guindon | A simple, fast and accurate algorithm to estimate large phylogenies by maximum likelihood[END_REF]. We used the EM algorithm mentioned above to estimate profile models, with K, the number of components of the mixture, ranging from 10 to 60, every 10. The corresponding model configurations will be referred to as C10, C20... C60. We used Xrate [START_REF] Holmes | An expectation maximization algorithm for training hidden substitution models[END_REF] and the same procedure as in Whelan and Goldman (2001) to estimate a single matrix model from the 1,030 HSSP alignments, which we named WAG HSSP.

The fit of the models was evaluated by estimating the maximum likelihood under 57 protein alignments of TreeBase [START_REF] Sanderson | The growth of phylogenetic information and the need for a phylogenetic database[END_REF]. For both matrix and profile models, the likelihood has to be maximized with respect to the topology, the branch lengths and the α parameter of the gamma distribution of rates across sites. The topology was optimized using the SPR search developed in PhyML [START_REF] Hordijk | Improving the efficiency of SPR moves in phylogenetic tree search methods based on maximum likelihood[END_REF], starting from the tree provided by TreeBase. In the case of profile mixture models, we can set the weights (w) associated to each component of the mixture equal to those inferred on the training database. In that case, all model configurations have the same number of parameters, and the likelihood score can then be directly compared between all models, including the WAG empirical matrix. Alternatively, we also tested the possibility of reoptimizing the weights on the test alignment, in which case the number of additional parameters, compared to WAG or to mixture models with fixed weights, is then equal to K -1. To account for such variations in model dimensionality, we used the Bayesian Information Criterion [START_REF] Schwartz | Estimating the dimension of a model[END_REF]:

BIC = -2 ln L + m ln N,
as well as the second order Akaike information criterion [START_REF] Akaike | A new look at the statistical model identification[END_REF]:

AIC = -2 ln L + 2m + 2m(m + 1) N -m -1 ,
where ln L is the maximum likelihood estimate under the model of interest, N is the number of aligned positions of the test dataset, and m is the number of parameters specific to the model under investigation (here, m = K -1).

The second order correction is necessary for small alignments. In any case, it automatically reduces to the standard first order AIC when N >> m. All AIC and BIC scores are displayed using WAG HSSP as the reference. Thus, for a given model M , the score is defined as:

∆AIC M = AIC W AG HSSP -AIC M ,
so that better models have higher scores. Significance of the difference in log likelihood was assessed, between alternative models with the same number of parameters, or between a priori specified alternative topologies under the same model, using the Kishino Hasegawa test [START_REF] Kishino | Evaluation of the maximum likelihood estimate of the evolutionary tree topology from DNA sequence data, and the branching order in Hominoidea[END_REF]. When two alternative models, M 1 and M 2 , return distinct ML trees T 1 and T 2 when applied on a given dataset, we may be interested in knowing whether the two trees are contained in the confidence sets of each model, which we assessed using the Shimodaira Hasegawa test [START_REF] Shimodaira | Multiple comparisons of log-likelihoods with applications to phylogenetic inference[END_REF]. Specifically, if we call L i1 and L i2 the log-likelihoods of the two trees under model M i , i = 1, 2, we assess the null-hypothesis that

E[L i2 ] = E[L i1 ].
Note that in the present case, where we test only two topologies, the Shimodaira Hasegawa test reduces to the Kishino Hasegawa test.

RESULTS

We used an Expectation-Maximization algorithm (see methods) to estimate by Maximum Likelihood a series of mixture models, with the number of components ranging from 10 to 60. We tried two alternative training databases: a subsample of HSSP [START_REF] Sander | Database of homology-derived protein structures and the structural meaning of sequence alignment[END_REF], and a subset of HOGENOM [START_REF] Dufayard | Tree pattern matching in phylogenetic trees: automatic search for orthologs or paralogs in homologous gene sequence databases[END_REF]. The resulting mixture models were then assessed by measuring their AIC and BIC score on a series of 57 alignments taken from TreeBase. We reasoned that testing on alignments that do not come from the same database as those that were used for estimating the parameters of the models avoids possible biases of the database. In addition, TreeBase contains alignments that have been produced especially for phylogenetic analyses, and thus, provide a priori a good benchmark for comparing models meant for phylogenetic reconstruction.

Reliability of the learning method

Focusing on the C20 model, we first evaluated the reliability and the reproducibility of the learning process. Five independent runs were performed, starting from random initial parameter configurations, and using the full training set (1,030 non-redundant alignments from HSSP). In addition, as a way of testing the effect of the finite size of the training set, we randomly split the training database into two subsets of 515 alignments, and estimated the parameters of the C20 model on each of them.

The independent runs performed on the full training database each time led to a distinct point of the parameter space, indicating that the likelihood surface has many local maxima. Such a complicated behavior was also observed in previous EM applications for Markov modulated models [START_REF] Holmes | An expectation maximization algorithm for training hidden substitution models[END_REF]. However, the profiles obtained from different runs are very similar. Specifically, we could unequivocally identify 18 of the 20 profiles across two runs taken at random (Tables 2 and3 of Supplementary Information). In addition, the differences in the AIC scores per site obtained when testing two independent runs against TreeBase are small (0.008), compared to the differences observed between the two independent training sets of 515 alignments (0.04), suggesting that effects due to the finite size of the database are dominant over problems related to the presence of local maxima. And more fundamentally, both sources of variability are smaller than the differences observed between different model configurations (see below).

Similar results were obtained under the C50 model: 41 of the 50 profiles could be unequivocally identified between two independent runs (Tables 4-7 of Supplementary Information), and a difference of 0.03 logarithmic units was observed between the resulting scores when testing on TreeBase. The difference was of 0.01 units when comparing the two independent training sets. Thus, for C50, the variations across independent runs and across independent subsets of the database are comparable, but still small compared to the differences between models.

The profiles obtained here confirm previous observations (Lartillot and [START_REF] Lartillot | A Bayesian mixture model for across-site heterogeneities in the amino-acid replacement process[END_REF]. Importantly, they are sparse, in that most of them give a high probability for only 2 or 3 amino-acids, generally similar in their biochemical or geometrical properties. In higher dimensional models, such as C50, they also display some overlap and redundancy, and some profiles found for smaller models tend to split into more specialized forms in higher dimensional models. Some correlation can be found between the profile preferred at each site and other known structural attributes, such as exposure to solvent or, more weakly, secondary structure (Figures 1 and2 of Supplementary Information). However, these correlations are weak, suggesting that the constraints experienced by each site in a protein are both more complex and more local than what can be predicted based on such simple structural and functional correlates.

Model comparison

We compared the profile mixture models C10-C60 with the best current time-reversible amino-acid replacement matrix, WAG (Whelan and Goldman, 2001). In addition, to compare models trained on the same database, we estimated a single matrix model from the 1,030 HSSP alignments (see Methods). This matrix, which we named WAG HSSP, was used as our reference in the following model comparisons.

¢¡ £ ¥¤ ¦ ¨ § © ¤ ¤ ¤ ¤ ¤ ¦ ¨ § © ¤ ! " $# % & '& ( ¤ ( ¤ ( ¤ ) ¤ ¤ ) ) 0
Fig. 1. AIC and BIC scores per site, using WAG HSSP as the reference, as a function of the number of profiles in the mixture. Weights were either fixed (Fix) or optimized (Opt).

The WAG HSSP matrix is clearly superior to JTT and WAG (Figure 1). The mean gain is 0.210 and 0.082 AIC/site, and a higher likelihood is obtained for 47 and 48 alignments out of 57, respectively. There are many reasons for the difference between these three matrix models among which the size and the quality of the training database, and the learning method [START_REF] Le | An improved general amino-acid replacement matrix[END_REF].

We first consider the mixture models with the weights of the components fixed to the values estimated along with the other parameters (Figure 1, Fix.). In this situation, all models have the same number of parameters, so that AIC and BIC scores are equivalent to twice the logarithm of the likelihood.

We see a clear improvement of the fit of the model as the number of categories increases (Figure 1). In fact, even with 60 categories, we did not reach the point where the fit would start to decrease. Among all mixture configurations, only C10 is worse than both WAG and WAG HSSP. Thus, C10 may be too simple a model compared to single matrix models such as WAG or WAG HSSP. In contrast, all other mixture models, from C20 to C60 clearly outperform both matrices.

Interestingly, the highest improvement in fit is accomplished when going from C10 to C20 (a gain of more than 0.3 AIC per site), whereas the difference between C20 and the best model tested here, C60, is of less than 0.2 AIC per site (Figure 1). This suggests that 20 profiles (but not 10) may be sufficient to capture the essential aspects of across site heterogeneities in amino-acid propensities.

If, instead of keeping the weights of the mixture fixed, we try to reoptimize them on the test alignments (Figure 1, Opt.) the conclusion are ambiguous: whereas AIC tends to favor weight reoptimization, BIC draws the opposite conclusion. In fact, according to BIC, not only do mixture models with reoptimized weights have a poorer fit, compared to their counterparts with fixed weights, but in addition, they are increasingly less fit as the number of components increases between 20 and 60. How to interpret this disagreement between AIC and BIC is not clear. A likelihood ratio test (LRT) can be performed, for assessing whether weights should be fixed or reoptimized for a given number of components. By this test, the fixed-weight null hypothesis is rejected in favor of weight reoptimization in about two thirds of the alignments at the 0.05 level (Table

SAT > 2 8 3-1 3-1 4-1 4-1 5-1 5-1 ¢¡ ¤£ ¦¥ ¨ § ©¡ ¤ £ ¥ § © ¤ ¢ ¢ ¢ ¢ ! £ ¦¥ § "¡ ¤ # $ # $ # $ ¢ ¢ % & $ % & $ % ¢' $ % & $ % ¢( ¢ % & $ ) 0 1 2 3 4
Fig. 2. Correlation between the AIC gain per site and the saturation level. We used as reference the mean AIC score per site obtained under WAG HSSP, and averaged over the subset of alignments with given saturation index 8, Supplementary Information), thus suggesting that BIC puts too strong a penalty on the additional parameters represented by the weights. On the other hand, according to LRT, one third of the alignments still prefer the weights to be fixed to their predefined values. For simplicity, in the following, we focus exclusively on mixture models with fixed weights. Finally, as an alternative to HSSP, we tried another database, HOGENOM [START_REF] Dufayard | Tree pattern matching in phylogenetic trees: automatic search for orthologs or paralogs in homologous gene sequence databases[END_REF], for estimating a complete series of profile mixtures. These mixtures where then tested against TreeBase, and were found to lead to a poorer fit in all cases (Figure 3 of Supplementary Information). In the following, we only consider the profiles obtained on HSSP.

Mixture models are better on saturated data

Not all alignments prefer mixture models with more than 20 profiles over matrices. For instance, C20 was found significantly better than WAG HSSP on 17 alignments out of 57 (Table 1), that is, 30% of the alignments. This proportion increases along with the number of components. Yet, even under the richest model proposed in the present work, at least 21% of the alignments seem to favor single matrix models over mixture models, among which 5 (9%) do so significantly. This raises the question of which aspects of the data are important in determining this preference.

Previous analyses have suggested that mixture models are especially better in the presence of saturation, i.e., when many sites of the protein have undergone a high level of multiple substitutions over the phylogenetic history. To investigate a potential relationship between fit and saturation in the present context, we computed a [START_REF] Lartillot | Suppressing long branch attraction artefacts in the animal phylogeny using a site-heterogeneous model[END_REF]. We observe that, as more and more saturated data are analyzed, the gain in fit obtained by profile mixtures is indeed higher (Figure 2, see also Figure 6 of Supplementary Information), and similarly, that the proportion of datasets for which profile mixtures are significantly better than WAG HSSP is larger (Table 1). In particular, mixture models with 20 profiles or more were found significantly worse on only one aligment with a saturation index greater than 2. Altogether, this experiment confirms that profile models are better at describing saturated data.

Inferring the topology of the phylogenetic tree

We implemented WAG HSSP and the profile mixture models in the program PhyML [START_REF] Guindon | A simple, fast and accurate algorithm to estimate large phylogenies by maximum likelihood[END_REF]. For each alignment of our test set, we searched for the maximum likelihood topology under WAG, WAG HSSP, C20 and C50.

A different topology is found by WAG and by C20 in 51 out of 57 cases, with significance achieved in 4 cases. Similar observations were made when comparing WAG and C50, WAG and WAG HSSP, or C20 and C50 (Table 9 of Supplementary Information). Note that the present statistical test is different from the previous one (Table 1). In Table 1, we were comparing the maximum likelihood scores obtained by the two models, each under its preferred topology. In the present case, we assess the two alternative trees under the same model, which is a more stringent assessment.

Previous analyses, based on phylogenomic datasets, have suggested that profile mixtures may be more robust to long branch attractions. To check that the profile mixtures developed here have similar properties, we analyzed a series of phylogenomic datasets taken from previously published phylogenomic studies [START_REF] Brinkmann | An empirical assessment of long-branch attraction artefacts in deep eukaryotic phylogenomics[END_REF][START_REF] Philippe | Multigene analyses of bilaterian animals corroborate the monophyly of Ecysozoa, Lophotrochozoa and Protostomia[END_REF]. These datasets have been built by concatenating single gene alignments, and display well documented phylogenetic artifacts. Based on the literature, we can thus define for each dataset two alternative topologies: one, obtained using most currently available phylogenetic models, including WAG, and likely to be artifactual; and another one, thought to be closer to the true topology. In the following, we will refer to them directly as the 'artifactual' and the 'true' topology (see Supplementary Information).

We first compared the models. As can be seen from Table 2, all phylogenomic datasets display a strong preference for profile mixtures (C20 and C50), over WAG HSSP (p < 10 -4 for the Kishino Hasegawa test in all cases). They also have a high saturation index 2), further indicating the correlation between high saturation and good fit for profile models.

Next, using the two alternative 'true' and 'artifactual' trees, we performed a Kishino Hasegawa test, successively under WAG HSSP, C20 and C50. In all cases (Table 3), the WAG HSSP matrix prefers the artifactual trees, and the true tree is significantly rejected at the 0.05 level in one case (microsporidians). In contrast, C50 always prefers the true topology, with significant rejection of the alternative tree achieved in two cases (microsporidians and Guillardia's nucleomorph). In one case, namely the position of microsporidians, there is a complete reversion from a significant rejection of the 'true' tree by WAG, to a significant rejection of the 'artifactual' tree under C50. A basal position of microsporidia, as favored by WAG, is probably the result of a long branch attraction artifact (LBA) [START_REF] Keeling | Microsporidia: biology and evolution of highly reduced intracellular parasites[END_REF][START_REF] Brinkmann | An empirical assessment of long-branch attraction artefacts in deep eukaryotic phylogenomics[END_REF]. Consequently, the reversion observed under C50 can be interpreted as a case of suppression of the LBA by the use of a more adequate model. Concerning C20, the results are similar: for all datasets but one (nematodes), the true topology is preferred, which suggests that C20 is a sufficiently rich model for phylogenetic reconstruction.

DISCUSSION AND CONCLUSION

In this work, we have developed and tested an empirical profile mixture model. Our aim was two-fold. First, we wanted to transpose some aspects of earlier work initially conducted in a Bayesian framework within a Maximum Likelihood paradigm. But also, we wanted to dispense with the need of inferring all the parameters of the profile mixture directly on the dataset under investigation, as it precludes the analysis of small alignments. Accordingly, we tried to provide phylogeneticists with a model whose parameters are predefined, and should be suitable for a large spectrum of alignments of phylogenetic interest.

Our analysis indicates that a moderate number of profiles, around 20, is enough to accomplish essential improvements in two directions, compared to standard amino-acid replacement matrices. First, it provides a better fit in many cases, and second, it reconstructs more reliable phylogenetic trees, with a lesser sensitivity to long-branch attraction artifacts. In particular, the striking reversion observed in the case of microsporidians illustrates how relaxing the assumption of substitutional homogeneity along the sequence can dramatically alleviate systematic errors in phylogenetic reconstruction. Similar results have been obtained previously in a Bayesian context, based on a more complex non-parametric profile mixture model. In this respect, the present work shows that those results are robust, with respect to the details of the model, the implementation, and the statistical framework.

We also observe that, in a non negligible fraction of test alignments, the empirical matrix remains the best alternative. In a theoretical perspective, this could mean that, if profile mixtures account for some fundamental aspects of the amino-acid replacement patterns that empirical matrices fail at capturing, on the other hand, they seem to miss something that these matrices correctly describe. As was previously suggested [START_REF] Lartillot | Suppressing long branch attraction artefacts in the animal phylogeny using a site-heterogeneous model[END_REF], and further confirmed in the present analysis, what profile mixtures are good at modelling is the fact that some sites undergo repeated replacements among very small subsets of the amino-acid alphabet, something which is not easily encoded into one single global matrix. Conversely, we observe here that empirical matrices still better descibe the amino-acid replacement process at shorter evolutionary distance. This suggests that mixtures of matrices could be even better than mixtures of profiles, provided that the equilibrium frequency vectors of some of those matrices are concentrated on small subsets of amino-acids.

Note that progress is still possible in the way empirical matrices are estimated. In this direction, a better empirical matrix was recently obtained by [START_REF] Le | An improved general amino-acid replacement matrix[END_REF]. In the present context, we observe that it provides a significantly better fit than WAG HSSP in most cases (averaged AIC score difference of 0.10 over Tree-Base). However, the essential result obtained here, namely, that empirical profile models are better on saturated data, is not fundamentally changed if this new matrix is used instead of WAG HSSP.

In a more practical perspective, all these observations raise the question of how to decide for the best model in individual cases. A straightforward method would consist in systematically running phylogenetic analyses under the two models, WAG HSSP and C20, in parallel, and then choosing the best model on a case by case basis. As we have shown above, a simple alternative to this systematic duplication would consist in computing the saturation index of the dataset based on a maximum parsimony reconstruction, and choose C20 whenever this saturation index is too high. Otherwise, WAG HSSP, or another empirical matrix, should be used.

An important feature of the profile mixtures is their overall complexity. We did not investigate mixtures with more than 60 profiles, yet the likelihood scores indicate that higher dimensional mixtures would still bring further improvement. On the other hand, the resulting profiles would probably display some redundancy, as is already observed under C60. A possible explanation of this phenomenon is that the variations across sites form a continuum, imperfectly described by a finite mixture. In such a situation, arbitrarily rich mixtures can be favored given a sufficient amount of data, and the decision to stop has to be made based on a tradeoff between tractability and usefulness [START_REF] Steel | Should phylogenetic models be trying to 'fit an elephant[END_REF]. In the present case, the C20 model represents such a pragmatic tradeoff, representing about as many parameters as 2 time-reversible empirical matrices, while achieving a substantial improvement over one matrix.

The profile mixture models developed in this paper are suitable for both ML and Bayesian analyses. Accordingly, we have developed them in the two frameworks in parallel, by simultaneously implementing them in two phylogenetic reconstruction programs, PhyML [START_REF] Guindon | A simple, fast and accurate algorithm to estimate large phylogenies by maximum likelihood[END_REF] and PhyloBayes (http://atgc.lirmm.fr/cat). In the ML framework, C20 probably provides the best practical compromise between fit and computational efficiency. Thanks to a recoding trick (possible only for F81 processes), the C20 model represents an overall 4 to 5 fold increase, both in computational time and memory requirements, instead of the 20-fold increase that would be obtained without the recoding.

In the Bayesian context, any model from C20 to C60 can be used, with only marginal differences in the resulting computational load, thanks to Monte Carlo parameter expansion methods. Empirical mixtures are especially suited to phylogenetic analyses using small alignments (e.g. single gene analyses), whereas large analyses should preferably be conducted under the more flexible non-parametric CAT model developed previously.

Altogether, our two implementations make it possible to perform extensive comparisons between matrix and vector models, under two distinct statistical frameworks, and using both singlegene alignments or large concatenations. Such a generalized crosstalk between alternative models and methods will greatly enlarge the experimental perspective, and will hopefully provide important clues as to what further improvements are still needed for converging towards more reliable phylogenetic trees.

Table 1 .

 1 Number of alignments with significantly better (first number) or significantly worse (second number) likelihood scores under each model, compared to WAG HSSP, as a function of saturation index. Significance is assessed by the Kishino-Hasegawa test, with a threshold p-value of 0.01.

		#Aln C10 C20 C30 C40 C50 C60
	all	57	5-13 17-6 18-6 19-6 24-5 27-5
	SAT < 1	35	1-10 6-4	6-4	7-4	8-3 11-3
	SAT > 1	22	4-3 11-2 12-2 12-2 16-2 16-2

Table 2 .

 2 Number of sites, saturation index (SAT) and AIC score per site on genomic data sets

		# sites SAT AIC / site
		C20 C50
	Metazoans (first half)	17,136 1.96 0.61 0.91
	Metazoans (second half) 17,135 2.14 0.60 0.91
	Nematodes	35,366 1.45 0.31 0.52
	Microsporidia	24,294 2.55 0.62 0.95
	Nucleomorph	24,294 2.51 0.63 0.89
	mean saturation index (SAT ) for each of the 57 alignments of the
	test database. This saturation index is essentially the average mini-
	mum number of convergence or reversions per site, as inferred by a
	Maximum Parsimony reconstruction (see

Table 3 .

 3 Topology comparison on genomic data sets. ∆ ln L is between true and artifactual topologies, under the specified model

		WAG HSSP	C20		C50
		∆ ln L	p	∆ ln L	p	∆ ln L	p
	Metazoans (1)	-27	0.26	14	0.29	16	0.26
	Metazoans (2)	-67	0.07	1	0.49	13	0.29
	Nematodes	-106 0.07	-2	0.48	7	0.43
	Microsporidia	-180 0.00	144	0.00	147	0.00
	Nucleomorph	-9	0.43	62	0.00	63	0.00
	(Table

ACKNOWLEDGEMENT AND FUNDING

We wish to thank Nicolas Rodrigue and the anonymous referees for their useful comments on the manuscript. This work was financially supported by the ANR-BIOSYS MITOSYS project and by the ACI-IMPBIO ModelPhylo project.


A Dereeper

V Guignon

G Blanc

S Audic

S Buffet

F Chevenet

J.-F Dufayard

S Guindon

V Lefort

M Lescot

J.-M Claverie

O Gascuel

Phylogeny.fr: Robust Phylogenetic Analysis for the Non-Specialist

Phylogenetic analyses are central to many research areas in biology and typically involve the identification of homologous sequences, their multiple alignment, the phylogenetic reconstruction and the graphical representation of the inferred tree. The Phylogeny.fr platform transparently chains programs to automatically perform these tasks. It is primarily designed for biologists with no experience in phylogeny, but can also meet the needs of specialists; the first ones will find up-todate tools chained in a phylogeny pipeline to analyze their data in a simple and robust way, while the specialists will be able to easily build and run sophisticated analyses. Phylogeny.fr offers three main modes. The "One Click" mode targets non-specialists and provides a ready-to-use pipeline chaining programs with recognized accuracy and speed: MUSCLE for multiple alignment, PhyML for tree building, and TreeDyn for tree rendering. All parameters are set up to suite most studies, and users only have to provide their input sequences to obtain a ready-to-print tree. The "Advanced" mode uses the same pipeline but allows the parameters of each program to be customized by users. The "A la Carte" mode offers more flexibility and sophistication, as users can build their own pipeline by selecting and setting up the required steps from a large choice of tools to suite their specific needs.

Prior to phylogenetic analysis, users can also collect neighbors of a query sequence by running BLAST on general or specialized databases. A guide tree then helps to select neighbor sequences to be used as input for the phylogeny pipeline.

INTRODUCTION

Reconstructing the evolutionary history of molecular sequences through phylogenetic analysis is at the heart of many biological research areas such as comparative genomics, functional prediction, detection of lateral gene transfer, or the identification of new micro-organisms. Starting from a sequence of interest, a typical phylogenetic analysis goes through successive steps that include the identification of homologous sequences, multiple alignment, phylogenetic reconstruction and graphical representation of the inferred tree. This process requires substantial computational resources depending on the number and length of the sequences, and on the methods being run.

A huge variety of models, approaches and computer programs are currently available, as can be seen from Joe Felsenstein's phylogeny software inventory 1 . The task of deciding which method to use, installing the corresponding programs, and running them on a local computer, is beyond the reach of most occasional users. Yet, phylogenetic trees have become a compulsory illustration (and referee's request) in most sequence-related studies. As a consequence, user-friendly but ancient programs are still widely used, though much improved methods now exist and are only used by the specialists, typically involved in molecular evolution or systematics.

In this context, Phylogeny.fr has been designed to provide a ready-to-use platform that transparently chains alignment and phylogeny programs in a comprehensive and flexible manner. Although phylogenetic specialists will be able to find up-to-date tools and run sophisticated analyses based on their favorite approaches and their own parameter settings, the primary philosophy of Phylogeny.fr is to assist biologists with no experience in phylogeny in analyzing their data in a simple and robust way, using methods corresponding to well accepted standards. Maximum-likelihood (ML) tree construction is the default option to infer phylogenies, which is commonly recognized (1) as the most accurate approach (along with Bayesian) in molecular phylogenetics.

Phylogeny.fr offers "one-stop-shopping" among a variety of leading methods for multiple sequence alignment, phylogenetic reconstruction and graphical representation of trees, and chains these methods into a pipeline that can be executed in three modes. The "One Click" mode is designed for biologists with no experience in bioinformatics; given a set of unaligned sequences, a predefined pipeline using MUSCLE (2), Gblocks (3), PhyML (4) and TreeDyn (5) outputs the corresponding phylogenetic tree in a ready-to-print format. The "Advanced" mode allows the settings of each "One Click" tool to be customized by users. The "A la Carte" mode offers flexible choices in the pipeline steps, the tools and their settings to suite the more specific needs of experts.

Several other websites propose related services. PhyloBuilder [START_REF] Glanville | Berkeley Phylogenomics Group web servers: resources for structural phylogenomic analysis[END_REF] and PhyloBlast [START_REF] Brinkman | PhyloBLAST: facilitating phylogenetic analysis of BLAST results[END_REF] are dedicated to proteins and gather homologs of a query sequence to build a phylogenetic tree using a distance or a parsimony method. POWER (8) infers a phylogenetic tree from a sequence set using a pipeline involving ClustalW [START_REF] Thompson | CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice[END_REF] and PHYLIP [START_REF] Felsenstein | PHYLogeny Inference Package) version 3.6a2, distributed by the author, department of genetics[END_REF] programs. Tarraga et al. recently proposed Phylemon [START_REF] Tarraga | Phylemon: a suite of web tools for molecular evolution, phylogenetics and phylogenomics[END_REF] which provides experts with a suite of on-line programs and a Java interface to build a phylogeny pipeline. The main specificity of Phylogeny.fr is the combination of an interface designed for the non-specialists with up-to-date programs that are often reserved to experts. Moreover, Phylogeny.fr is able to analyze both DNA and protein sequences.

OVERVIEW

The Phylogeny.fr platform proposes three major components:

(i) a pipeline to reconstruct a phylogenetic tree from a set of sequences through an automated process that successively performs multiple sequence alignment, alignment refinement, phylogenetic reconstruction and, finally, a graphical representation of the resulting tree;

(ii) a fast parallel Blast (12) module to search for similar sequences of a query sequence. Sequence selection is facilitated by a "quick-and-dirty" guide tree based on Blast results, and by an estimate of final alignment length;

(iii) a suite of stand-alone phylogenetic programs.

All these components are offered in a user-friendly tabulated menu facilitating navigation, together with examples to familiarize users with the correct input and expected results. All programs used in the pipeline and their parameters and references are displayed in an "Analysis overview" that provides the main information to be indicated in a publication. In addition, the web interface enables users to visualize, manipulate, and locally edit their results through applet viewers such as Jalview [START_REF] Clamp | The Jalview Java alignment editor[END_REF] for alignment, or ATV [START_REF] Zmasek | ATV: display and manipulation of annotated phylogenetic trees[END_REF] for phylogenetic tree.

PHYLOGENY PIPELINE

At the core of the system lies the pipeline allowing a phylogenetic tree to be built at once from a set of sequences. This pipeline consists of a succession of Perl modules wrapping different external software programs, which can be executed through the three modes described below. This process requires many additional programs that are not detailed here: [START_REF] Kuhner | A simulation comparison of phylogeny algorithms under equal and unequal evolutionary rates[END_REF] to make the junction between the main programs (format compatibility), [START_REF] Edgar | MUSCLE: a multiple sequence alignment method with reduced time and space complexity[END_REF] to output the results in several formats (e.g. PNG or PDF), or (3) to rearrange and modify the tree (e.g. Retree from PHYLIP is used to root the tree with the midpoint method or using an outgroup).

"One Click" Mode

The "One Click" mode targets users who do not wish to deal with program and parameter selection. By default, the pipeline is already set up to run and connect well recognized programs: MUSCLE for multiple alignment, Gblocks for automatic alignment curation, PhyML for tree building and TreeDyn for tree drawing.

Several studies showed that these programs are both fast and accurate. MUSCLE was assessed using several alignment databases, including the BAliBASE benchmark [START_REF] Thompson | BAliBASE: a benchmark alignment database for the evaluation of multiple alignment programs[END_REF], on which it achieved the highest ranking of any method at the time of publication. PhyML was shown to be at least as accurate as other existing phylogeny programs using simulated data, while being one order of magnitude faster. Both programs have continuously been improved. These two programs are widely used (~650 and ~1000 citations in Web of Science, for MUSCLE and PhyML, respectively). It follows that they are robust, stable and bug free, and that the best way to use them has been reported by their numerous users. Notably, the default options and parameter values have been selected with care; they are used in the "One Click" mode and should be suited for most studies.

PhyML is run with the aLRT statistical test ( 16) of branch support. This test is based on an approximation of the standard Likelihood Ratio Test, and is much faster to compute than the usual bootstrap procedure. Both methods output the same trees; branch supports may differ but are generally highly correlated. Altogether, this mode makes routine phylogenetic analyses simple and fast. Users only need to load their set of sequences and submit their request. A few minutes later, the system will display a ready-toprint phylogenetic tree in a publication quality image. One has still the possibility to disable or not the alignment refinement performed by the Gblocks program, which eliminates poorly aligned positions and divergent regions. It is also possible to specify an outgroup (midpoint rooting is the default option).

"Advanced" Mode

In the "Advanced" mode, the pipeline is the same as the one used in the "One Click" mode but users can now freely edit the settings of each program. The pipeline is flexible enough to enable users to select which steps to perform. In this manner, input data can be a set of non-aligned sequences, a sequence alignment or a tree in Newick format. Various file formats are supported as shown in Table 1. Furthermore, the system offers the possibility to inspect the results of each step before launching the next program, so that expert users can detect problems and properly adjust parameters accordingly.

"A la Carte" Mode

This mode presents the same interface as the "Advanced" mode (e.g., step-by-step execution, parameter settings) but offers the possibility of running, testing and comparing the efficiency of a larger panel of methods and programs. It proposes a user-friendly customizable pipeline on which users can select the steps to perform, and the program to use for each of them.

For the alignment step, this mode proposes T-Coffee [START_REF] Notredame | T-Coffee: a novel method for fast and accurate multiple sequence alignment[END_REF] which ranks high in accuracy but is slower than MUSCLE, and its derivative 3DCoffee (18) that offers the unique possibility of incorporating protein structure information to improve alignments.

For phylogeny reconstruction, a choice is offered among the following approaches: maximum likelihood method using PhyML, parsimony using NONA (for DNA, [START_REF] Goloboff | NO NAME) Version 2[END_REF] or Protpars (for proteins, 10), or distance-based methods using BioNJ [START_REF] Gascuel | BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data[END_REF] or Neighbor (NJ, 21) from PHYLIP [START_REF] Felsenstein | PHYLogeny Inference Package) version 3.6a2, distributed by the author, department of genetics[END_REF]. With BioNJ and NJ, evolutionary distances are calculated using Protdist [START_REF] Felsenstein | PHYLogeny Inference Package) version 3.6a2, distributed by the author, department of genetics[END_REF] for proteins and Fastdist [START_REF] Elias | Fast computation of distance estimators[END_REF] for DNA, while bootstrap analysis is achieved by combining Seqboot and Consense [START_REF] Felsenstein | PHYLogeny Inference Package) version 3.6a2, distributed by the author, department of genetics[END_REF]. Distance methods are very fast and should be preferred for large scale analyses, or when performing bootstrap studies. Parsimony and ML methods are roughly as fast in practice (at least when using the programs selected here); they are able to analyze relatively large datasets, but become quite slow with bootstrap. ML is commonly reported as the most accurate approach, but parsimony does not require any substitution model to be selected, which appears preferable to some users. In practice, it is always a good idea to run several programs with various options and verify the stability of the results. Phylogeny.fr greatly facilitates this conservative approach as users can easily modify the choices in the pipeline and run the server again.

The "A la Carte" mode makes it easy to integrate and compare both old reference programs -such as the basic tree viewers Drawgram and Drawtree, from PHYLIP -and new ones which are not yet well diffused in the scientific community -such as TreeDyn that offers sophisticated tools and options to draw trees. For example, this mode still provides ClustalW which is still the most popular alignment tool to date even though it has been shown to be less accurate than modern alignment programs [START_REF] Wallace | Multiple sequence alignments[END_REF].

Tree visualization and drawing

Following phylogenetic tree reconstruction by any of these three modes, an image of the tree is produced to facilitate its interpretation. Phylogeny.fr allows the tree image to be drawn and modified using a variety of TreeDyn options. These include: change of the tree shape (rectangular or radial), font style, color of the text and branches, display of branch supports, edition of taxon names, and root selection. The tree image can then be obtained in PNG and PDF formats and eventually incorporated in any artwork for publication.

SEQUENCE SEARCHING

Prior to a phylogenetic analysis, users can also build their initial dataset by running Blastall against several public sequence databases (e.g. 16S Database, GenBank NT, Swissprot). This search is performed in parallel on a 25-node Linux cluster, and allows for a quick exploration of the neighbors of a query sequence. A rough multiple alignment of the query and hit sequences is generated on the fly by piling the individual Blast pairwise alignments. This alignment serves to construct and display a NJ tree using uncorrected p-distances, so that sequences can be selected and incorporated into the dataset by simply clicking on the tree image. To facilitate this process, sequence names are colored according to the taxonomic group of the species they belong to. Furthermore, to help in the choice of sequences maximizing the number of gap-free sites available for phylogenetic reconstruction, an estimation of the length of the final multiple alignment is dynamically computed each time a sequence is selected. Once the selection step is completed, users can run the phylogeny pipeline or retrieve the selected sequences in FASTA format.

SERVER FEATURES AND LIMITATIONS

The platform currently runs on a dedicated server (PowerEdge 2850-Xeon 2.8GHz/2x2MB Dual Core), except for the Blast module which is parallelized on a 25-CPU cluster. Input limitations depend on the selected program and its computational speed, as shown in Table 1. MUSCLE and ClustalW are limited to 200 sequences, while T-Coffee and 3D-Coffee limitations are <50 sequences and <2000 sites. Distance-based phylogeny programs (i.e. NJ and BioNJ) have no limitation, while all other phylogeny programs are limited to: (sequence size) x (number of taxa) 2 < 10,000,000 (e.g. 100 sequences with length 1000). The number of bootstrap replicates -a highly sensitive parameter in terms of computational burden -is limited to 100 with all phylogeny programs, except NJ and BioNJ which are allowed for 500 replicates. These limitations should be relaxed in the near future, as the platform should be implemented on more powerful servers and multiple sites.

CONCLUSION AND FUTURE DEVELOPMENTS

To our knowledge, Phylogeny.fr is the first web server designed for both non-specialists and experts, which provides a complete automated phylogenetic analysis -from FASTA file to tree image -thanks to a pipeline of well trusted tools. Phylogeny.fr has been advertised in 2007 on a few diffusion lists. It currently processes about 2500 submissions per month from about 800 different users, who beta-tested the initial platform and now seem to be satisfied. Current waiting times with the "One Click" mode (without bootstrap) are ~2 and ~6 minutes for alignments of 50 DNA sequences with 1463 sites and of 40 proteins with 430 sites, respectively. When a bootstrap analysis is performed ("Advanced" and "A la carte" modes) waiting times are much longer (~100 times with 100 replicates) and the results are sent by email.

The modular architecture of the Phylogeny.fr pipeline will facilitate the addition of new features and new programs according to the evolution of the field. Phylogeny.fr is not committed to a particular set of programs. The methods constituting the backbone of the "One Click" mode today will be replaced by newer ones, if considered better by a large consensus in the field. The dethroned methods will then become part of the "A la carte" mode, to ensure compatibility with previous studies.

Among the further developments to come shortly, we will add the possibility of selecting the best substitution model to reconstruct the phylogeny, using MODELTEST (DNA, [START_REF] Posada | MODELTEST: testing the model of DNA substitution[END_REF] or ProtTest (proteins, 25).

Additional tools facilitating the selection of sequences after a BLAST search will be integrated, and the final output tree will become clickable to help recovering information on the sequences, e.g. functional annotations from the Gene Ontology. Moreover, we plan to provide password-protected accounts to enable users to store and load their results, and to refine their previous analyses; researchers will thus be able to create and share on-line projects with their collaborators.

Table 1 : Main tools available on Phylogeny.fr.

 1 : With fast programs Phylogeny.fr displays the results within a few minutes or less. Slow programs require about one hour or so, but their results are still displayed on-line, while the results of very slow program are sent by email.

	Program	Function	Input	Output	Speed Current limitations Use
	Blastall 2.2.17	Sequence searching	Raw, FASTA	FASTA	Fast	None	Advanced and A la Carte
	MUSCLE 3.6	Multiple alignment	FASTA,	FASTA, Clustal,	Fast	<200 sequences	All modes
			EMBL/Uniprot,	PHYLIP			Large dataset
			GenBank
	T-Coffee 4.97	Multiple alignment	FASTA,	FASTA, Clustal,	Very	<50 sequences	A la Carte
			EMBL/Uniprot,	PHYLIP, others	slow	<2000 sites	Small dataset
			GenBank
	3DCoffee 4.97	Multiple alignment	FASTA,	FASTA, Clustal,	Very	<50 sequences	A la Carte
		using structural	EMBL/Uniprot,	PHYLIP, others	slow	<2000 sites	Small dataset
		information	GenBank
	ClustalW 1.83	Multiple alignment	FASTA,	FASTA, Clustal,	Fast	<200 sequences	A la Carte
			EMBL/Uniprot,	PHYLIP			Large dataset
			GenBank
	Gblocks 0.91b	Alignment refinement	FASTA	FASTA, Clustal,	Fast	None	All modes
				PHYLIP			Large dataset
	PhyML 3.0	Phylogeny using	FASTA, Clustal,	Newick	Fast to	(sequence size) ×	All modes
		maximum likelihood	PHYLIP, EMBL,		Slow	(number of taxa) 2	Medium to large dataset
			PAUP*/Nexus			< 10,000,000
	NONA 2.0	DNA phylogeny using	FASTA, Clustal,	Newick	Fast to	(sequence size) ×	A la Carte
		parsimony	PHYLIP, EMBL,		slow	(number of taxa) 2	Medium to large datasets
			PAUP*/Nexus			< 10,000,000
	Protpars 3.66	Protein phylogeny	FASTA, Clustal,	Newick	Very	(sequence size) ×	A la Carte
		using parsimony	PHYLIP, EMBL,		slow	(number of taxa) 2	Small to medium dataset
			PAUP*/Nexus			< 10,000,000
	Dnadist Protdist	Phylogeny using	FASTA, Clustal,	Newick	Fast	None	A la Carte
	BioNJ	distances	PHYLIP, EMBL,				Large dataset
	NJ PHYLYP 3.66		PAUP*/Nexus
	Bootstrap with	Estimations of clade	FASTA, Clustal,	Newick	Very	<100 replicates	Advanced and A la Carte
	PhyML, NONA,	supports	PHYLIP, EMBL,		slow		Small to medium datasets
	Protpars		PAUP*/Nexus
	Bootstrap with	Estimations of clade	FASTA, Clustal,	Newick	Slow	<500 replicates	A la Carte
	distance methods	supports	PHYLIP, EMBL,				Medium to large datasets
			PAUP*/Nexus
	TreeDyn 198	Tree rendering	Newick	PNG, PDF, TGF	Fast	None	All modes
	Drawgram 3.66	Various tree shapes	Newick	PNG, PDF	Fast	None	A la Carte
		rendering
	Drawtree 3.66	Unrooted tree rendering Newick	PNG, PDF	Fast	None	A la Carte

Note

1: http://evolution.genetics.washington.edu/phylip/software.html

Abascal F, Zardoya R, Posada D. ProtTest: selection of best-fit models of protein evolution. Bioinformatics. 2005 May 1;21(9):2104-5.

ACKNOWLEDGEMENTS

Sincere thanks to the authors of the many tools used by Phylogeny.fr. This project was supported by fundings from the "Réseau National des Génopoles" (RNG).

CONTRIBUTIONS

IGS and MAB-LIRMM teams equally contributed to this work. AD, VG: conceived and implemented the server, wrote the manuscript; GB, SA, SB, FC, JFD, SG, VL, ML: helped in the design of the server and/or implemented specific tools; JMC: initiated the project, managed the alignment side, wrote the manuscript; OG: managed the phylogeny side, wrote the manuscript.