Nicolas Vergne
email: nvergne@genopole.cnrs.fr

Drifting Markov Models with Polynomial Drift and Applications to DNA Sequences *

Keywords: drifting Markov models, Markov models, DNA sequences, heterogeneity, rare Drifting Markov Models

In this article, we introduce the drifting Markov models (DMMs) which are inhomogeneous Markov models designed for modeling the heterogeneities of sequences (in our case DNA or protein sequences) in a more flexible way than homogeneous Markov chains or even hidden Markov models (HMMs). We focus here on the polynomial drift: the transition matrix varies in a polynomial way. To show the reliability of our models on DNA, we exhibit high similarities between the probability distributions of nucleotides obtained by our models and the frequencies of these nucleotides computed by using a sliding window. In a further step, these DMMs can be used as the states of an HMM: on each of its segments, the observed process can be modeled by a drifting Markov model. Search of rare words in DNA sequences remains possible with DMMs and according to the fits provided, DMMs turn out to be a powerful tool for this purpose. The software is available on request from the author. It will soon be integrated on seq++ library (http://stat.genopole.cnrs.fr/seqpp/).

Introduction

Modeling DNA sequences with stochastic models and developing statistical methods to analyze the enormous set of data that results from the multiple projects of DNA sequencing are challenging questions for statisticians and biologists. The most popular model in this domain is the Markov model on the nucleotides that gives a description of the local behaviour of the sequence (see [START_REF] Almagor | A Markov analysis of DNA sequences[END_REF], [START_REF] Blaisdell | Markov chain analysis finds a significant influence of neighboring bases on the occurrence of a base in eucaryotic nuclear DNA sequences both protein-coding and noncoding[END_REF], [START_REF] Phillips | The effect of codon usage on the oligonucleotide composition of the e. coli genome and identification of overand underrepresented sequences by Markov chain analysis[END_REF], [START_REF] Gelfand | Extendable words in nucleotide sequences[END_REF]).

Thanks to the statistical properties of these Markov models, we can enlighten different biological properties of DNA or protein sequences. Different Markov models may be proposed. First of all, classical homogeneous with some order k Markov chains provide a general description of a sequence (for instance, the different frequencies of the dinucleotides). [START_REF] Simons | Global Markov models for eukaryote nucleotide data[END_REF] provides a good discussion about this global Markov model. See also [START_REF] Almagor | A Markov analysis of DNA sequences[END_REF] or [START_REF] Blaisdell | Markov chain analysis finds a significant influence of neighboring bases on the occurrence of a base in eucaryotic nuclear DNA sequences both protein-coding and noncoding[END_REF] on this point. [START_REF] Schbath | Exceptional motifs in different Markov chain models for a statistical analysis of DNA sequences[END_REF] identifies exceptional motifs in sequences using Markov models. [START_REF] Reinert | Compound Poisson and Poisson process approximations for occurrences of multiple words in Markov chains[END_REF] gives another way to detect rare words in biological sequences and [START_REF] Nuel | Grandes déviations et chaînes de Markov pour l'étude des occurrences de mots dans les séquences biologiques[END_REF] proposes a comparison between the most used methods for discovering relevant patterns in sequences modeled by classical Markov chains.

The main drawback in considering classical Markov models for the analysis of sequences is that it supposes the homogeneity of sequences, whereas it turns out that long biological sequences are inhomogeneous. A way to take into account this heterogeneity is the use of hidden Markov models (HMMs). HMM is largely used for modeling biological sequences. For instance [START_REF] Churchill | Stochastic models for heterogeneous DNA sequences[END_REF] analyzes the heterogeneity of DNA sequences using HMMs. See for example [START_REF] Stanke | Gene prediction with a hidden Markov model and a new intron submodel[END_REF] or [START_REF] Krogh | A hidden Markov model that finds genes in escherichia coli DNA[END_REF] for applications to gene prediction. Thanks to HMMs, one can detect coding or non-coding regions, exons or introns, but also homologies between sequences or discover horizontal transfers [START_REF] Nicolas | Mining bascillus subtilis chromosome heterogeneity using hidden Markov models[END_REF]. HMM corresponds to the biological fact that some signals succeed one another along the sequence. For example, on a DNA strand, a gene may be followed by a non-coding region, then by a promoter, an other gene and so on. Proteins are often composed of various "domains" separated by hinge. It is natural to think that the way the letters succeed differs from one of these regions to the others, and this explains the success of HMM in the search of regions with different biological roles.

Nevertheless, it is common to observe gradual variations along a biological sequence, either at a global level, either within one of the regions we just mentioned. For example, the gc-richness of a sequence varies according to the position. A first model refers to two kinds of behaviours: high percentage of gc (denoted by H) against low percentage of gc (denoted by L). Then a refined model has been developed, introducing H1, H2, H3, and L1, L2 regions. But there is a broad consensus about the simplifying aspect of this model, in particular its inability of a sharp determination of the limits of the regions belonging to one of the two (or five) categories: a soft transition from a gc-richness to another one is always observed. As an example, we will use our model on Phage lambda complete genome (see [START_REF] Wu | Nucleotide sequence analysis of DNA. II. Complete nucleotide sequence of the cohesive ends of bacteriophage lambda DNA[END_REF]). Figures 3(a), 3(b), 3(c) and 3(d) show an estimation of the richness of each of the four nucleotides as a function of the position along this genome (this estimation was obtained using a sliding window of width 2000). The figure shows that "at each position" at least one of the 4 curves has a soft variation. Even around positions 22000, where 3 richness curves seem to have a discontinuity (smoothed by the usage of a window), the fourth one (corresponding to the a nucleotide) has a continuous variation. Even inside genes, for example, this type of behaviour is observed (see [START_REF] Nicolas | Mining bascillus subtilis chromosome heterogeneity using hidden Markov models[END_REF]).

It is then necessary to develop mathematical tools to account for such gradual changes and we propose such a model, the drifting Markov model (DMM, see precise definition below). It can be seen as a competitive model to the HMM one: a DMM can be adjusted to a whole sequence; and it turns out that the classical problem of the search of rare words remains tractable with this model. But it over all can be understood as a complementary tool: the hidden models of an HMM, usually fixed Markov chains, can be replaced by DMM. This second approach will be treated in a further paper, the present one presenting the necessary tool and first results about its ability for the modeling of biological sequences.

Walking Markov models (WMM), introduced by [START_REF] Fickett | Base compositional Structure of Genomes[END_REF] were the first models with a continuous change of base composition. They want to model gc and at composition in a DNA sequence as we just discussed above isochores. For example, they cut a sequence in 1000-base windows and estimate a Markov model on all the windows containing between 300 and 400 at, 400 and 500 at, 500 and 600 at, 600 and 700 at, to have four Markov models. Then for any value w (the at-content), a Markov model M w is defined by linear interpolation of these primary processes. At last, WMM is defined by a random walk on w: they choose an initial value for w between 1/3 and 2/3 (that changes according to the studied sequence), and to choose each succeeding base, they add or subtract (with probability 0.5) 0.0015 from w and use M w to generate the next base. We use a totally different way to define our DMMs. First, we do not use a random walk to choose our transition matrix: our models are based on the sequence. Second, our models are adapted for any size of state space without a lot of preliminary treatments such as the estimation of some Markov models. It would be difficult to adapt WMM to state space of size 20. Of course, WMMs, just as well as DMMs, do not model detailed local structure, such as the local structure of genes. They are intended to model the large-scale background variation of base composition in the genome.

At now, let us explain the principle of a DMM. Instead of fitting a transition matrix on a whole sequence (homogeneous Markov model) or different transition matrices on different homogeneous parts of the sequence (HMM), we allow the transition matrix to vary (to drift) from the beginning to the end of the sequence. At each position, we obtain a different transition matrix. Our models are thus constrained heterogeneous Markov models. In this paper, we focus on a polynomial drift. The use of such models, where the transition matrix on DNA-alphabet or protein-alphabet (state space) may vary along the genome is a completely new approach.

In the second part of this work, the correct adjustment of probability distributions of nucleotides in DMMs to nucleotide frequencies computed on real sequences shows that our new models provide a more flexible, higher-dimensional parameterization of the data that can be hoped to result in better fits than homogeneous Markov models or HMMs (see Figures from 3(a) to 3(d) and from 4(a) to 4(h)). Then, we compute some model selection criteria (AIC and BIC) to compare different models. Two applications of our models are proposed here. Relying on the compositional asymmetries between the leading and the lagging strand of replication, the program ORILOC [START_REF] Lobry | Oriloc: prediction of replication boundaries in unannotated bacterial chromosomes[END_REF] helps to predict replication origins in bacterial genomes. We propose an alternative method based on our modeling to detect replication origins which present the advantage of being able to compute analytically a maximum. At last, we discuss a new application for the search of rare words in sequences modeled by a DMM. We offer a simple example with the Chi (gctggtgg motif) of Escherichia coli and we give different classifications of words according to different models. Many papers treat of rare words and patterns in biological sequences modeled by Markov chains [START_REF] Schbath | Exceptional motifs in different Markov chain models for a statistical analysis of DNA sequences[END_REF][START_REF] Reinert | Compound Poisson and Poisson process approximations for occurrences of multiple words in Markov chains[END_REF]), but all of them are based on Markov chains and their homogeneity. We offer the possibility to study rare words with a model which better correspond to the real sequence, so we can assume the reliability of our result in a better way than before.

This paper is organized in the following way. In Section 2, we describe the drifting Markov models with polynomial drift. Different methods of estimation are proposed and explained. In Section 3, we give first results concerning these new models. We establish reliability of DMMs by adjusting probability distributions of nucleotides and nucleotide frequencies (Section 3.1). We compare different models using AIC and BIC. We propose an alternative to the software ORILOC [START_REF] Lobry | Oriloc: prediction of replication boundaries in unannotated bacterial chromosomes[END_REF] for detecting replication origins (Section 3.3) and another application to the search of rare words in DNA sequences (Section 3.4). At last, in Section 4, we discuss our results and offer perspectives about these models.

Drifting Markov models

A sequence modeled by an order k drifting Markov model is a sequence of random variables X t : X = (X t ) t∈{0,...,n}

where n + 1 is the length of the sequence and where instead of fitting only one transition matrix on the whole sequence, we fit a possibly different transition matrix at each position in the sequence. Hence, the distribution of X t is defined in the following way:

P (X t = v|X t-k . . . X t-1 = u) = Π t n (u, v) with u = u 1 u 2 . . . u k a k-word and (u 1 , u 2 , . . . , u k , v) ∈ A k+1
where A is the state space (the alphabet A = {a,c,g,t} for example). Drifting Markov models are inhomogeneous Markov models and without constraints they cannot be estimated. Thus we propose a polynomial evolution of the transition matrix, according to the position in the sequence. We begin by using a linear drift (and later we will more generally use polynomial drifts).

Drifting Markov models: linear drift

We fix a transition matrix Π 0 at the beginning of the sequence and a transition matrix Π 1 at the end of the sequence and we allow the transition matrix to vary linearly from Π 0 to Π 1 :

Π t n = 1 - t n Π 0 + t n Π 1 .
Polynomials (1 -t/n) and t/n are chosen to establish the stochasticity of Π 0 to Π 1 . Obviously, role of Π 0 to Π 1 is artificial as any model parameters but stochastic matrices make easier the understanding of the model. We want to estimate these two matrices in order to build the model. In the case of a simple Markov model, the method of maximum likelihood is successfully used but, because of numerical complexity, we cannot use it here. Hence, we propose two different methods to estimate Π 0 and Π 1 : a matrix regression method and a point by point method. We describe these two methods in the following subsections. We just give the example of an attempt of likelihood maximisation for a DMM of order 1 and degree 1, to conclude that numerical complexity precludes the use of the estimation by maximum likelihood. The likelihood ℓ of a DMM of order 1 and degree 1 is the product over t of the transition matrices Π t n (X t-1 , X t ). Then the log-likelihood L is:

L(X, Π 0 , Π 1 ) = log ν 0 (X 0 ) + n t=1 u∈A 1 {X t-1 =u} v∈A 1 {Xt=v} log Π t n (u, v) .
In order to obtain the maximum likelihood, we look for the zero of the derivative of L. We obtain a system of 2|A|(|A| -1) equations with 2|A|(|A| -1) variables. In fact, it reduces to |A| systems of 2(|A| -1) equations with 2(|A| -1) variables. In Appendix A, we give an example of one of these systems with alphabet A = {a,c,g,t}, in order to see that all the variables are in the denominator and that it is unthinkable to solve these systems by analytical or numerical methods. Obviously the same problem exists with polynomial DMMs of higher order or higher degree. It is sheer madness to envisage to solve numerically such a system.

Estimation by a matrix regression method.

A first idea to obtain estimators of the matrices Π 0 and Π 1 is to divide the sequence in N segments of size m (size m will be chosen later). The idea of this method is to use an "approximated homogeneity" on each segment. Then, we fit on each segment S ℓ , for 1 ≤ ℓ ≤ N, a Markov model Π S ℓ classically estimated by maximum likelihood estimation. In order to fit our heterogeneous model on the whole sequence, we choose one point in each segment. We choose the N centers τ ℓ of the N segments S ℓ because E( Π S ℓ ) tends to Π τ ℓ as m goes to infinity. We could choose more than one point by segment but that induces numerical complexities without improving the estimation. We want our matrix Π t n to be the nearest possible to each Π S ℓ at the center of each segment S ℓ . Thus, for the matrix regression, we minimize the sum of distances between the estimated matrices on each segment S ℓ and the transition matrices Π t n at the center τ ℓ of the ℓ th segment:

N ℓ=1 d Π S ℓ , (1 -τ ℓ )Π 0 + τ ℓ Π 1 .
We choose a quadratic distance for d. Hence, we minimize with respect to Π 0 (u, v) and Π 1 (u, v) the following function:

For all u in A k and v in A, we obtain the following estimators by Lagrange minimization:

Π 0 (u, v) = b 1 c 2 (u, v) -b 2 c 1 (u, v) a 2 b 1 -a 1 b 2 Π 1 (u, v) = a 2 c 1 (u, v) -a 1 c 2 (u, v) a 2 b 1 -a 1 b 2
where

a 1 = N ℓ=1 (1 -τ ℓ ), a 2 = N ℓ=1 τ ℓ (1 -τ ℓ ), b 1 = N ℓ=1 τ ℓ , b 2 = N ℓ=1 τ ℓ 2 , c 1 (u, v) = N ℓ=1 Π S ℓ (u, v), c 2 (u, v) = N ℓ=1 τ ℓ Π S ℓ (u, v).
Matrices Π 0 and Π 1 are stochastic. In some cases, for small values of N, it is possible to obtain negative terms in the estimated matrices. This problem is solved by a proportional rescaling of the values. Note that we do not obtain a homogeneous model on the segments and that this assumption is only used to get preliminary estimators Π S ℓ . Size m of the segments is chosen in order to minimize the variance of the estimators. Simulations led us to conclude that the value of m minimizing the variance is √ n, where n is the length of the sequence. Variance of estimators is analytically obtained using expectation and variance of estimators on each segment.

Estimation by a point by point method.

Another way to estimate Π 0 and Π 1 is a least squares method. We minimize a quadratic form of the different parameters which is the sum of "prediction errors". At each position t in the sequence, knowing the k-word u = X t-k . . . X t-1 preceding X t , we want Π t n (u, v) to be the nearest possible to 1 if X t = v or the nearest possible to 0 if X t = v. We minimize the sum of error squares:

1 {X t-k ...X t-1 =u} Π t n (u, v) -1 {X t-k ...Xt=uv} . Let us note 1 u for 1 {X t-k ...X t-1 =u} and 1 uv for 1 {X t-k ...X t-1 =u,Xt=v} .
We choose a quadratic distance and then we minimize the following function:

n t=1 u∈A k v∈A 1 {X t-k ...X t-1 =u} Π t n (u, v) -1 {Xt=v} 2 .
For all u in A k and v in A, we obtain the following estimators by Lagrange minimization

Π 0 (u, v) = B 2 (u)C 1 (u, v) -B 1 (u)C 2 (u, v) A 1 (u)B 2 (u) -A 2 (u)B 1 (u) Π 1 (u, v) = A 1 (u)C 2 (u, v) -A 2 (u)C 1 (u, v) A 1 (u)B 2 (u) -A 2 (u)B 1 (u) with A 1 (u) = 2 n t=1 1 u 1 - t n 2 , A 2 (u) = 2 n t=1 1 u 1 - t n t n , B 1 (u) = 2 n t=1 1 u 1 - t n t n , B 2 (u) = 2 n t=1 1 u t n 2 , C 1 (u, v) = 2 n t=1 1 uv 1 - t n , C 2 (u, v) = 2 n t=1 1 uv t n .
Once again, matrices Π 0 and Π 1 are stochastic except in rare cases where negative terms appear. But they are then modified by a proportional rescalement.

Drifting Markov models: polynomial drift

Up to now, we have only considered a linear variation of the transition matrix (DMM of degree 1), but we can generalize to DMMs of higher degree. Thus DMMs have two order parameters: the order k of the Markov model and the degree d of the polynomial drift. To describe such a polynomial model of degree d, we need d + 1 points of support. For linear drift (d = 1), the model was based on the only two matrices of parameters Π 0 and Π 1 . Now, we base our model on d + 1 matrices

Π i d , for 0 ≤ i ≤ d. We choose Π i d
uniformly spaced along the sequence. Any other choice would not be penalizing. Indeed, simulations show that the obtained transition matrices Π t n are similar. The drifting transition matrix has the following form

Π t n (u, v) = d i=0 p i (t)Π i d (u, v),
where p i are the polynomial functions of degree d such that

∀(i, j) ∈ {0, . . . , d} 2 , p i nj d = 1 {i=j} .
Polynomials p i are chosen to have stochastic matrices

Π i d . Hence, for t = ni/d, we have Π t n = Π i d
and for all integer 0 ≤ t ≤ n, we have v∈A Π t n (u, v) = 1.

For d = 1, we intuitively obtain p 0 (t) = 1 -t/n and p 1 (t) = t/n in order to have

Π t n = (1 -t/n) Π 0 + (t/n)Π 1 .
We give their expression for degree d = 2 to illustrate that polynomial functions p i have not a so simple expression than for degree 1. Indeed,

Π t n = p 0 (t)Π 0 + p 1 (t)Π1 2 + p 2 (t)Π 1 leads to Π t n = 2 t 2 n 2 -3 t n + 1 Π 0 + -4 t 2 n 2 + 4 t n Π1 2 + 2 t 2 n 2 - t n Π 1 .
Note that such a system is easy to solve for any degree because it is a simple linear system of (d + 1)(d + 1) independent equations with (d + 1)(d + 1) variables. Nonetheless, we cannot give a general explicit expression for p i with any degree d.

At degree 3, we have

Π t n = - 9 2 t 3 n 3 + 9 t 2 n 2 - 11 2 t n + 1 Π 0 + 27 2 t 3 n 3 - 45 2 t 2 n 2 + 9 t n Π1 3 + - 27 2 t 3 n 3 + 18 t 2 n 2 - 9 2 t n Π2 3 + 9 2 t 3 n 3 - 9 2 t 2 n 2 + t n Π 1 .

Estimation by a matrix regression method.

As for the linear drift, we minimize the following function

N ℓ=1 u∈A k v∈A Π S ℓ (u, v) - d i=0 p i (nτ ℓ )Π i d (u, v) 2 .
Hence, for each (u, v) in A k × A, solving system AX = B where A, X and B are defined below, gives us Once again, as in the linear drift case, we need to minimize the following function

Π i d , estimators of Π i d . • A ij = N ℓ=1 p i (nτ ℓ )p j (nτ ℓ ), 0 ≤ i, j ≤ d; • X i = Π i d (u, v), 1 ≤ i ≤ d; • B i = N ℓ=1 p i (nτ ℓ ) Π S ℓ (u, v), 1 ≤ i ≤ d.
n t=1 u∈A k v∈A 1 {X t-k ...X t-1 =u} Π t n (u, v) -1 {Xt=v} 2 . ( 1
)
As in the matrix regression method, we need to solve for each (u, v) in A k × A, a system AX = B where

• A ij = n t=k 1 {X t-k ...X t-1 =u} p i (t)p j (t), 0 ≤ i, j ≤ d; • X i = Π i d (u, v), 1 ≤ i ≤ d; • B i = n t=k p i (t)1 {X t-k ...X t-1 =u,Xt=v} , 1 ≤ i ≤ d. Hence, we obtain Π i d , estimators of Π i d .

Comparison of the methods

There are some differences between matrix regression method and point by point method. Matrix regression method uses preliminary estimations on each segment S ℓ and the global estimators are computed at a unique point of each segment (the center τ ℓ ). Point by point method enables a direct estimation on all the points of the sequence.

We use the log-likelihood to compare the two methods of estimation of drifting Markov models (see Table 1). We estimate models on the phage Lambda complete genome (see [START_REF] Wu | Nucleotide sequence analysis of DNA. II. Complete nucleotide sequence of the cohesive ends of bacteriophage lambda DNA[END_REF]) and we consider these models as the true models. Then, we simulate a sequence with each one of these models and we compute the log-likelihood for the two estimation methods. We notice that whatever the order, point by point method always gives better likelihood than the regression method. That is a little bit more apparent when we compute the log-likelihood on the real phage Lambda complete genome (see Table 2). Thus, point by point method presents the advantages of a more practical implementation and better results and it is the method which we will privilege thereafter.

Consistence of the estimators

We assume that our estimators are asymptotically unbiased and their variances converge to zero (these theoretical results are not presented here). Thus, our estimators are consistent. In order to show this consistence, we simulate some data where the true model is known. Firstly, we estimate a model on the phage Lambda complete genome and we consider this model as the true model. Then, we simulate some sequences with this model and estimate a mean model on all these sequences. In Appendix B, we give an example with both true and estimated models. For each parameter of the model, we note the absolute value of the difference between the true parameter and the estimated one. In Table 3, we give the mean of these differences for some different drifting Markov models to assume consistence.

Marginal distributions

DMMs offer models which describe faithfully real sequences. This fact is particularly highlighted by the study of the probability distributions of nucleotides in the present section. Indeed, analyzing µ t , the probability distribution of nucleotides at position t associated with our models, is the main way to evaluate their quality. At order 1, the distribution µ t is recursively defined as follows:

µ t+1 (v) = u∈A µ t (u)Π t n (u, v) ∀v ∈ A.
There are similar definitions for order greater than 1. We recall that an ergodic Markov chain Π on a finite state space has a unique stationary probability distribution ν (such that νΠ = ν). The transition matrix Π 0 is ergodic thus, we choose µ 0 as the stationary probability distribution ν 0 of Π 0 . We compute the probability distribution µ t for each position t to analyze the composition of phage Lambda complete genome. First, we draw the evolution, with respect to t, of these distributions of a, c, g and t to observe the differences of composition in the sequence.

We present in Figures 1(a), 1(b), 1(c) and 1(d) these distributions for a modeling of the phage Lambda sequence by a DMM respectively of degree 2, 4, 6 and 8. For degree d = 2, we already notice that gc-rate decreases with respect to the position in the sequence. For degree d = 8, we observe the first gc-rich segment already obtained by an HMM algorithm developed by [START_REF] Muri | Comparaisons d'algorithmes d'identification de chaînes de Markov cachées et application à la détection de régions homogènes dans les séquences d'ADN[END_REF].

Moreover, comparing the HMM segmentation in [START_REF] Muri | Comparaisons d'algorithmes d'identification de chaînes de Markov cachées et application à la détection de régions homogènes dans les séquences d'ADN[END_REF] and the DMM evolution of distributions, we observe similarities. Looking at Figure 1(d), we observe very reliable probability distributions of letters which correspond to the HMM segmentation (see Figure 2). This comparison is interesting because it shows the limits of HMM. Although the first long gc-rich segment is well known and provided by HMM, other parts of the HMM segmentation are not really convincing in view of the evolution of transition probability. Moreover, DMMs are more tractable numerically and they provide a soft evolution contrary to the sudden segmentation of HMMs.

To establish reliability of drifting Markov models, we compare evolutions of probability distributions with nucleotide frequencies. To compute these frequencies, we use sliding windows of size 2000 nucleotides. Figures 3(a), 3(b), 3(c) and 3(d) show that probability distributions of nucleotides of our degree 8 models are very close to the real distribution of nucleotides in the sequence (respectively nucleotide a, c, g and t). Degree 8 is sufficient to observe a good similarity between the curves. In order to compare our polynomial DMM to other Markov models, we draw in Figures 3(a), 3(b), 3(c) and 3(d) the evolutions of probability distributions under an order 1 DMM of degree 0 (it corresponds to a classical homogeneous order 1 Markov model). It turns out that the distance between the two curves is smaller in the case of degree 8 DMM. In the HMM case, we do not observe only one constant probability for each letter as in the Markov model, but few regions with constant probability corresponding to the HMM segmentation. From Figure 4(a) to Figure 4(h), we compare evolutions of probability distributions for an order 1 HMM with 3 hidden states and an order 1 degree 3 DMM on the Phage T4 complete genome (see [START_REF] Miller | Bacteriophage T4 genome[END_REF]). We also compute a distance d df between these evolutions and nucleotide frequencies:

d df = v∈A t (f t (v) -µ t (v)) 2 ,
where f t (v) is the frequency of the nucleotide v at position t and µ t (v) the probability distribution of v at position t. In the way to avoid long computations, we do not take into account all positions t. Then we have no more than 10000 positions, uniformly distributed. It is sufficient in order to compute d df . An order 1 HMM with 3 hidden states and an order 1 DMM with degree 3 have approximately the same number of parameters (in fact, this number is 42 for the HMM and 48 for the DMM). However, we already note that d df is lightly smaller for the DMM: 5.865 versus 5.873. Obviously, this distance is still smaller for a degree 8 DMM (d df = 3.391). In that sense, we show that DMMs represent a new class of flexible models for DNA sequences that can be hoped to provide better fits than HMMs in many cases.

In order to illustrate this fact in another way, we draw in Figure 5, the frequency of gc in the complete genome of phage Lambda. As we said in the introduction, biologists are very concerned in the gc-percent because it may induce presence of genes. They consider five families of isochores: two gc-poor families (L1 and L2) and three gc-rich families (H1, H2 and H3) [START_REF] Bernardi | The vertebrate Genome: Isochores and Evolution[END_REF][START_REF] Oliver | Isochore chromosome maps of eukaryotic genomes[END_REF]. But the transition between two families is often judged to be too sudden when modeled by HMMs. DMM, with its continuous evolution, is a good way to model these transitions. For instance, in Figure 5, between the position 26000 and 32000, we observe a linear increase of the gc content that we model with a degree 1 DMM. Thus DMMs are useful for modeling of heterogeneous phenomena, in particular the linear evolution of gc content, whereas HMM would predict a constant evolution or an abrupt change.

AIC and BIC: comparisons between different models.

In order to analyze drifting Markov models, we compute AIC and BIC values of these models. First of all, we recall the definition of AIC and BIC values (Akaike Information Criterion and Bayesian Information Criterion, introduced respectively by [START_REF] Akaike | A new look at the statistical identification model[END_REF] and [START_REF] Schwarz | Estimating the dimension of a model[END_REF]): where L(θ) is the log-likelihood of the model, K the number of parameters and n the sample size. The model which has the smallest AIC or BIC is considered as the "best" model according to this criterion.

AIC = -(2L(θ) -2K) BIC = -(2L(θ) -K log n)
AIC and BIC are usually built by adding a penalization to the log-likelihood evaluated at the maximum of likelihood. Although least squares estimator is the only disposal, we use here adapted criteria obtained by a penalization of the loglikelihood taken in this estimator. This can be justified by the fact that for Markov chain the mean-square estimation (based on a formula similar to our formula (1)) is asymptotically equivalent to the maximum of likelihood one.

In Table 4, we compute BIC for order 0, 1, 2 and 3 a DMM of degrees from 0 to 5, estimated by the point by point method. These results have been obtained on the Haemophilus influenzae complete genome (see [START_REF] Fleischmann | Whole-genome random sequencing and assembly of haemophilus influenzae rd[END_REF]). Whereas AIC prefers models with a lot of parameters (results are not presented here, but AIC-values generally decrease with order and degree), BIC prefers models with a small number of parameters. That is why a DMMs with high degrees are partially ignored by BIC. Indeed, for an order k DMM of degree d, the number of parameters K is equal to (d + 1)|A| k (|A| -1). You can choose to select order and degree of a DMM with BIC, but higher-dimensional parameterization of DMM provides better fits to the real sequence, as you can see in the precedent section.

Moreover, in order to compare DMMs with other currents methods, the main way is to see in Figures 3(a), 3(b), 3(c) and 3(d) that variations of nucleotides are continuous whatever the position. According to BIC, DMMs are better models than classical Markov models (whatever the order and the degree), but HMMs are better models than DMMs. As DMM can be adjusted to a whole sequence, DMM can be seen as competitive model to HMM. However note that a perspective of this work is to introduce DMMs in HMMs: a DMM could be a hidden state of an HMM. Even if this difference of BIC is not very large, the essential thing to remember is that we provide the first models including the possibility of a continuous variation of the transition matrix. Combined to the quality of HMM, DMM provides powerful tools for sequences analysis.

Replication origin

An example of application of our new models is the search of replication origins on the bacteria. This application draws its inspiration from the program ORILOC [START_REF] Lobry | Oriloc: prediction of replication boundaries in unannotated bacterial chromosomes[END_REF], which has been developed for the prediction of bacterial replications origins. DNA replication is the process of copying a double-stranded DNA strand in a cell, prior to cell division. The two resulting double strands are identical, and each of them consists of one original and one newly synthesized strand. The replication origin is a unique DNA sequence at which DNA replication is initiated. DNA replication may proceed from this point bidirectionally or unidirectionally. Based on the compositional asymmetries between the leading and the lagging strand of replication, the program performs a DNA walk (see [START_REF] Lobry | Genomic landscapes[END_REF]) to obtain the position of the replication origin. A curve is drawn by this program and the peak of this curve corresponds to the replication origin. The values allowing to draw the curve are computed as follows. The first value is 0, and during the walk along the DNA sequence, ORILOC adds 1 each time letter g is found and subtracts 1 each time letter c is found. Thus ORILOC does not rely on a probabilistic model, it draws a curve by running along the real sequence.

We use the same properties of asymmetries in bacterial genomes to perform a detection of the replication origins based on DMMs. Indeed, thanks to the computation of probability distributions of nucleotides at each position t in the sequence, we draw a curve similar to ORILOC. The values of our curve are computed as follows. The first value is 0, and at each position in the sequence, we add the probability of letter g and subtract the probability of letter c.

This work was done on the complete genome of Chlamydia trachomatis (see [START_REF] Stephens | Genome sequence of an obligate intracellular pathogen of humans: Chlamydia trachomatis[END_REF]). Note on Figure 6 the great similarity between the curve obtained by the software ORILOC and the one obtained by DMMs. Note also that our curve is softer than the one of ORILOC because the aim of DMM is to model soft transitions. Although search of replication origins is a break-point detection problem, our method works in the sense that it offers to biologists a window which permits to find the replication origin "in vivo". Then soft transitions do not prevent us to locate the origin of replication. The advantage of our method is to be able to compute analytically a maximum.

Rare words

A second and important example of application of DMMs is the search of rare words in DNA sequences. Many DNA sequence analysis are based on the distribution of the occurrences of patterns having some special biological function. An important problem is to determine the statistical significance of a word frequency in a DNA sequence. [START_REF] Nicodème | Proteome analysis based on motif statistics[END_REF] discuss this relevance of finding over-or under-represented words. The naive idea is the following: a word may have a significant low frequency in a DNA sequence because it disrupts replication or gene expression, whereas a significantly frequent word may have a fundamental activity with regard to genome stability. Well-known examples of words with exceptional frequencies in DNA sequences are biological palindromes corresponding to restriction sites avoided for instance in E. coli [START_REF] Karlin | Statistical analyses of counts and distributions of restriction sites in dna sequences[END_REF]), the Cross-over Hotspot Instigator sites in several bacteria, in E. coli for example [START_REF] Smith | Structure of chi hotspots of generalized recombination[END_REF][START_REF] El Karoui | Characteristics of Chi distribution on different bacterial genomes[END_REF], and uptake sequences [START_REF] Smith | DNA uptake signal sequences in naturally transformable bacteria[END_REF] or polyadenylation signals [START_REF] Van Helden | Statistical analysis of yeast genomic downstream sequences reveals putative polyadenylation signals[END_REF]. The most popular approach consist in fitting a Markov model on the sequence and computing the p-value which is P(N > N obs ) for an over-represented word or P(N < N obs ) for an under-represented word, where N is the random variable of the number of occurrences of the studied word and N obs the number of observed occurrences. We define the pattern statistic associated to any number N obs by:

S = -log 10 P(N > N obs ) if N ≥ E(N) + log 10 P(N < N obs ) if N < E(N) .
This way, a pattern has a positive statistic if it is seen more than expected, a negative statistic if seen less than expected, in both cases, the corresponding p-value is given (in log scale) by the magnitude of the statistic. See [START_REF] Nuel | Numerical Solutions for Patterns Statistics on Markov Chains[END_REF] for a review of the methods available to compute pattern statistics on text generated by a Markov source.

As these probabilities are computed under a model, small p-value can be provided for some words without biological interest if the model is not reliable. That is why it is preferable to rely on a background model the most possible close to the real sequence. DMMs provide it. It always will be more convincing to obtain pvalues for the most realistic models. In that way, considering a DMM for searching rare words in sequences seems to be a better approach than using Markov models (see on Figures 3(a), 3(b), 3(c) and 3(d) that DMM offers a model closer to the reality than Markov model).

Numerical complexities appear when we want to compute exact p-value of inhomogeneous Markov models but a new approach proposed by [START_REF] Nuel | Effective p-value computations using Finite Markov Chain Imbedding (FMCI): application to local score and to pattern statistics[END_REF], using finite Markov chain imbedding (FMCI, see [START_REF] Lou | On runs and longest run tests: A method of finite markov chain imbedding[END_REF]), provides solutions to this problem. We refer to [START_REF] Nuel | Analyse statistique des séquences biologiques: modélisation markovienne, alignements et motifs[END_REF] for a detailed description of this method. We just give here one example of search of rare words. We choose the most popular word in this domain, the Chi site of Escherichia coli K12 (see [START_REF] Blattner | The complete genome sequence of escherichia coli k-12[END_REF]). We consider the complete sequence of the bacteria where the Chi site gctggtgg appears 499 times. As can be seen in Table 6, the Chi site was expected to appear 70.10 times by an order 1 DMM of degree 0 and 175.31 times by an order 2 DMM of order 8. In a more realistic model such as DMM, Chi sites are more expected than in other models. As already said, we cannot compare p-value of different models between them. But we could compare the different classification provided by the different models. Which classification do you prefer? That one given by HMM and its segmentation or that one given by DMM and its soft evolution? Obviously, it is more reliable to consider p-values in the model which provide a better fit to the data even if it is higher-dimensional parameterized. Thus, polynomial DMMs are very useful for the search of rare words in DNA sequences. In Table 5, we give classification of words of size 5, for classical Markov model, 3-states HMM and degree 1 DMM, at order 1.

Discussion and conclusion

We introduce a new class of inhomogeneous Markov models, the drifting Markov models. These new models allow the transition matrix to vary along the sequence.

Notwithstanding the fact that classical Markov models are homogeneous, hidden Markov models cannot model every heterogeneity structures. Heterogeneity of sequences encourages us to consider more flexible models such as drifting Markov models and the continuous variation of their transition matrix. An important illustration of these models concerns the gc-content of a DNA sequence. It is commonly accepted that a high gc-content may induce presence of genes [START_REF] Zoubak | The gene distribution of the human genome[END_REF].

Since they provide a soft evolution and a different transition matrix at each position in the sequence, DMM provides a better fit to the gc-content than HMM with its sudden changes of state. Other applications such as the search of replication origins and especially the search of rare words are very relevant examples of the possibilities of DMMs. We conclude that DMMs are convenient tools for the statistical analysis of sequences. They provide detailed description of the sequence and can be used for structural analysis or direct biological applications. Moreover, it would be interesting not to limit our studies to polynomial drift. Future prospects are to fit new models with co-variables such as the gc-content, the degree of hydrophobicity or an indicator of the protein structure (α-helix, β-sheet. . . ).

Figure 2 :

 2 Figure2: Probability distributions of nucleotides for a DMM of degree 8 in Phage Lambda genome and HMM segmentation with three hidden states (0, 1 and 2, marked on the vertical right-hand side axis).

  17

Vergne

Figure 3 :

 3 Figure 3: Frequency f and probability distribution µ of a for degrees d = 0 and d = 8 DMM in Phage Lambda genome.

Figure 3 :

 3 Figure 3: Frequency f and probability distribution µ of c for degrees d = 0 and d = 8 DMM in Phage Lambda genome.

Figure 3 :

 3 Figure 3: Frequency f and probability distribution µ of g for degrees d = 0 and d = 8 DMM in Phage Lambda genome.

Figure 3 :Figure 4 :Figure 4 :Figure 4 :Figure 4 :

 34444 Figure 3: Frequency f and probability distribution µ of t for degrees d = 0 and d = 8 DMM in Phage Lambda genome.

Figure 4 :

 4 Figure 4: Frequency f and probability distribution µ of t for a degree 3 DMM in Phage T4 genome.

30Figure 5 :

 5 Figure 5: Frequency f and probability distribution µ of gc for a degree 1 DMM between the position 26000 and 32000 in Phage Lambda genome.

Figure 6 :

 6 Figure 6: Search of replication origin in C. trachomatis. Comparison between ORILOC (Lobry 2000) and Drifting Markov Models (DMM).

Table 1 :

 1 Log-likelihood of drifting Markov models computed on sequence simulated by each one of these models (R means regression method and P means point by point method).

	Degree	0	1	2	3	4	5
	Order R -67191 -66999 -66962 -66910 -66909 -66907 0 P -67191 -66999 -66962 -66910 -66909 -66907 Order R -66718 -66504 -66448 -66382 -66376 -66368 1 P -66710 -66501 -66445 -66380 -66374 -66366 Order R -66706 -66482 -66407 -66321 -66295 -66275 2 P -66693 -66477 -66402 -66317 -66290 -66270 Order R -66630 -66331 -66186 -66038 -65938 -65883 3 P -66612 -66320 -66169 -66014 -65898 -65817

Table 2 :

 2 Log-likelihood of drifting Markov models on phage Lambda (R means regression method and P means point by point method).

	Degree	0	1	2	3	4	5
	Order R -67191 -66973 -66934 -66873 -66760 -66680 0 P -67191 -66973 -66934 -66873 66760 -66680 Order R -66743 -66500 -66439 -66362 -66234 -66146 1 P -66714 -66483 -66419 -66345 -66220 -66135 Order R -66052 -65657 -65577 -65438 -65281 -65160 2 P -66005 -65631 -65544 -65410 -65255 -65139 Order R -65661 -65168 -65033 -64809 -64597 -64432 3 P -65579 -65116 -64951 -64746 -64497 -64329

Table 3 :

 3 Comparison between true models and estimated ones. We give the mean of absolute values of differences between true parameters and estimated ones. The number of simulated sequences is given by N.

	Degree	0	2	4	6
	N
	1 0.0026691 0.00943604 0.00920864 0.00996372
	Order 1 10 0.0018065 0.00323508 0.00381477 0.00338738
	100 0.0001906 0.00081519 0.00092941 0.00114088
	1 0.0356951 0.0499792 0.0484183 0.0549545
	Order 2 10 0.0346129 0.0447842 0.0421838 0.0500320
	100 0.0333321 0.0442818 0.0406379 0.0482310
	3 Implementation and results

We developed a program, called DRIMM (as drifting Markov model), dedicated to the estimation of drifting Markov models. This software is written in ANSI C++ and developed on x86 GNU/Linux systems with GCC 3.4, and successfully tested with GCC latest versions on Sun and Apple Mac OSX systems. It relies on seq++ library

[START_REF] Miele | seq++ : analyzing biological sequences with a range of Markov-related models[END_REF]

) and will soon be integrated on seq++ library. Compilation and installation are compliant with the GNU standard procedure. It is available on request from the author. The software is licensed under the GNU General Public License (http://www.gnu.org).

  Figure 1: Probability distributions µ of the 4 nucleotides a, c, g and t for degree 2 DMM in Phage Lambda genome.Figure 1: Probability distributions µ of the 4 nucleotides a, c, g and t for degree 4 DMM in Phage Lambda genome.Figure 1: Probability distributions µ of the 4 nucleotides a, c, g and t for degree 6 DMM in Phage Lambda genome.Figure 1: Probability distributions µ of the 4 nucleotides a, c, g and t for degree 8 DMM in Phage Lambda genome.

	Submission to Statistical Applications in Genetics and Molecular Biology Vergne: Drifting Markov Models Vergne: Drifting Markov Models
		0.34 0.36 0.36 0.38 0.34 0.36 0.36 0.4 0.38 0.35						HMM segmentation µ(a) µ(c) µ(g) µ(t)	µ(a) µ(c) µ(g) µ(t) µ(a) µ(c) µ(g) µ(t) µ(a) µ(c) µ(g) µ(t) µ(a) µ(c) µ(g) µ(t)	2
		0.32 0.34 0.32 0.34
		0.3 0.32 0.3 0.32 0.3
	Distributions of a,c,g,t Distributions of a,c,g,t Distributions of a,c,g,t Distributions of a,c,g,t Distributions	0.24 0.26 0.28 0.26 0.28 0.3 0.24 0.26 0.28 0.24 0.26 0.28 0.3 0.2 0.25								1	States
		0.22 0.24 0.22 0.22
		0.2 0.22 0.2 0.2 0.15
		0.18 0.2 0.18 0.18
		0.1								0
		0.16 0.18 0.16 0.16	0	5000	10000 15000 20000 25000 30000 35000 40000 45000 50000
			0 0 0 0	5000 5000 5000 5000	10000 10000 10000 10000	15000 15000 15000 15000	20000 20000 20000 20000 Position in the sequence 25000 30000 25000 30000 25000 30000 25000 30000	35000 35000 35000 35000	40000 40000 40000 40000	45000 45000 45000 45000	50000 50000 50000 50000
							Position in the sequence (a) Degree 2 Position in the sequence (b) Degree 4 Position in the sequence (c) Degree 6 Position in the sequence (d) Degree 8
	Published by The Berkeley Electronic Press, 2008 Published by The Berkeley Electronic Press, 2008					13 15
	http://www.bepress.com/sagmb					12

14 Submission to Statistical Applications in Genetics and Molecular Biology http://www.bepress.com/sagmb 16 Submission to Statistical Applications in Genetics and Molecular Biology http://www.bepress.com/sagmb

  Frequency f and probability distribution µ of g for a 3-states HMM in Phage T4 genome. Frequency f and probability distribution µ of g for a degree 3 DMM in Phage T4 genome. Frequency f and probability distribution µ of t for a 3-states HMM in Phage T4 genome.

			Vergne: Drifting Markov Models Submission to Statistical Applications in Genetics and Molecular Biology Vergne: Drifting Markov Models Submission to Statistical Applications in Genetics and Molecular Biology
		0.22 0.22 0.4 0.38								f(g) f(g) f(t) f(t)
										µ(g) µ(g) µ(t) µ(t)
		0.21 0.21 0.37
		0.38
		0.2 0.2 0.36
		0.35
		0.19 0.19 0.36
	Distribution of g Distribution of g Distribution of t Distribution of t	0.17 0.18 0.17 0.18 0.34 0.32 0.33 0.34
		0.16 0.16 0.32
		0.31
		0.15 0.15 0.3
		0.3
		0.14 0.14 0.29
		0.13 0.13 0.28 0.28
		0 0 0 0	20000 20000 20000 20000	40000 40000 40000 40000	60000 60000 60000 60000	80000 80000 80000 80000	100000 100000 100000 100000	120000 120000 120000 120000	140000 140000 140000 140000	160000 160000 160000 160000	180000 180000 180000 180000
						Position in the sequence Position in the sequence Position in the sequence Position in the sequence
					(e) Nucleotide g for HMM (f) Nucleotide g for DMM (g) Nucleotide t for HMM
	25 Figure 4: 26 Published by The Berkeley Electronic Press, 2008 http://www.bepress.com/sagmb Figure 4: 27 Published by The Berkeley Electronic Press, 2008 Figure 4: 28 http://www.bepress.com/sagmb

(h) Nucleotide t for DMM

Table 4 :

 4 BIC of drifting Markov models on Haemophilus influenzae.

	Degree	0	1	2	3	4	5
	Order 0	4970473 4970494 4969534 4969471 4969472 4969358
	Order 1	4907845 4907947 4907011	4907051 4907108 4907117
	Order 2	4892907 4893442 4892907 4890996 4893807 4894224
	Order 3 4868040 4870422 4871721 4874079 4876395 4878650

Table 5 :

 5 Classification of words of size 5 of Phage lambda complete genome, for different models according to their pattern statistic S. N obs is the observed number of occurrences of the word. Exp means Expected value. We only give the five first under-represented words and the five first over-represented words.

		MM	HMM 3 states	DMM degree 1
	Words N obs	Exp	S Words N obs	Exp	S Words N obs	Exp	S
	aattg ttggg ttgga acttg taggg gccgg 114 53.97 32 88.22 -11.41 20 65.12 -10.33 21 66.70 -10.29 13 50.74 -9.59 3 29.60 -9.21 tcgag aattg acttg tctag ttgga 12.13 gctgg 127 65.44 32 83.38 -10.07 13 47.59 -8.57 2 24.60 -8.19 21 59.47 -8.15 9 39.01 -8.11 tcgag aattg ttgga ttggg acttg 14.23 gctgg 127 64.80 32 86.53 -10.94 21 64.94 -9.76 20 62.94 -9.66 13 50.27 -9.44 9 40.69 -8.68 11.77
	ctgaa 124 61.02	12.16 ctgaa 124 61.34	14.90 ctgaa 124 60.85	12.21
	tccgg 100 39.98	15.08 ccgga 112 44.00	20.58 tccgg 100 38.81	16.18
	ccgga 112 43.11	17.93 tccgg 100 36.50	20.65 ccgga 112 43.57	18.10
	gcaga 141 57.51	20.20 gcaga 141 58.35	22.66 gcaga 141 57.59	20.31

Table 6 :

 6 Pattern statistic S (log p-value) of the over-represented word gctggtgg for DMMs of different orders and degrees: the Chi of E. coli which appears 499 times in the sequence. Note that a DMM of degree 0 corresponds to a classical Markov model.

	Order Degree Expected value	S
	1	0	70.10 240.814
	1	1	70.26 240.398
	1	2	71.88 238.766
	1	3	71.87 238.774
	1	8	71.94 238.605
	2	0	173.84 88.902
	2	1	174.03 88.747
	2	2	175.16 87.837
	2	3	175.10 87.881
	2	8	175.31 87.717

Submission to Statistical Applications in Genetics and Molecular Biologyhttp://www.bepress.com/sagmb

Vergne: Drifting Markov ModelsPublished by The Berkeley ElectronicPress, 2008

http://www.bepress.com/sagmb

Vergne: Drifting Markov ModelsPublished by The Berkeley ElectronicPress, 2008

Vergne: Drifting Markov ModelsPublished by The Berkeley ElectronicPress, 2008

Vergne: Drifting Markov ModelsPublished by The Berkeley ElectronicPress, 2008

Submission to Statistical Applications in Genetics and Molecular Biologyhttp://www.bepress.com/sagmb

We give here an example of the systems we would need to solve to provide estimation by maximum likelihood. It corresponds to one of the 4 systems of 6 equations with 6 variables obtained for an order 1 DMM of degree 1 for the nucleotide alphabet A = {a, c, g, t}. Knowing that n is very high (as the length of the DNA sequence) and that all the parameters Π 0 (u, v) and Π 1 (u, v) are in the denominator of each equality, you can note the numerical complexity which precludes the use of this natural method. Obviously, complexity still is a problem for DMM of higher order and degree.

Submission to Statistical Applications in Genetics and Molecular Biology

Appendix B: Consistence of the estimators

In order to show consistence of our estimators, we give here one example of the results of some simulations. In this example, we estimate an order 1 degree 2 drifting Markov model on the phage Lambda complete genome (see [START_REF] Wu | Nucleotide sequence analysis of DNA. II. Complete nucleotide sequence of the cohesive ends of bacteriophage lambda DNA[END_REF]) and we consider this model as the true model. Then, we simulate 10 sequences with this model and estimate a mean model on all these sequences. Both true and estimated models are shown in


Olivier David

Catherine Larédo

Raphael Leblois

Brigitte Schaeffer

Nicolas Vergne

Coalescent-based DNA barcoding: multilocus analysis and robustness

Keywords: Bayesian inference, classification, coalescent, DNA barcoding, species assignment

published or not. The documents may come

Introduction

DNA barcoding is the assignment of individuals to species or higher taxonomic levels using standardized genetic data observed on the target individuals and samples from each species [START_REF] Frézal | Four years of DNA barcoding: Current advances and prospects[END_REF][START_REF] Valentini | DNA barcoding for ecologists[END_REF]. The DNA barcode project is conceived as a standard system for fast and accurate identification of all eukaryotic species [START_REF] Hebert | Barcoding animal life: cytochrome c oxidase subunit 1 divergences among closely related species[END_REF][START_REF] Miller | DNA barcoding and the renaissance of taxonomy[END_REF]. The DNA barcode itself consists of a 648 bp region of the cytochrome c oxidase 1 (COI) gene. Additionally to the mitochondrial COI gene, nuclear loci are sometimes also considered to improve assignment performance [START_REF] Austerlitz | DNA barcode analysis: a comparison of phylogenetic and statistical classification methods[END_REF][START_REF] Elias | Limited performance of DNA barcoding in a diverse community of tropical butterflies[END_REF].

DNA barcoding is a classification problem rather than a clustering one since the classes (species) are predefined and do not have to be inferred from the data (but see [START_REF] Pons | Sequence-based species delimitation for the DNA taxonomy of undescribed insects[END_REF] for an application of clustering to barcoding). Barcoding assignment methods can be divided into similarity methods based on the match between the query sequence and the reference sequences such as BLAST search, phylogenetic approaches [START_REF] Hebert | Barcoding animal life: cytochrome c oxidase subunit 1 divergences among closely related species[END_REF][START_REF] Elias | Limited performance of DNA barcoding in a diverse community of tropical butterflies[END_REF], classification algorithms with no underlying biological models such as the nearest-neighbour method and methods based on population genetics [START_REF] Matz | A likelihood ratio test for species membership based on DNA sequence data[END_REF]. Two Bayesian methods based on models have recently been developed. In the method of Munch et al. (2008a), species are assumed to evolve according to a phylogenetic model while the within-species variation is not modelled. Conversely, in TheAssigner, the method of [START_REF] Abdo | A step toward barcoding life: A model-based, decisiontheoretic method to assign genes to preexisting species groups[END_REF], species are assumed to evolve independently and the dependence between sequences within species is modelled using a classical population genetics model called the coalescent. The latter is a model for the genealogical tree of a random sample of genes drawn from a large panmictic population (Chapter 10, [START_REF] Ewens | Mathematical population genetics[END_REF]Kingman, 1982a,b;[START_REF] Tajima | Evolutionary relationship of DNA sequences in finite populations[END_REF]. Model-based barcoding methods raise various issues. Current methods assume that the data are mitochondrial and cannot cope with nuclear data. Moreover, their robustness to departures from model assumptions has been little studied.

The main objective of the present paper is to study how to take account of nuclear data in coalescent-based classification and to study the robustness of this type of classification to departures from model assumptions. First a coalescent-based classification for assigning individuals to species using mitochondrial data is developed (Section 3). Then this method is extended to take account of nuclear data (Section 4). Finally the performance and robustness of coalescent-based classification are studied using simulated and real data sets (Sections 5 and 6).

Bayesian classification

First we briefly review some basic material on Bayesian classification. In this method, individuals are assumed to belong to c classes. A data set y is available that includes measurements observed on reference individuals whose class is known. The objective is to predict the class z ∈ {1, . . . , c} of a test individual given its data x and the reference data y.

In Bayesian classification, a test individual is assigned to the class with the largest posterior probability of membership [START_REF] Abdo | A step toward barcoding life: A model-based, decisiontheoretic method to assign genes to preexisting species groups[END_REF]Munch et al., 2008a;[START_REF] Ripley | Pattern recognition and neural networks[END_REF]. The assignment may be considered as ambiguous if the latter probability does not exceed some specified threshold. According to Bayes theorem, the posterior probability that a test individual belongs to class i is equal to P (z = i|y, x) = P (z = i, x|y)/P (x|y) = r i / k r k , where:

r i = P (z = i|y)P (x|y, z = i).
(1)

In this equation, P (z = i|y) is the probability that the test individual belongs to class i given the reference data y prior to the knowledge of x and plays the role of a prior probability of membership. The probability P (x|y, z = i) is the conditional probability that an individual sampled in class i has data x. Bayesian classification is optimal for the 0-1 loss function (Chapter 2, [START_REF] Ripley | Pattern recognition and neural networks[END_REF] and provides a measure of assignment confidence.

3 Species assignment with mitochondrial data

We now apply Bayesian classification to DNA barcoding. In this section, the data consist of mitochondrial DNA sequences. The assumed demographic model is a set of isolated and panmictic species with a common ancestry at a given time in the past (i.e., the divergence time). This demographic model is the same as the one of [START_REF] Abdo | A step toward barcoding life: A model-based, decisiontheoretic method to assign genes to preexisting species groups[END_REF]. The mitochondrial locus is assumed to evolve according to the coalescent model within each species independently. Following the standard coalescent, it is assumed that species sizes do not vary over time, that there is no migration between species and that all alleles are neutral. All individuals are assumed to be sampled at the same time and the species of any test individual is assumed to be represented in the reference data y. In this model, mutations occur on each ancestral lineage of species i according to a Poisson process with parameter θ i /2. The assumed mutation model is the infinitely many-sites model (ISM), in which a gene is considered as an infinitely long DNA sequence and each new mutant site is sampled uniformly and independently along the sequence (Chapter 9, [START_REF] Ewens | Mathematical population genetics[END_REF].

Finally, it is assumed that at each site it is known which base is the mutant base or the ancestral base (Section 7) and that there are no missing data or errors in the data.

The mutation parameters θ i are first assumed to be known. Then, under the assumption that species evolve independently, the probability P (x|y, z = i) in ( 1) is equal to:

P (x|y, z = i) = P (x|y i , z = i),
where y i denotes the data of species i in the reference data base. This probability will be written for simplicity as P (x|y i ) in what follows. Generally it cannot be calculated explicitly under the ISM but it can be estimated as follows. It is equal to (p. 420, De Iorio

P (x|y i ) = n i (x) + 1 n i + 1 P 0 (x, y i ) P 0 (y i ) , (2)
where P 0 is the probability of an unordered sample, n i is the number of genes in the sample of species i and n i (x) is the number of genes with sequence x in the sample of species i.

The probabilities P 0 (x, y i ) and P 0 (y i ) can be estimated using importance sampling (IS)

(De Iorio and Griffiths, 2004, supplementary materials A). Note that the probability P 0 (y i )

needs to be estimated only once if there are several individuals to assign.

Mutation processes are generally unknown for most species and the vector θ of mutation parameters is thus usually not known. In this case, the posterior probabilities of membership can be estimated by plug-in, that is by assuming that θ is known and equal to an estimate θ, computed from the reference data set as in [START_REF] Abdo | A step toward barcoding life: A model-based, decisiontheoretic method to assign genes to preexisting species groups[END_REF]. The vector θ may be estimated, for example, using the method of [START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF], by coalescent-based maximum likelihood or Bayesian methods [START_REF] Bahlo | Inference from gene trees in a subdivided population[END_REF][START_REF] Kuhner | Estimating effective population size and mutation rate from sequence data using Metropolis-Hastings sampling[END_REF].

Alternatively, a predictive approach can be used in which the dependence of probabilities on θ is removed by integration (Chapter 2, Ripley, 1996, supplementary materials A).

4 Species assignment with mitochondrial and nuclear data Individuals are now assumed to be genotyped at l diploid nuclear loci in addition to the mitochondrial locus. The two sequences of an individual at a nuclear locus are assumed to be known (Section 7). The genetic data of a test individual are denoted by x = (x 0 , . . . , x l ), where x 0 is the mitochondrial sequence and x j (j ≥ 1) is the pair of sequences at nuclear locus j. Each locus is assumed, as in the previous section, to evolve according to the coalescent model within each species independently. All the loci are assumed to evolve independently [START_REF] Hudson | Gene genealogies and the coalescent process[END_REF][START_REF] Nordborg | Coalescent theory[END_REF] and there is no recombination within a locus. Mutations are assumed to occur according to the ISM with parameter θ ij /2 for species i and locus j. Finally, for simplicity, all parameters θ ij are assumed to be known in this section.

With independent loci, the quantity r i in (1) is equal to (Chapter 8, [START_REF] Ripley | Pattern recognition and neural networks[END_REF]:

r i = P (z = i|y) l j=0 P (x j |y ij , z = i), (3)
where y ij is the reference data for species i and locus j. This equation allows us to easily combine the mitochondrial and the nuclear informations. For a nuclear locus in a diploid species, (2) becomes:

P (x j |y ij , z = i) = (n ij (x j1 ) + 1)(n ij (x j2 ) + 1 + δ j ) (n ij + 2)(n ij + 1) P 0 (x j , y ij ) P 0 (y ij ) ,
where δ j = 0 if the test individual is heterozygote at locus j and δ j = 1 if the test individual is homozygote at locus j. In this equation, x j1 and x j2 denote the two test sequences at locus j, n ij denotes the number of genes sampled for species i at locus j, n ij (x j1 ) denotes the multiplicity of allele x j1 in the sample of species i and locus j. The probabilities P 0 (x j , y ij )

Simulation study

Simulations were carried out to assess the methods described above. In these simulations, one ancestral species split T generations ago into two new species with effective size N e and mutation parameter θ. There were n reference individuals in each species. First sequences were simulated for a mitochondrial locus and diploid nuclear locus to study the effect of adding nuclear data. Then to test the robustness of the methods developed, mitochondrial sequences were simulated assuming that species size varied over time or that each species was divided into several populations exchanging migrants. To mimic extreme sampling strategies that can be done in structured populations, we considered an "extended" sampling, in which the reference individuals were sampled in all populations for each new species, and a "clustered" sampling, in which all reference individuals were sampled from a single population in each new species. Details on these simulations are presented in supplementary materials B.

The simulated data were analysed with the nearest-neighbour classification (1NN) and the developed Bayesian assigner (BA) (supplementary materials B). The 1NN method was used because it had been found to be efficient compared with other barcoding methods [START_REF] Austerlitz | DNA barcode analysis: a comparison of phylogenetic and statistical classification methods[END_REF] and it was expected to be robust since it was not based on a specific biological model. This method was implemented with bagging in order to obtain a measure of confidence for an assignment [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF], supplementary materials B).

Assignment performance was quantified using sensitivity and specificity (Munch et al., 2008a,b). Specificity is the fraction of non-ambiguous assignments (Section 2) that are correct. Sensitivity is the fraction of all the assignments that are correct.

The simulations with nuclear data first showed that performances were the best for the combination of the mitochondrial and the nuclear data, intermediate for the mitochondrial data (Fig. 1), and the least good for the nuclear data alone (Fig. S2 andS3). The poor results for the nuclear data alone were probably due to the larger effective size we used for the nuclear locus, leading to smaller scaled divergence times T /N e and thus lower levels of differentiation between the two new species. Nevertheless adding nuclear data clearly increased sensitivity (Fig. 1). This was mainly due to a reduction of the number of ambiguous assignments since specificity did not increase much (Fig. 1). Our simulations also showed that 1NN and BA had similar performances, except for the nuclear data alone for which 1NN had a low sensitivity (Fig. S3). However, we can note that BA had more ambiguous assignments than 1NN but made fewer errors among the non-ambiguous assignments (Fig. 1). Another important result was that the estimation of mutation parameters did not change the BA performance much (Fig. S1). Finally and as it was expected, increasing the values of θ, T or n improved the performance of both methods as in [START_REF] Austerlitz | DNA barcode analysis: a comparison of phylogenetic and statistical classification methods[END_REF].

For past population size variations, the main results were that past expansions strongly increased specificity, sensitivity and the rate of non-ambiguous assignments, whereas past contractions had the opposite effect of decreasing specificity and sensitivity (Fig. 2). Our simulations also showed that past expansions affected both methods similarly, but 1NN always showed a slightly better performance than BA. On the contrary, it is interesting to note that the effect of past contractions was more pronounced for 1NN than for BA, resulting in much better performances for BA. Finally, the effect of past population size variations was found to be important for all the growth rate values we used and to be stronger for expansions than for contractions.

The effect of population structure was more complex because it depended on the sampling strategy. Compared with the unstructured species results, a population structure with a weak migration mainly affected sensitivity and the rate of non-ambiguous assignments, that both increased for the "clustered" samples and decreased for the "extended" samples (Fig. 2). This result was unexpected as the population of origin of a test individual was represented by two individuals in the reference samples for the "extended" samples but not for the "clustered" samples. Finally, we note that population structure affected both methods similarly and that the effect of population structure became noticeable only when migration was weak enough. θ=3, T=500, n=5

Probability threshold Specificity q q q q q q q q q 1NNm BAm 1NNmn BAmn 0.55 0.65 0.75 0.85 0.4 0.6 0.8 θ=3, T=500, n=5

Probability threshold Sensitivity q q q q q q q q 0.55 0.65 0.75 0.85 0.85 0.90 0.95 1.00 θ=20, T=500, n=5

Probability threshold Specificity q q q q q q q q 0.55 0.65 0.75 0.85 0.4 0.6 0.8 θ=20, T=500, n=5

Probability threshold Sensitivity q q q q q q q q Figure 1: Effect of adding nuclear data on the performance of coalescent-based barcoding. Specificity is the fraction of non-ambiguous assignments that are correct. Sensitivity is the fraction of all the assignments that are correct. The probability threshold is the threshold used to decide if an assignment is ambiguous. 1NN and BA are the nearest-neighbour classification and the developed Bayesian assigner with a known value of θ. The subscripts m and mn denote the mitochondrial data and the combination of mitochondrial and nuclear data, respectively. Adding nuclear data increases sensitivity and reduces the ambiguity of assignments. rep(cproba, times = 2) q q q q q q q q q 1NN BA 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.0 0.2 0.4 0.6 0.8 1.0 rep(cproba, times = 2) q q q q q q q q q 1NN BA 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.65 0.75 0.85 0.95 rep(cproba, times = 2) q q q q q q q q 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.0 0.2 0.4 0.6 0.8 1.0 rep(cproba, times = 2) q q q q q q q q 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.65 0.75 0.85 0.95 rep(cproba, times = 2) q q q q q q q q 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.0 0.2 0.4 0.6 0.8 1.0 rep(cproba, times = 2) q q q q q q q q 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.65 0.75 0.85 0.95 q q q q q q q q 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.0 0.2 0.4 0.6 0.8 1.0 rep(cproba, times = 2) q q q q q q q q Probability threshold Specificity Sensitivity θ=3, T=500 Nm=0.01 n=8 "clustered" θ=3, T=500 Nm=0.01 n=8 "extended" θ=3, T=500 Gf=0.01 n=5 θ=3, T=500 Gf=10 n=5

Astraptes

Probability threshold Specificity q q q q q q q q q 1NN TA BA 0.55 0.65 0.75 0.85 0.5 0.6 0.7 0.8 0.9 1.0

Astraptes

Probability threshold Sensitivity q q q q q q q q 0.55 0.65 0.75 0.85 0.75 0.85 0.95

Ithomiinae

Probability threshold Specificity q q q q q q q q q 1NNm TAm BAm BAmn 0.55 0.65 0.75 0.85 0.5 0.6 0.7 0.8 0.9 1.0

Ithomiinae

Probability threshold Sensitivity q q q q q q q q Figure 3: Performance of coalescent-based barcoding with real data. 1NN, TA and BA are the nearest-neighbour classification, TheAssigner and the developed Bayesian assigner. The subscripts m and mn denote the mitochondrial data and the combination of mitochondrial and nuclear data. Adding nuclear data increases sensitivity and reduces the ambiguity of BA assignments. No method has the best specificity for both data sets.

Discussion

Classification inputs. Bayesian classification requires prior probabilities of membership.

When these probabilities are not known, they may be estimated from the reference data provided that these data can be considered as a random sample among all the species considered (page 53, [START_REF] Ripley | Pattern recognition and neural networks[END_REF] or they may be fixed to 1/c.

The developed methods require the ancestral sequence of each locus. If this sequence is not known, it can be inferred from the data [START_REF] Bahlo | Inference from gene trees in a subdivided population[END_REF][START_REF] Gascuel | Inferring ancestral sequences in taxon-rich phylogenies[END_REF] or posterior probabilities of membership may be estimated using unrooted trees (Section 5, [START_REF] Tavaré | Lectures on probability theory and statistics : Ecole d'été de probabilités de Saint-Flour XXXI -2001[END_REF][START_REF] Bahlo | Inference from gene trees in a subdivided population[END_REF]. Moreover, many sequences from the barcoding reference database could be used as outgroups and thus greatly facilitate the inference of the ancestral sequence.

Finally, both alleles of an individual at a nuclear locus were assumed to be known.

Current genotyping technologies are able to determine which two bases are present at each site of a nuclear locus but not the two sequences of the locus. It is a general problem for most nuclear sequence analysis methods, and statistical methods, known as phasing methods, can infer these two sequences from unphased data together with missing data [START_REF] Scheet | A fast and flexible statistical model for large-scale population genotype data: applications to inferring missing genotypes and haplotypic phase[END_REF].

Classification assumptions. The mutation model considered in this paper was the ISM, a model that requires fewer computations than models with a finite number of sites. However it assumes that a particular mutation can only occur once so that in particular there is no homoplasy. It is more adapted to situations where species are closely related since the assumption of absence of homoplasy is more likely to be satisfied in this case. This does not seem to be a problem for DNA barcoding since species that are distantly related to a test individual can be discarded using simpler methods [START_REF] Austerlitz | DNA barcode analysis: a comparison of phylogenetic and statistical classification methods[END_REF]Munch et al., 2008b). In our study, classification methods were compared using data sets that were compatible with the ISM so that all the methods had the same amount of information.

Species classification based on the ISM could be extended to account for different mutation rates for transitions and transversions.

The species of a test individual was assumed to be represented in the reference data.

The conditional probabilities of an allele P (x|y i ) can be used to check if this assumption is satisfied: low probabilities are an indication that this assumption may not be satisfied.

The developed methods are based on various simplifying assumptions. It would be interesting to relax some of these assumptions to improve classification performance. The program genetree can perform likelihood estimations with varying population size and population structures under the ISM [START_REF] Bahlo | Inference from gene trees in a subdivided population[END_REF]. Divergence models and models that combine phylogenetic and population genetics models do not assume that species are independent [START_REF] Matz | A likelihood ratio test for species membership based on DNA sequence data[END_REF][START_REF] Pons | Sequence-based species delimitation for the DNA taxonomy of undescribed insects[END_REF].

Performance of the developed methods. The method developed to combine mitochondrial and nuclear informations appeared satisfactory. Adding nuclear data reduced the ambiguity of assignments in our analyses.

We showed that coalescent-based classification was robust to departures from demographic stability and panmixia and to designs that did not sample the within-species variation efficiently. It performed similarly to a model-free method (1NN) in the robustness study. Demographic expansion was found to increase the power of barcoding. This is an expected result, however, considering that speciation events are probably often associated with founder events followed by demographic expansions or selective sweeps on the mitochondria, it may highlight the reasons why DNA barcoding works so well with a limited sequence information.

Finally, no assignment method was found to be always the best in our analyses. Similar results were obtained by [START_REF] Austerlitz | DNA barcode analysis: a comparison of phylogenetic and statistical classification methods[END_REF] when comparing phylogenetic and statistical methods. However the developed Bayesian assigner generally appeared more cautious than the other methods in the sense that it assigned fewer individuals but made fewer errors

Figure 2 :

 2 Figure2: Robustness of coalescent-based barcoding to past population size changes and population structures. 1NN and BA are the nearest-neighbour classification and the developed Bayesian assigner with estimated mutation parameters. Results for past population size changes are presented on the first two lines, with Gf being the growth factor. A growth factor larger than one indicates a population expansion from divergence to present, whereas a growth factor smaller than one indicates a population decline. Results for population structures are presented on the last two lines, with N m being the number of migrants exchanged between adjacent populations in one generation. BA appears robust since its performance is similar to the one of 1NN that is model-free. 14

among the assigned individuals.

The supplementary materials referenced in Sections 3, 5 and 6 are available at arxiv.org.

Acknowledgements

This study was funded by the Agence Nationale de la Recherche (IFORA ANR-06-BDIV-014 and EMILE NT09-611697 projects). We thank F. Austerlitz for helpful comments.

We chose to test our method on two different data sets that contained both differentiated and undifferentiated species. The first data set used came from the study of [START_REF] Hebert | Ten species in one: Dna barcoding reveals cryptic species in the neotropical skipper butterfly astraptes fulgerator[END_REF] on Astraptes species and consisted of mitochondrial sequences (CO1 locus). The second data set used came from the study of [START_REF] Elias | Limited performance of DNA barcoding in a diverse community of tropical butterflies[END_REF] on Ithomiinae species and consisted of mitochondrial (CO1 locus) and nuclear data (EF1α locus). The data were analysed with 1NN, BA and TheAssigner [START_REF] Abdo | A step toward barcoding life: A model-based, decisiontheoretic method to assign genes to preexisting species groups[END_REF]. The performance of each method was quantified using a leave-one-out analysis in which each haplotype was used as a test sequence after reducing its multiplicity by one in the reference data. Details on these data sets and their analyses are given in supplementary materials C.

The results first showed that adding nuclear data reduced the ambiguity of the BA assignments (Fig. 3). The analyses also showed that no method had the highest specificity in all cases (Fig. 3). Moreover BA had a lower sensitivity than the other methods and thus assigned fewer individuals (Fig. 3), except for the nuclear Ithomiinae data alone (Fig. S4).

Another result of our analyses was that some posterior probabilities of membership were sensitive to the choice of the ancestral bases (supplementary materials C). Finally a few conditional probabilities were estimated with the predictive method (supplementary materials A) and the corresponding estimates were close to the plug-in estimates.

Disclosure statement

No competing financial interests exist.