T O Gallouet

G Mijoule

Y Swan

REGULARITY OF SOLUTIONS OF THE STEIN EQUATION AND RATES IN THE MULTIVARIATE CENTRAL LIMIT THEOREM

Keywords: Berry-esseen bounds, Stein's method, Elliptic regularity, AMS subjects classification

Regularity of solutions of the Stein equation and rates in the multivariate central limit theorem

Thomas Gallouët, Guillaume Mijoule, Yvik

1. Introduction 1.1. Multivariate Stein's method. Stein's method is a powerful tool for estimating distances between probability distributions. It first appeared in [START_REF] Stein | A bound for the error in the normal approximation to the distribution of a sum of dependent random variables[END_REF], where the method was introduced for the purpose of comparison with a (univariate) Gaussian target. The idea, which still provides the backbone for the contemporary instantiations of the method, is as follows. If Z is a standard Gaussian random variable, then [START_REF] Barbour | Stein's method for diffusion approximations[END_REF] E[f ′ (Z) -Zf (Z)] = 0 for all absolutely continuous functions f with E|f ′ (Z)| < +∞. Let X be another random variable, and consider the integral probability distance between the laws of X and Z given by

(2) d H (X, Z) = sup h∈H E[h(X) -h(Z)],
for H a class of tests functions which are integrable with respect to the laws of both X and Z.

Many classical distances admit a representation of the form [START_REF] Bentkus | On the dependence of the Berry-Esseen bound on dimension[END_REF], including the Kolmogorov (with H the characteristic functions of half-lines), total variation (with H the characteristic functions of Borel sets), and 1-Wasserstein a.k.a. Kantorovitch (with H the 1-Lipschitz real functions) distances; see e.g. [START_REF] Müller | Integral probability metrics and their generating classes of functions[END_REF]. Letting ω denote the standard Gaussian pdf, we define for every h ∈ H the function

(3) f h (x) = 1 ω(x) x -∞ (h(y) -Eh(Z)) ω(y)dy.
This function is a solution to the ODE (called a Stein equation) [START_REF] Bobkov | Entropic approach to E. Rio's central limit theorem for w2 transport distance[END_REF] f ′ h (x)xf h (x) = h(x) -Eh(Z), x ∈ R, which allows to rewrite integral probability metrics [START_REF] Bentkus | On the dependence of the Berry-Esseen bound on dimension[END_REF] as

d H (X, Z) = sup h∈H E[f ′ h (X) -Xf h (X)].
Stein's intuition was to exploit this last identity to estimate the distance between the laws of X and Z. One of the reasons for which the method works is the fact that the function f h defined in (3) enjoys many regularity properties. For instance, one can show (see e.g. [17, pp. 65-67]) that if h is absolutely continuous then

(5)

f h ∞ ≤ 2 h ′ ∞ , f ′ h ∞ ≤ 2/π h ′ ∞ and f ′′ h ∞ ≤ 2 h ′ ∞ , .
∞ holding for the supremum norm. This offers a wide variety of handles on d H (X, Z) -typically via low order Taylor expansion arguments -for all important choices of test functions H and under weak assumptions on X. This observation has been used, for instance, to obtain Berry-Esseen-type bounds in the classical central limit theorem in 1-Wasserstein distance, Kolmogorov or total variation distances, see [START_REF] Chen | Normal approximation by Stein's method[END_REF][START_REF] Nourdin | Normal approximations with Malliavin calculus: from Stein's method to universality[END_REF].

Consider now a d-dimensional Gaussian target Z ∼ N (0, I d ). The d-dimensional equivalent to identity [START_REF] Barbour | Stein's method for diffusion approximations[END_REF] was identified in [START_REF] Barbour | Stein's method for diffusion approximations[END_REF][START_REF] Götze | On the rate of convergence in the multivariate CLT[END_REF] as

E[∆f (Z) -Z • ∇f (Z)] = 0,
which holds for a "large class" of functions f : R d → R (x • y denotes the usual scalar product between vectors x, y ∈ R d ). We will define the "large class" of functions precisely in Proposition 2.1 below. For h a function with finite Gaussian mean, the multivariate Stein equation then reads [START_REF] Bonis | Rates in the central limit theorem and diffusion approximation via Stein's method[END_REF] ∆f (x)

-x • ∇f (x) = h(x) -Eh(Z), x ∈ R d .
Note that ( 6) is a second order equation in the unknown function f

; in dimension d = 1, (6) reduces to f ′′ (x) -xf ′ (x) = h(x) -Eh(Z)
which is obtained by applying (4) to f ′ . Barbour [START_REF] Barbour | Stein's method for diffusion approximations[END_REF] identified a solution of (6) to be (7)

f h (x) = - 1 0 1 2t E[h( √ tx + √ 1 -tZ) -h(Z)] dt,
and the same argument as in the 1-dimensional setting leads to the identity

(8) d H (X, Z) = sup h∈H E [∆f h (X) -X • ∇f h (X)] ,
which is the starting point for multivariate Gaussian approximation via Stein's method. The explicit representation [START_REF] Chatterjee | Multivariate normal approximation using exchangeable pairs[END_REF] is suitable to obtain regularity properties of f h in terms of those of h; for instance (see e.g. [START_REF] Reinert | Multivariate normal approximation with Stein's method of exchangeable pairs under a general linearity condition[END_REF]Lemma 2.6]) it is known that if h is n times differentiable then f h is n times differentiable and (9)

∂ k f h (x) k j=1 ∂x i j ≤ 1 k ∂ k h(x) k j=1 ∂x i j ,
for every x ∈ R d . Hence, contrarily to the univariate case where first order assumption on h was sufficient to deduce second order regularity for f h (recall (5)), a bound such as [START_REF] Chernozhukov | Central limit theorems and bootstrap in high dimensions[END_REF] only shows the same regularity for h and f h . In most practical implementations of the method, however, Taylor expansiontype arguments are used to obtain the convergence rates from the rhs of (8); hence regularity of f h is necessary in order for the argument to work. This restricts the choice of class H in which the statements are made and therefore weakens the strength -be it only in terms of the choice of distance -of the resulting statements. An important improvement in this regard is due to Chatterjee and Meckes [START_REF] Chatterjee | Multivariate normal approximation using exchangeable pairs[END_REF] who obtained (among other regularity results) that [START_REF] Courtade | Existence of Stein kernels under a spectral gap, and discrepancy bound[END_REF] sup

x∈R d ∇ 2 f h (x) H.S ∇h ∞ ,
( M H.S. stands for the Hilbert-Schmidt norm of a matrix M and ∇ 2 f h for the Hessian of f h ); see also [START_REF] Reinert | Multivariate normal approximation with Stein's method of exchangeable pairs under a general linearity condition[END_REF]. Gaunt [START_REF] Gaunt | Rates of convergence in normal approximation under moment conditions via new bounds on solutions of the Stein equation[END_REF] later showed a generalization of this result, namely a version of [START_REF] Chernozhukov | Central limit theorems and bootstrap in high dimensions[END_REF] where the derivatives of order k of f h can be bounded by derivatives of order k -1 of h. This still does not concur with the univariate case as we know that, in this case and when h ′ is bounded, one can bound one higher derivative of f h : indeed, it holds |f

h | 2|h ′ | (here the function f h is the solution [START_REF] Chatterjee | Multivariate normal approximation using exchangeable pairs[END_REF] to the univariate version of the second order equation ( 6)). This loss of regularity is, however, not an artefact of the method of proof but is inherent to the method itself: Raič [START_REF] Raič | A multivariate CLT for decomposable random vectors with finite second moments[END_REF] exhibits a counterexample, namely a Lipschitz-continuous function such that the second derivative of f h is not Lipschitz-continuous. We will discuss this example in detail later on.

1.2. Multivariate Berry-Esseen bounds. Let (X i ) i 1 be an i.i.d. sequence of random vectors in R d , and for simplicity take them centered with identity covariance matrix. Let W = n -1/2 n i=1 X i , Z ∼ N (0, I d ) and consider the problem of estimating D(W, Z), some probability distance between the law of Z and that of W . According to [START_REF] Götze | On the rate of convergence in the multivariate CLT[END_REF], the earliest results on this problem in dimension d ≥ 2 concern distances of the form (2) with H indicator functions of measurable convex sets in R d (which is a multivariate generalization of the Kolmogorov distance). The best result for this choice of distance is due to [START_REF] Bentkus | On the dependence of the Berry-Esseen bound on dimension[END_REF] where an estimate of the form d H (W, Z) ≤ 400

d 1/4 n -1/2 E[|X 1 | 3 ] is shown (| • | is the Euclidean norm)
; the dependence on the dimension is explicit and the best available for these moment assumptions and this distance. More recently, a high dimensional version of the same problem was studied in [START_REF] Chernozhukov | Central limit theorems and bootstrap in high dimensions[END_REF], with H the class of indicators of hyper-rectangles in R d ; we also refer to the latter paper for an extensive and up-to-date literature review on such results.

Another important natural family of probability distances are the Wasserstein distances of order p (a.k.a. Mallows distances) defined as [START_REF] Fathi | Stein kernels and moment maps[END_REF] W

p (W, Z) = (inf E [|X 1 -Y 1 | p ]) 1/p
where the infimum is taken over all joint distributions of the random vectors X 1 and Y 1 with respective marginals the laws of W and Z. Except in the case p = 1, such distances cannot be written under the form (2); as previously mentionned, when p = 1 the distance

W := W 1 in (
W p (•, •) ≥ W 1 (•,
•) for all p ≥ 1. We refer to [START_REF] Villani | Topics in Optimal Transportation[END_REF] for more information on p-Wasserstein distances. CLT's in Wasserstein distance have been studied, particularly in dimension 1, where we refer to the works [START_REF] Bobkov | Entropic approach to E. Rio's central limit theorem for w2 transport distance[END_REF][START_REF] Rio | Upper bounds for minimal distances in the central limit theorem[END_REF][START_REF] Rio | Asymptotic constants for minimal distance in the central limit theorem[END_REF] as well as [START_REF] Bobkov | Berry-Esseen bounds and Edgeworth expansions in the central limit theorem for transport distances[END_REF] (and references therein) for convergence rates in p-Wasserstein for all p ≥ 1 under the condition of existence of moments of order 2 + p; in all cases the rate obtained is of optimal order O(1/ √ n). In higher dimensions, results are also available in 2-Wasserstein distance, under more stringent assumptions on the X i . For instance, Zhai [START_REF] Zhai | A multivariate CLT in Wasserstein distance with near optimal convergence rate[END_REF] shows that when X i is almost surely bounded, then a near-optimal rate of convergence in O(log n/ √ n) holds (this improves a result by Valiant et al. [START_REF] Valiant | A CLT and tight lower bounds for estimating entropy[END_REF]). More recently, Courtade et al. [START_REF] Courtade | Existence of Stein kernels under a spectral gap, and discrepancy bound[END_REF] attained the optimal rate of convergence O(n -1/2 ), again in Wasserstein distance of order 2, under the assumption that X i satisfies a Poincarétype inequality; see also [START_REF] Fathi | Stein kernels and moment maps[END_REF] for a similar result under assumption of log-concavity. Finally we mention the work of Bonis [START_REF] Bonis | Rates in the central limit theorem and diffusion approximation via Stein's method[END_REF] where similar estimates are investigated (in Wasserstein-2 distance) under moment assumptions only; dependence of these estimates on the dimension is unclear (see [10, page 12]). One of the key ingredients in many of the more recent above-mentioned references is the multivariate Stein's method. Rates of convergence in the multivariate CLT were first obtained Stein's method by Barbour in [START_REF] Barbour | Stein's method for diffusion approximations[END_REF] (see also Götze [START_REF] Götze | On the rate of convergence in the multivariate CLT[END_REF]) whose methods (which rest on viewing the normal distribution as the stationary distribution of an Ornstein-Uhlenbeck diffusion, and using the generator of this diffusion as a characterizing operator) led to the so-called generator approach to Stein's method with starting point equation ( 6) and its solution given by the classical formula [START_REF] Chatterjee | Multivariate normal approximation using exchangeable pairs[END_REF]. Such an approach readily provides rates of convergence in smooth-k-Wasserstein distances, i.e. integral probability metrics of the form [START_REF] Bentkus | On the dependence of the Berry-Esseen bound on dimension[END_REF] with H (= H (k) ) a set of smooth functions with derivatives up to some order k bounded by 1. Of course, the smaller the order k, the stronger the distance; in particular the case k = 1 coincides with the classical 1-Wasserstein distance (and therefore also [START_REF] Fathi | Stein kernels and moment maps[END_REF] with p = 1). In [7, Theorem 3.1] it is proved that if X i has a finite moment of order 4 then, for any smooth h,

(12) E[h(W ) -h(Z)] 1 √ n 1 2 E|X i | 4 -d ∇h ∞ + √ 2π 3 E|X i | 3 sup x∈R d ∇ 2 h(x) op ,
where M op denotes the operator norm of a matrix M and |x| the Euclidean norm of a vector x ∈ R d . The rate O(1/ √ n) is optimal. The fourth moment conditions are not optimal, nor is the restriction to twice differentiable test functions which implies that [START_REF] Gaunt | Rates of convergence in normal approximation under moment conditions via new bounds on solutions of the Stein equation[END_REF] does not lead to rates of convergence in the 1-Wasserstein distance. Similarly, the bounds on 2-Wasserstein distance recently obtained in [START_REF] Courtade | Existence of Stein kernels under a spectral gap, and discrepancy bound[END_REF][START_REF] Fathi | Stein kernels and moment maps[END_REF] are inspired by concepts related to Stein's method which were introduced in [START_REF] Ledoux | Stein's method, logarithmic Sobolev and transport inequalities[END_REF]; such an approach necessarily requires regularity assumptions on the density of the X i . Hence no simple extension of their approach can lead to rates of convergence in Wasserstein distance with only moment conditions on the X i (and in particular no smoothness assumptions on the densities). In other words, no optimal rates of convergence in Wasserstein distance are available under moment assumptions, and they seem out of reach if based on current available regularity results of Stein's equation.

Contribution.

In this paper, we study Barbour's solution [START_REF] Chatterjee | Multivariate normal approximation using exchangeable pairs[END_REF] to the Stein equation ( 6) and prove new regularity results: namely, if h is α-Hölder for some 0 < α 1, then for all i, j, ∂ 2 f h ∂x i ∂x j is β-Hölder for 0 < β < α. Actually, we show the stronger estimate

(13) ∂ 2 f h ∂x i ∂x j (x) - ∂ 2 f h ∂x i ∂x j (y) = O (|x -y| α log |x -y|) ,
for |x-y| small. A precise statement, with explicit constants (which depend on α and on the dimension d), is given in Proposition 2.2. Note that from Shauder's theory, in the multivariate case (and contrary to the univariate one), one cannot hope in general for the second derivative of f h to inherit the Lipschitz-regularity of h. Actually, Raič [START_REF] Raič | A multivariate CLT for decomposable random vectors with finite second moments[END_REF] gives a counter-example: if

h(x, y) = max{min{x, y}, 0},
then f h defined by ( 7) is twice differentiable but ∂ 2 f h ∂x∂y is not Lipschitz (whereas h is). We study this example in more detail in Proposition 2.4, which shows that (at least for α = 1), the regularity ( 13) cannot be improved in general.

In a second step, we apply those regularity results to estimate the rate of convergence in the CLT, in Wasserstein distance.

Theorem 1.1. Let (X i ) i 1 be an i.i.d. sequence of random vectors with unit covariance matrix, and Z ∼ N (0, I d ). Assume that there exists δ ∈ (0, 1)

such that E[|X i | 2+δ ] < ∞. Then W n -1/2 n i=1 X i , Z 1 n δ 2 (K 1 + 2 (1 -δ) -1 )E|X i | 2+δ + (K 2 + 2d (1 -δ) -1 )E|X i | δ ,
where W stands for the 1-Wasserstein distance, and

K 1 = 2 3/2 2d + 1 d Γ( 1+d 2 ) Γ(d/2) K 2 = 2 2 π √ d.
Note that the rate in O(n -δ/2 ) is optimal when only assuming moments of order 2 + δ; see [START_REF] Bikjalis | Estimates of the remainder term in the central limit theorem[END_REF] or [START_REF] Nefedova | On nonuniform convergence rate estimates in the central limit theorem[END_REF]. As mentioned previously, the bounds in Theorem 1.1 are to our knowledge the first optimal rates in 1-Wasserstein distance in the multidimensional case when assuming finite moments of order 2 + δ only.

From the previous Theorem is easily derived the following Corollary, which gives a near-optimal rate of order O(log n/ √ n) when X i has finite moment of order 3.

Corollary 1.2. Let (X i ) i 1 be an i.i.d. sequence of d-dimensional random vectors with unit covariance matrix. Assume that

E[|X i | 3 ] < ∞. Then for n 3 , W n -1/2 n i=1 X i , Z e C(d) + 2(1 + d) log n √ n E|X i | 3 , where C(d) = 2 3/2 2d+1 d Γ( 1+d 2 ) Γ(d/2) + 2 2 π √ d.
Compared to [START_REF] Zhai | A multivariate CLT in Wasserstein distance with near optimal convergence rate[END_REF], our assumption on the distribution of X i is much weaker; however the distance used in [START_REF] Zhai | A multivariate CLT in Wasserstein distance with near optimal convergence rate[END_REF] is stronger and the constants are sharper ( [START_REF] Zhai | A multivariate CLT in Wasserstein distance with near optimal convergence rate[END_REF] obtains a constant in O( √ d)). [START_REF] Courtade | Existence of Stein kernels under a spectral gap, and discrepancy bound[END_REF] has the advantage of stronger rate of convergence (it is optimal when ours is near-optimal) and stronger distance, but the drawback of a less tractable assumption on the distribution of X i (it should satisfy a Poincaré or weighted Poincaré inequality).

Regularity of solutions of Stein's equation

Throughout the rest of the paper, for x, y ∈ R d , we denote by x • y the Euclidean scalar product between x and y, and |x| the Euclidean norm of x. For a matrix M of size d × d, its operator norm is defined as

M op = sup x∈R d ; |x|=1 |M x|.
Define the α-Hölder semi-norm, for α ∈ (0, 1], by

[h] α = sup x =y |h(x) -h(x)| |x -y| α .
For a multi-index i = (i 1 , . . . , i d ) ∈ N d , the multivariate Hermite polynomial H i is defined by

H i (x) = (-1) |i| e |x| 2 /2 ∂ |i| ∂x i 1 1 . . . ∂x i d d e -|x| 2 /2 ,
where |i| = i 1 + . . . + i d .

Let h : R d → R, and f h be defined by (when the integral makes sense)

(14) f h (x) = - 1 0 1 2t E h(Z x,t ) dt, where h(x) = h(x) -E h(Z), and
Z x,t = √ t x + √ 1 -t Z.
Recall that, when h is smooth with compact support, then [START_REF] Ledoux | Stein's method, logarithmic Sobolev and transport inequalities[END_REF] defines a solution to the Stein equation ( 6), see [START_REF] Barbour | Stein's method for diffusion approximations[END_REF][START_REF] Chatterjee | Multivariate normal approximation using exchangeable pairs[END_REF]. We shall prove that this is still the case when only assuming Hölder-regularity of h. Proposition 2.1. Let h : R d → R be a α-Hölder function; that is, [h] α < ∞. Let f h be the function given by [START_REF] Ledoux | Stein's method, logarithmic Sobolev and transport inequalities[END_REF]. Then:

• f h is twice differentiable and for i = (i 1 , . . . , i d ) ∈ N d such that 1 |i| 2, (15)
∂ |i| f h ∂x i 1 1 . . . ∂x i d d = - 1 0 t |i| 2 -1 2(1 -t) |i| 2 E[H i (Z) h(Z x,t )]dt.
• f h is a solution to the Stein equation (6).

Proof. Fix t ∈ (0, 1). Recall ω(x) = (2π) -d/2 e -|x| 2 /2 is the density of the standard d-dimensional gaussian measure. Since

E h(Z x,t ) = R d h( √ t x + √ 1 -t z) ω(z) dz = (-1) d t d/2 (1 -t) -d/2 R d h(u) ω u - √ t x √ 1 -t du,
we have, from Lebesgue's derivation theorem, and another change of variable, ( 16)

∂ |i| ∂x i 1 1 . . . ∂x i d d E h(Z x,t ) = t |i| 2 (1 -t) |i| 2 E[H i (Z) h(Z x,t )].
Now note that by α-Hölder regularity, and using the fact that EH i (Z 1 , . . . , Z d ) = 0,

|EH i (Z) h(Z x,t )| = |EH i (Z)( h(Z x,t ) -h( √ t x))| E [|H i (Z)| |Z| α ] (1 -t) α/2 . ( 17
)
Thus we can apply again Lebesgue's derivation theorem to obtain [START_REF] Müller | Integral probability metrics and their generating classes of functions[END_REF].

Now let ω t (x) = t -d/2 ω x √ t ; ω 1-t is the density of √ 1 -t Z.
It is well known (and can be easily checked) that ω t solves the heat equation

∂ t ω t = 1 2 ∆ ω t .
We deduce (again applying Lebesgue's derivation theorem, valid since h has polynomial growth at infinity) that

∂ t E h(Z x,t ) = ∂ t R d h(u) ω 1-t (u - √ t x) du = - R d h(u) ∂ t ω 1-t (u - √ t x) du - 1 2 √ t R d h(u)∇ω 1-t (u - √ t x) • x du = - 1 2 R d h(u)∆ ω 1-t (u - √ t x) du - 1 2 √ t R d h(u)∇ω 1-t (u - √ t x) • x du = - 1 2t ∆ x R d h(u) ω 1-t (u - √ t x) du + 1 2t ∇ x R d h(u) ω 1-t (u - √ t x) du • x = - 1 2t (∆ -x • ∇)E h(Z x,t ).
Finally,

h(x) = 1 0 ∂ t E h(Z x,t )dt = - 1 0 1 2t (∆ -x • ∇)E h(Z x,t ) dt = (∆ -x • ∇)f h ,
the last equality being justified by ( 16) and the bound [START_REF] Nourdin | Normal approximations with Malliavin calculus: from Stein's method to universality[END_REF].

Before stating our main regularity results, let us give the idea behind the proof. Starting from (15), we have that

∂ 2 f h ∂x i ∂x j (x) - ∂ 2 f h ∂x i ∂x j (y) = - 1 0 1 2(1 -t) E[(Z i Z j -δ ij )( h(Z x,t ) -h(Z y,t) )]dt.
Using the α-Hölder regularity of h, the modulus of the integrand in last integral can be bounded by

1 2(1 -t) E[|Z i Z j -δ ij | | h(Z x,t ) -h(Z y,t )|] C ij t α/2 1 -t |x -y| α ,
C ij being some constant. However, the function in the right hand-side is not integrable for t close to 1. Thus, for η > 0, we split the integral between 0 and 1η on the one hand (where we can use lour previous bound), and between 1η and 1 on the other hand. To bound the second integral, we remark that since E[Z i Z jδ ij ] = 0, we have that

E[(Z i Z j -δ ij ) h(Z x,t )] = E[(Z i Z j -δ ij )( h(Z x,t ) -h( √ tx))],
which, in modulus, is less than (again using the regularity of h)

E[|Z i Z j -δ ij | Z α ](1 -t) α/2 .
The power (1t) α/2 that is gained makes the integral converge. Finally, we optimize in η > 0.

We are concerned, however, in obtaining the best constants possible (seen as functions of the dimension d and α); this tends to make the proofs more technical that needed if one is only concerned with showing regularity. For this reason, the detailed exposition of the proof in full detail is deferred to Section 5.

We start with the regularity in terms of the operator norm of the Hessian of f h .

Proposition 2.2. Let h : R d → R be a α-Hölder function for some α ∈ (0, 1]. Then the solution f h (7) of the Stein equation (6) satisfies:

(18) ∇ 2 f |x -∇ 2 f |y op [h] α |x -y| α (C 1 (α, d) -2 log |x -y|) , if |x -y| 1 C 1 (α, d) [h] α if |x -y| > 1,
where

(19) C 1 (α, d) = 2 α 2 +1 α + 2d αd Γ( α+d 2 ) Γ(d/2)
.

In particular, for all 0 < β < α , ∂ 2 f h ∂x i ∂x j is globally β-Hölder:

(20) ∇ 2 f |x -∇ 2 f |y op C 1 (α, d) + 2 α -β |x -y| β [h] α .
It also holds the (1 + log) α-Hölder regularity

(21) ∇ 2 f |x -∇ 2 f |y op |x -y| α (C 1 (α, d) + | log |x -y| |) [h] α .
Now we turn to the regularity of the Laplacian.

Proposition 2.3. Let h : R d → R be a α-Hölder function for some α ∈ (0, 1]. Then the solution f h (7) of the Stein equation (6) satisfies:

∆f |x -∆f |y [h] α |x -y| α (C 2 (α, d) -2 d log |x -y|) , if |x -y| 1 C 2 (α, d) [h] α if |x -y| > 1,
where

(22) C 2 (α, d) = 2 α 2 +1 (α+2d) Γ( α+d 2 ) α Γ(d/2) if α ∈ (0, 1), C 2 (1, d) = 2 2 π √ d.
In particular, for all 0 < β < α,

(23) |∆f |x -∆f |y | C 2 (α, d) + d 2 α -β |x -y| β [h] α .
Note that Proposition 2.2 implies that, for α = 1, when |x -y| is small, then

∂ 2 f h ∂x i ∂x j (x) - ∂ 2 f h ∂x i ∂x j (y) = O (|x -y| log |x -y|) .
The example given by Raič [START_REF] Raič | A multivariate CLT for decomposable random vectors with finite second moments[END_REF] shows that this rate is optimal; indeed, we have the following result.

Proposition 2.4. Let h : R 2 → R be the Lipschitz function defined by h(x, y) = max(0, min(x, y)).

Then ∂ 2 f h ∂x∂y (u, u) - ∂ 2 f h ∂x∂y (0, 0) ∼ u→0 + 1 √ 2π u log u.
The proof can be found in the Appendix.

Multivariate Berry-Esseen bounds in Wasserstein distance

As anticipated, we apply the regularity results obtained in previous section to obtain Berry-Esseen bounds in the CLT, in 1-Wasserstein distance.

Let X 1 , X 2 , ... be an i.i.d. sequence of centered, square-integrable and isotropic random vectors; that is,

E[X 1 ] = 0 and E[X 1 X T 1 ] = I d . Let W = n -1/2 n i=1 X i .
We are interested in W α (W, Z) for α ∈ (0, 1], where the α-Wasserstein distance is defined as

W α (X, Y ) = sup {h∈C(R n ,R)|[h]α 1} Eh(X) -Eh(Y ).
As in the introduction, for α = 1, the resulting distance is

W := W 1 , the classical 1-Wasserstein distance (that is, W(X, Y ) = sup H Eh(X) -Eh(Y ) with H the collection of 1-Lipschitz functions).
We are now in a position to prove our main Theorem. We first give a more general version of it in α-Wasserstein distances; Theorem 1.1 is just the following Theorem applied to α = 1. Theorem 3.1. Let α ∈ (0, 1] and (X i ) i 1 be an i.i.d. sequence of d-dimensional random vectors with unit covariance matrix. Assume that there exists δ ∈ (0, α)

such that E[|X i | 2+δ ] < ∞. Then W α n -1/2 n i=1 X i , Z 1
n δ 2 C 1 (α, d) + 2 α -δ E|X i | 2+δ + C 2 (α, d) + d 2 α -δ E|X i | δ ,
where C 1 (α, d) and C 2 (α, d) are respectively defined in [START_REF] Reinert | Multivariate normal approximation with Stein's method of exchangeable pairs under a general linearity condition[END_REF] and [START_REF] Stein | A bound for the error in the normal approximation to the distribution of a sum of dependent random variables[END_REF].

Remark 3.2. From Stirling's formula, C 1 (α, d) = O( √ d), C 2 (α, d) = O(d 1+α/2 ) for α ∈ (0, 1) and C 2 (1, d) = O( √ d).
Proof of Theorems 1.1 and 3.1 . Let h be α-Hölder (with [h] α 1) and f h be the solution of the Stein equation defined by Proposition 2.1. Then,

E[h(W ) -h(Z)] = E[∆f h (W ) -W • ∇f h (W )] = 1 n n i=1 E[∆f h (W ) - √ nX i • ∇f h (W ) .
The following calculations already appeared in the literature (see e.g. [START_REF] Raič | A multivariate CLT for decomposable random vectors with finite second moments[END_REF]), we include them here for completeness.

Let W i = W -X i / √ n = 1 √ n j =i X j
. By Taylor's formula, we have for some uniformly distributed in [0, 1] (and independent of everything else) θ

E[X i • ∇f h (W )] = 1 √ n E X T i ∇ 2 f h W i + θ X i √ n X i , leading to E[h(W ) -h(Z)] = 1 n n i=1 E ∆f h (W ) -X T i ∇ 2 f h W i + θ X i √ n X i .
Let X i,j be the jth coordinate of X i . Since W i is independent of X i , and X i has unit covariance matrix, we have

E X T i ∇ 2 f h (W i ) X i = d j,k=1 E X i,j X i,k ∂ 2 f h ∂x j ∂x k (W i ) = d j=1 E ∂ 2 f h ∂x 2 j (W i ) = E[∆f h (W i )].
Finally,

E[h(W ) -h(Z)] = 1 n n i=1 E ∆f h (W ) -∆f h (W i ) -X T i ∇ 2 f h W i + θ X i √ n -∇ 2 f h (W i ) X i .
Note then that

X T i ∇ 2 f h W i + θ X i √ n -∇ 2 f h (W i ) X i |X i | 2 ∇ 2 f h W i + θ X i √ n -∇ 2 f h (W i ) op .
We now use Propositions 2.2 and 2.3 with β = δ α to obtain:

E[h(W ) -h(Z)] 1 n n i=1 C 2 (α, d) + d 2 α -δ E[|X i | δ θ δ ] n δ/2 + C 1 (α, d) + 2 α -δ E[|X i | 2+δ θ δ ] n δ/2 .
Noting that E[θ β ] 1 and rearranging, we obtain the result.

Corollary 3.3. Let α ∈ (0, 1], and (X i ) i 1 be an i.i.d. sequence of d-dimensional random vectors with unit covariance matrix. Assume that

E[|X i | 2+α ] < ∞. Then for n > exp(2/α) , W α n -1/2 n i=1 X i , Z e C 1 (α, d) + C 2 (α, d) + 2(1 + d) log n n α 2 E|X i | 2+α ,
where C 1 (α, d) and C 2 (α, d) are respectively defined in [START_REF] Reinert | Multivariate normal approximation with Stein's method of exchangeable pairs under a general linearity condition[END_REF] and [START_REF] Stein | A bound for the error in the normal approximation to the distribution of a sum of dependent random variables[END_REF].

Proof. By Hölder's and the Cauchy-Schwarz inequalities, for any δ α,

E|X i | δ E|X i | 2+α δ/(2+α) . But by Jensen's inequality, E|X i | 2+α (E|X i | 2 ) (2+α)/2 = d (2+α)/2
1, so that, since δ/(2 + α) < 1,

E|X i | 2+α δ/(2+α) E|X i | 2+α . Similarly, E|X i | 2+δ E|X i | 2+α 1-α-δ 2+α E|X i | 2+α
. Note now that the bound of Theorem 3.1 holds for any 0 < δ < α. Choosing αδ = 2/ log n achieves the proof since n -1 log n = 1/e.

When applied to α = 1, previous corollary leads to Corollary 1.2, which we recall here: as long as

E|X i | 3 < ∞, W n -1/2 n i=1 X i , Z e C(d) + 2(1 + d) log n √ n E|X i | 3 ,
where

C(d) = 2 3/2 (2d+1) Γ( d+1 2 ) dΓ( d 2 ) + 2 2 π √ d.
[25] also obtains a near-optimal rate of convergence in O(log n/ √ n), but under the much stronger assumption that |X i | β almost surely; nevertheless, the distance used in [START_REF] Zhai | A multivariate CLT in Wasserstein distance with near optimal convergence rate[END_REF] (the quadratic Wasserstein distance) is stronger than ours, the behaviour of the constant on the higher order term is O( √ d), here we obtain O(d).

Extension to higher order derivatives

The regularity result easily extends to higher order derivatives.

∂ p+2 f p+2 j=1 ∂x i j (x) - ∂ p+2 f p+2 j=1 ∂x i j (y) [h] α,p |x -y| α (A -2 log |x -y|) , if |x -y| 1 A [h] α,p if |x -y| > 1,
where

A = 2 α/2+1 α + d + 1 α Γ( α+d 2 ) Γ(d/2)
.

In particular, all derivatives of the order p + 2 of f h are β-Hölder for any 0 < β < α and we have

∂ p+2 f p+2 j=1 ∂x i j (x) - ∂ p+2 f p+2 j=1 ∂x i j (y) [h] α,p A + 2 α -β |x -y| β .
Proof. Taking derivatives in [START_REF] Chatterjee | Multivariate normal approximation using exchangeable pairs[END_REF], we have

∂ p+2 f p+2 j=1 ∂x i j (x) = 1 0 t p 2 E ∂ p+2 h p+2 j=1 ∂x i j (Z x,t ) dt.
Next perform two Gaussian integration by parts against two indices i p+1 and i p+2 , say, to get

∂ p+2 f p+2 j=1 ∂x i j (x) = 1 0 t p 2(1 -t) E (Z i p+1 Z i p+2 -δ i p+1 i p+2 ) ∂ p h p j=1 ∂x i j (Z x,t ) dt.
Then, using the same method as in the proof of Proposition 2.2 (we do not give all the details here), we have

∂ p+2 f p+2 j=1 ∂x i j (x) - ∂ p+2 f p+2 j=1 ∂x i j (y) -[h] α,p log η + [h] α,p 2 α/2+1 α + d + 1 α Γ( α+d 2 ) Γ(d/2) η α/2 .
Choose η = |x-y| if |x-y| 1, 1 otherwise, to get the first result, and the fact thatlog u 1 α-β u β-α if u 1 for the second one.

We stress that one possible application of this Proposition would be a multivariate Berry-Esseen bound in the CLT with matching moments (i.e. assuming that the underlying random variables X i share the same first k moments with the Gaussian). In this case, faster rates of convergence are expected, see [START_REF] Gaunt | Rates of convergence in normal approximation under moment conditions via new bounds on solutions of the Stein equation[END_REF].

The remaining proofs

Proof of Proposition 2.2. Recall that

∂ 2 f h ∂x i ∂x j = - 1 0 1 2(1 -t) E[(Z i Z j -δ ij ) h(Z x,t )]dt. Since E[Z i Z j -δ ij ] = 0, we have E[(Z i Z j -δ ij ) h( √ tx)] = 0, so that E[(Z i Z j -δ ij ) h(Z x,t )] = E[(Z i Z j -δ ij )( h(Z x,t ) -h( √ tx))]. (24)
Thus,

∇ 2 f h (x) = - 1 0 1 2(1 -t) E[(ZZ T -I d )( h(Z x,t ) -h( √ tx))] dt,
where Z T denotes the transpose of Z. Let a = (a 1 , . . . , a d ) T ∈ R d and assume that |a| = 1. We have

a T ∇ 2 f h (x)a = - 1 0 1 2(1 -t) E[a T (ZZ T -I d ) a ( h(Z x,t ) -h( √ tx))]dt = - 1 0 1 2(1 -t) E[((Z • a) 2 -1) ( h(Z x,t ) -h( √ tx))]dt. Since | h(Z x,t ) -h( √ tx)| [h] α (1 -t) α/2 Z α , we also have |E[((Z • a) 2 -1) ( h(Z x,t ) -h( √ tx))] | [h] α E[|(a • Z) 2 -1| Z α ](1 -t) α/2 . (25) Let us bound E[|(a • Z) 2 -1| Z α ].
Let (a, e 2 , . . . , e d ) be an orthonormal basis and Z = (a

• Z, e 2 • Z, . . . , e d • Z) T . Then Z ∼ N (0, I d ). Moreover, E[|(a • Z) 2 -1| Z α ] = E[| Z2 1 -1| Z α ]. Thus, E[|(a • Z) 2 -1| Z α ] = E[| Z2 1 -1| Z α ] E[( Z2 1 + 1) Z α ] = 1 d d i=1 E[( Z2 i + 1) Z α ] = 1 d E[( Z 2 + d) Z α ]. For all β > 0, E Z β = 2 β 2 Γ( β+d 2 )
Γ(d/2) . We define

(26) C = 1 d E[( Z 2 + d) Z α ] = 2 α 2 α + 2d d Γ( α+d 2 ) Γ(d/2)
.

This shows in particular that ∇ 2 f h (x) op is bounded. Now we consider a T ∇ 2 f h (x) -∇ 2 f h (y) a and split the integral into two parts. Let η ∈ [0, 1]. We have

|a T ∇ 2 f h (x) -∇ 2 f h (y) a| = 1 0 1 2(1 -t) E a T (ZZ T -I d ) a( h(Z x,t ) -h(Z y,t ) dt 1-η 0 1 2(1 -t) E |a T (ZZ T -I d ) a| | h(Z x,t ) -h(Z y,t )| dt + 1 1-η 1 2(1 -t)
E a T (ZZ T -I d ) a( h(Z x,t ) -h(Z y,t )) dt.

Using the α-Hölder regularity of h for the first part of the integral and (24) twice in the second part together with [START_REF] Zhai | A multivariate CLT in Wasserstein distance with near optimal convergence rate[END_REF] and (26), we can bound the previous quantity by

[h] α |x -y| α E |(a • Z) 2 -1| 1-η 0 t α/2 2(1 -t) dt + [h] α C 1 1-η (1 -t) -1+α/2 dt (27) [h] α -|x -y| α log η + 2C α η α/2 , ( 28
)
where to obtain (28), we used the facts that E |(a • Z) 2 -1| 2 and t α/2 1. Choose η = |x -y| 2 if |x -y| 1, η = 1 otherwise to get [START_REF] Raič | A multivariate CLT for decomposable random vectors with finite second moments[END_REF]. Equation ( 21) is a straightforward reformulation since 1 + | log(u)| 1. To get [START_REF] Rio | Upper bounds for minimal distances in the central limit theorem[END_REF], simply note that for 0 < β < α and 0 < u 1,log u 1 α-β u β-α and for 1 u, 1 u β .

Proof of Proposition 2.3. The regularity of the Laplacian is proved in a similar manner as for the operator norm of the Hessian; we do not detail the computations here. Let α ∈ (0, 1). We have whose modulus can be thus bounded by

|∆f
2 √ 1 -t E[ Z ] = √ 1 -t 2 √ 2 √ π √ d.
Using this bound in the integral between 1η and 1, and choosing η as in Proposition 2.2, we obtain the results.

Proof of Proposition 2.4. Let u > 0. Denote Z t,u i = √ t u + √ 1t Z i . We have

∂ 2 f h ∂x∂y (u, u) = - 1 0 1 2(1 -t) E[Z 1 Z 2 h(Z t,u 1 , Z t,u 2 )] dt = - 1 0 1 2(1 -t) E Z 1 Z 2 (1 Z t,u 2 Z t,u 1 0 Z t,u 1 + 1 Z t,u 1 Z t,u 2 0 Z t,u 2 ) dt = - 1 0 1 1 -t E Z 1 Z 2 1 Z t,u 2 Z t,u 1 0 Z t,u 1 dt = - 1 √ 2π 1 0 1 1 -t E Z 1 e -Z 2 1 2 1 Z t,u 1 0 ( √ t u + √ 1 -t Z 1 ) dt,

Proposition 4 . 1 .

 41 Let h : R d → R be a smooth, compactly supported function, and denote by [h] α,p a common α-Hölder constant for all derivatives of order p of h. Then the solution f h (7) of the Stein equation (6) satisfies, for all (i 1 , . . . , i p+2 ) ∈ {1, . . . , d} p+2 :

  [START_REF] Fathi | Stein kernels and moment maps[END_REF]) is of the form[START_REF] Chen | Normal approximation by Stein's method[END_REF] with H the class of Lipschitz function with constant 1. Because W p (•, •) ≥ W p ′ (•, •) for all p ≥ p ′ , bounds in p-Wasserstein distance are stronger than those in p ′ -Wasserstein distance; in particular

  |x -∆f |y | -1)( h(Z x,t ) -h(Z y,t ) dt -1)( h(Z x,t ) -h(Z y,t )) dt [h] α |x -y| α E Z 2 + d -d |x -y| α log η + 2E[( Z 2 + d) Z α ] α η α/2 .and choose again η = |x -y| 2 if |x -y| 1, η = 1 otherwise.We can obtain better constants in the case α = 1. Indeed, note that by using only one integration by parts,

	= 0	1 i 1-η 0 d 1 2(1 -t) E i=1 (Z 2 d 1 2(1 -t) E i=1 (Z 2 i -1) | h(Z x,t ) -h(Z y,t )| dt
	+	1 1-η	1 2(1 -t)	E	d i=1	(Z 2 i 1-η 0	t α/2 2(1 -t)	dt + [h] α E[( Z 2 + d) Z α ]	1 1-η	(1 -t) -1+α/2 dt
	[h] Note that
		E[( Z 2 + d) Z α ] =	2	α 2 +1 Γ( α+d 2 + 1) + d 2 Γ(d/2)	α 2 Γ( α+d 2 )	= 2	α 2 (α + 2d)	Γ( α+d 2 ) Γ(d/2)	,
			d
		E	(Z 2
			i=1

α i -1)( h(Z x,t ) -h(Z y,t )) = √ 1 -tE d i=1 Z i (∂ i h(Z x,t ) -∂ i h(Z y,t )) = √ 1 -tE [Z • (∇h(Z t,x ) -∇h(Z t,y )] ,

Acknowledgements

We would like to thank Guillaume Carlier for suggesting the problem to us and subsequent useful discussions. We also thank Max Fathi for useful discussions. The three authors were supported by the Fonds de la Recherche Scientifique -FNRS under Grant MIS F.4539.16.

It is readily checked that the last integral is equivalent to -2 log u, when u → 0 + . On the other hand, by Fubini's theorem, we have

which is a O(u) as u → 0 + . This achieves the proof.


Alexandre Francisco

João Carriço

Maxime Crochemore

Solon Pissis

Bruno Ribeiro-Gonçalves

Ph.D Cátia Vaz

Keywords: computational biology, phylogenetic inference, Hamming distance

Distance-based phylogeny algorithms usually assume a matrix of pairwise distances between different taxa as input. However, there are algorithms that need only distances betweeen taxa that are sufficiently small, say smaller than some given k. The question is how to construct such a restricted similarity matrix as fast as possible. The submitted article tackles this question with Hamming distance as the measure between the genotype profile sequence representing each taxa. An O(md) average time optimal algorithm is given for small k, where d is the number of taxa and m is their length. The algorithm is a simple application of suffix arrays enhanced with LCP information and RMQ data structure.

The main observation is that one can afford to output all length L>=m/k matching substring pairs between all profiles, and check which ones lead to real matches in O(k) time using constant time longest common extension queries implemented by RMQ on LCP array. The probability of false positives is small enough for small enough k, so that the running time is dominated by that of producing the output. Similar analyses have been conducted earlier for approximate string matching (Fredriksson and Navarro. Average-Optimal Multiple Approximate String Matching. CPM 2003).

The resulting algorithm is plugged into an existing phylogeny tool that can exploit restricted matrices. Experiments show that the speed-up is significant and the simulations also confirm that the average case analysis assumptions are not too optimistic.

A new application for the technique has been included to this extended journal version compared to the original conference paper. Namely, the application of querying typing databases is considered. In this application, a query pattern is searched in a database for approximate matches. The specialization of the all-pairs algorithm for this application is implemented and incorporated to INNUENDO Platform.

Minor revision requests:

Please add a reference to the approximate pattern matching literature, e.g. [Fredriksson and Navarro. Average-Optimal Multiple Approximate String Matching. CPM 2003], is one candidate, but some earlier work already contains similar analyses. With the added application of searching for query patterns this connection is even more evident. The analyses are really pretty much the same including the limitations on k when the approach works, so this connection should really be made visible, to give credit to the earlier work. R: We agree with the reviewer and we added the suggested reference, and another one, stating the relationship between the problems addressed in our work and approximate pattern matching. Both the statement and cited references may be found in Conclusions. Page 12, first line: Add "which is" or something like this to make the sentence complete.

R: Thanks, we fixed the sentence adding "which is" as suggested.

References: ??? in many places R: Thanks, some info was missing in our bib file. We revised all references and we believe that all are complete now.

With my best regards,

Alexandre Francisco

Powered by Editorial Manager® and ProduXion Manager® from Aries Systems Corporation

RESEARCH

Fast phylogenetic inference from typing data

João A Carriço 1 , Maxime Crochemore 2 , Alexandre P Francisco 3,4* , Solon P Pissis 2 , Bruno Ribeiro-Gonçalves 1 and Cátia Vaz 3,5 Background

Introduction

The evolutionary relationships between different species or taxa are usually inferred through known phylogenetic analysis techniques. Some of these techniques rely on

Manuscript

Click here to download Manuscript paper.tex

Click here to view linked References the inference of phylogenetic trees, which can be computed from DNA or Protein sequences, or from allelic profiles where the sequences of defined loci are abstracted to categorical indexes. The most popular method is MultiLocus Sequence Typing (MLST) [START_REF] Maiden | Multilocus sequence typing: a portable approach to the identification of clones within populations of pathogenic microorganisms[END_REF] that typically uses seven 450 to 700 bp fragments of housekeeping genes for a given species. Phylogenetic trees are also used in other contexts, such as to understand the evolutionary history of gene families, to allow phylogenetic footprinting, to trace the origin and transmission of infectious diseases, or to study the co-evolution of hosts and parasites [START_REF] Huson | Phylogenetic Networks: Concepts, Algorithms and Applications[END_REF][START_REF] Robinson | Bacterial Population Genetics in Infectious Disease[END_REF].

In traditional phylogenetic methods, the process of phylogenetic inference starts with a multiple alignment of the sequences under study that is then corrected using models of DNA or Protein evolution. Tree-building methodologies can then be applied on the resulting distance matrix. These methods rely on some distancebased analysis of sequences or profiles [START_REF] Saitou | Introduction to Evolutionary Genomics[END_REF].

Distance-based methods for phylogenetic analysis rely on a measure of genetic evolution distance, which is often defined directly or indirectly from the fraction of mismatches at aligned positions, with gaps either ignored or counted as mismatches. A first step of these methods is to compute this distance between all pairs of sequences. The simplest approach is to use the Hamming distance, also known as observed p-distance, defined as the number of positions at which two aligned sequences differ. Note that the Hamming distance between two sequences underestimates their true evolutionary distance and, thus, a correction formula based on some model of evolution is often used [START_REF] Huson | Phylogenetic Networks: Concepts, Algorithms and Applications[END_REF][START_REF] Saitou | Introduction to Evolutionary Genomics[END_REF]. Although distance-based methods not always produce the best tree for the data, usually they also incorporate an optimality criterion into the distance model for getting more plausible phylogenetic reconstructions, such as the minimum evolution criterion [START_REF] Desper | Fast and accurate phylogeny reconstruction algorithms based on the minimum-evolution principle[END_REF], the least squares criterion [START_REF] Pardi | Distance-based methods in phylogenetics[END_REF] or the clonal complexes expansion and diversification [START_REF] Feil | Recombination within natural populations of pathogenic bacteria: short-term empirical estimates and long-term phylogenetic consequences[END_REF]. Nevertheless, this category of methods are much faster than Maximum Likelihood or Bayesian

Inference Methods [START_REF] Yang | Molecular phylogenetics: principles and practice[END_REF], making them excellent choices for the primary analysis of large data sets.

Most of the distance-based methods are agglomerative methods. They start with each sequence being a singleton cluster and, at each step, they join two clusters. The iterative process stops when all sequences are part of a single cluster, resulting in a phylogenetic tree. At each step the candidate pair is selected taking into account the distance among clusters as well as the optimality criterion chosen to adjust it.

The computation of a distance matrix (2D array containing the pairwise distances between the elements of a set) is a common first step for distance-based methods, such as eBURST [9], goeBURST [START_REF] Francisco | Global optimal eBURST analysis of multilocus typing data using a graphic matroid approach[END_REF], Neighbor Joining [START_REF] Saitou | The neighbor-joining method: a new method for reconstructing phylogenetic trees[END_REF] and UPGMA [START_REF] Sokal | A statistical method for evaluating systematic relationships[END_REF]. This particular step dominates the running time of most methods, taking Θ(md 2 ) time in general, d being the number of sequences or profiles and m the length of each sequence or profile. For large-scale datasets this running time may be quite problematic. And nowadays, with the mainstream use of High Throughput Sequencing, the amount of data being accumulated in typing databases is huge. It is common to find databases storing thousands of different profiles for a single microbial species, with each profile having thousands of loci [START_REF] Sergean | EnteroBase[END_REF][START_REF] Jolley | BIGSdb: Scalable analysis of bacterial genome variation at the population level[END_REF].

However, depending on application, on the underlying model of evolution and on the optimality criterion, it may not be strictly necessary to be aware of the complete distance matrix. There are methods that continue to provide optimal solutions without a complete matrix. For such methods, one may still consider a truncated distance matrix and several heuristics, combined with final local searches through topology rearrangements, to improve the running time [START_REF] Pardi | Distance-based methods in phylogenetics[END_REF]. The goeBURST algorithm, one of our use cases in this article, is an example of a method that can work with truncated distance matrices by construction, i.e., one needs only to know which pairs are at Hamming distance at most k.

Our results

We propose here an average-case O(md)-time and O(md)-space algorithm to compute the pairs of sequences, among d sequences of length m, that are at distance at

most k, when k < (m-k-1)•log σ log md
, where σ is the size of the sequences alphabet. We support our result with both a theoretical analysis and an experimental evaluation on synthetic and real datasets of different data types (MLST, cgMLST, wgMLST and SNP). We further show that our method improves goeBURST, and that we can use it to speedup querying local phylogenetic patterns over large typing databases.

A preliminary version of this paper was presented at the Workshop on Algorithms in Bioinformatics (WABI) 2017 [START_REF] Crochemore | Towards Distance-Based Phylogenetic Inference in Average-Case Linear-Time[END_REF]. We address this problem by indexing all profiles P using the suffix array (denoted by SA) and the longest common prefix (denoted by LCP) array [START_REF] Manber | Suffix arrays: a new method for on-line string searches[END_REF]. We rely also on a range minimum queries (RMQ) data structure [START_REF] Bender | The LCA problem revisited[END_REF][START_REF] Bender | Lowest common ancestors in trees and directed acyclic graphs[END_REF] over the LCP array (denoted by RMQ LCP ). The problem is then solved in three main steps:

1 Index all profiles using the SA data structure.

2 Enumerate all candidate profile pairs given the maximum Hamming distance k.

3 Verify each candidate profile pair by checking if the associated Hamming distance is no more than k.

Table 1 summarizes the data structures and strategies followed in each step. Profiles are concatenated and indexed using SA. Depending on the strategy to be used, we further process the SA and build the LCP array and pre-process it for fast RMQ.

This allows for enumerating candidate profile pairs and computing distances faster.

In what follows, we detail the above steps and show how the data structures are used to improve the overall running time.

Step 1: Profile indexing Profiles are concatenated and indexed in an SA in O(md) time and space [START_REF] Kärkkäinen | Linear work suffix array construction[END_REF][START_REF] Ko | Space efficient linear time construction of suffix arrays[END_REF].

Let us denote this string by s. Since we only need to compute the distances between profiles that are at Hamming distance at most k, we can conceptually split each profile into k non-overlapping blocks of length L = ⌊ m k+1 ⌋ each. It is then folklore knowledge that if two profiles are within distance k, they must share at least one such block of length L. Our approach is based on using the SA of s to efficiently Another solution relies on computing the LCP array: the longest common prefix between each pair of consecutive elements within the SA. This information can also be computed in O(n) time and space [START_REF] Kasai | Linear-time longest-common-prefix computation in suffix arrays and its applications[END_REF]. Since SA is an ordered structure, for the contiguous suffixes s i , s i+1 , s i+2 of s, with 0 ≤ i < n -2, we have that the common prefix between s i and s i+1 is at least as long as the common prefix of s i and s i+2 .

By construction, it is possible to get the position of each suffix in the corresponding profile in constant time. Then, we cluster the corresponding profiles of contiguous pairs if they have an LCP value greater than or equal to L and they are also aligned.

This clustering procedure can be done in O(kd 2 ) time.

Step

3: Pairs verification

After getting the set of candidate profile pairs, a naïve solution would be to compute the distance for each pair of profiles by comparing them in linear time, i.e., O(m) time. However, if we compute the LCP array of s, we can then perform a sequence of O(k) RMQ over the LCP array for checking if a pair of profiles is at distance at most k. These RMQ over the LCP array correspond to longest common prefix queries between a pair of suffixes of s. Since after a linear-time pre-processing over the LCP array, RMQ can be answered in constant time per query [START_REF] Bender | The LCA problem revisited[END_REF], we obtain a faster approach for computing the distances. This alternative approach takes O(k) time to verify each candidate profile pair instead of O(m) time.

Average-case analysis

Algorithm 1 below details the solution based on LCP clusters; and Theorem 1 shows that this algorithm runs in linear time on average using linear space. We rely here on well-known results concerning the linear-time construction of the SA [START_REF] Kärkkäinen | Linear work suffix array construction[END_REF][START_REF] Ko | Space efficient linear time construction of suffix arrays[END_REF] and the LCP array [START_REF] Kasai | Linear-time longest-common-prefix computation in suffix arrays and its applications[END_REF], as well as the linear-time pre-processing for the RMQ data structure [START_REF] Bender | Lowest common ancestors in trees and directed acyclic graphs[END_REF].

In what follows, LCP[i], i > 0, stores the length of the longest common prefix of suffixes s i-1 and s i of s, and RMQ LCP (i, j) returns the index of the smallest element in the subarray LCP[i . . . j] in constant time [START_REF] Bender | Lowest common ancestors in trees and directed acyclic graphs[END_REF]. We rely also on some auxiliary subroutines; let L = ⌊ m k+1 ⌋: Aligned(i) Let ℓ = i mod m, i.e., the starting position of the suffix s i within a profile. Then this subroutine returns ℓ/L if ℓ is multiple of L, and -1 otherwise.

HD(p i , p j , ℓ) Given two profiles p i and p j which share a substring of length L, starting at index ℓL, this subroutine computes the minimum of k and the Hamming distance between p i and p j . This subroutine relies on RMQ LCP to find matches between p i and p j and, hence, it runs in O(k) time since it can terminate after k mismatches. Let us denote by B the total number of blocks over s and by L the block length.

We set L = ⌊ m k+1 ⌋ and thus we have that B = d⌊ m L ⌋. Let us also denote by C a maximal set of indices over x satisfying the following:

1 the length of the longest common prefix between any two suffixes of s starting at these indices is at least L; Algorithm 1: Algorithm using LCP clusters.

1 Input: A set P of d profiles of length m each; an integer threshold 0 < k < m.

2 Output: The set X of distinct pairs of profiles that are at Hamming distance at most k, i.e., X = {(u, v) ∈ P × P | u < v and H(u, v) ≤ k}.

3 Initialization: Let s = s[0 . . . n -1] be the string of length n = md obtained after concatenating the d profiles, and L = ⌊ m k+1 ⌋. Construct the SA S for s, the LCP array for s and RMQ LCP . Initialize a hash table HT to track verified pairs. 4 Candidate pairs enumeration: δ := HD(p, q, t)

5 X := ∅; ℓp := -1; Ct := ∅, for 0 ≤ t ≤ k 6 foreach 1 ≤ i < n do 7 ℓ := LCP[i] 8 if ℓ ≥ L then 9 p i := ⌊S[i]/m⌋ 10 x := Aligned(i) 11 if x = -1 then 12 Cx := Cx ∪ {p i } 13 if ℓp = -1 then 14 p i-1 := ⌊S[i -1]/m⌋ 15 x := Aligned(i -1) 16 if x = -1 then 17 Cx := Cx ∪ {p i-1 }
26 if δ ≤ k then 27 X := X ∪ {(p, q)} 28 ℓp := -1; Ct := ∅, for 0 ≤ t ≤ k 29 Finalize: Return the set X.
2 both of these suffixes start at the starting position of a block; This can be done in O(md) time using the LCP array (lines 7-17 where PROC i,j is the time required to process a pair i, j of elements of a set C, and Pairs is the sum of |C| 2 over all such sets C. We have that PROC i,j = O(k) by using RMQ over the LCP array. Additionally, by the stated assumption on the d profiles, the expected value for Pairs is no more than Bd σ L : we have B blocks in total and each block can only match at most d other blocks by the conditions above.

Hence, the algorithm requires on average the following running time

O(md + k • Bd σ L ).
Let us analyze this further to obtain the relevant condition on k. We have the following:

k • Bd σ L = k • ⌊ m ⌊m/(k+1)⌋ ⌋ • d 2 σ ⌊ m k+1 ⌋ ≤ k • ( m ⌊m/(k+1)⌋ ) • d 2 σ m k+1 -1 .
Since 0 < k < m by hypothesis, we have the following:

k • ( m ⌊m/(k+1)⌋ ) • d 2 σ m k+1 -1 ≤ (md) 2 σ m k+1 -1 .
By some simple rearrangements we have that:

(md) 2 σ m k+1 -1 = (md) 2 (md) log σ log md ( m k+1 -1) = (md) 2-(m-k-1) log σ (k+1) log md .
Consequently, in the case when

k < (m -k -1) • log σ log md
the algorithm requires O(md) time on average. The extra space usage is clearly O(md).

Use case 1: goeBURST algorithm

The distance matrix computation is a main step in distance-based methods for phylogenetic inference. This step dominates the running time of most methods, taking Θ(md 2 ) time, for d sequences of length m, since it must compute the distance among all sequence pairs. But for some methods, or when we are only interested in local phylogenies for sequences or profiles of interest, one does not need to know all pairwise distances for reconstructing a phylogenetic tree. The problem addressed in this article was motivated by the goeBURST algorithm [START_REF] Francisco | Global optimal eBURST analysis of multilocus typing data using a graphic matroid approach[END_REF], our use case 1.

goeBURST is one of such methods for which one must know only the pairs of sequences that are at Hamming distance at most k. The solution proposed here can however be extended to other distance-based phylogenetic inference methods, that rely directly or indirectly on Hamming distance computations. Note that most methods either consider the Hamming distance or its correction accordingly to some formula based on some model of evolution [START_REF] Huson | Phylogenetic Networks: Concepts, Algorithms and Applications[END_REF][START_REF] Saitou | Introduction to Evolutionary Genomics[END_REF]. In both cases we must start by computing the Hamming distance among sequences, but not necessarily all of them [START_REF] Pardi | Distance-based methods in phylogenetics[END_REF].

The underlying model of goeBURST is as follows: a given genotype increases in frequency in the population as a consequence of a fitness advantage or of random genetic drift, becoming a founder clone in the population; and this increase is accompanied by a gradual diversification of that genotype, by mutation and recombination, forming a cluster of phylogenetic closely-related strains. This diversification of the "founding" genotype is reflected in the appearance of genetic profiles differing only in one housekeeping gene sequence from this genotype -single locus variants (SLVs). Further diversification of those SLVs will result in the appearance of variations of the original genotype with more than one difference in the allelic profile, e.g., double and triple locus variants (DLVs and TLVs).

The problem solved by goeBURST can be stated as a graphic matroid optimization problem and, hence, it follows a classic greedy approach [START_REF] Papadimitriou | Combinatorial Optimization: Algorithms and Complexity[END_REF]. Given the maximum Hamming distance k, we can define a graph G = (V, E), where V = P (set of profiles) and

E = {(u, v) ∈ V 2 | H(u, v) ≤ k}.
The main goal of goeBURST is then to compute a minimum spanning forest for G taking into account the distance H and a total order on links. It starts with a forest of singleton trees (each sequence/profile is a tree). Then it constructs the optimal forest by adding links connecting profiles in different trees in increasing order accordingly to the total order, similarly to what is done in the Kruskal's algorithm [START_REF] Kruskal | On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF]. In the current implementation, a total order for links is implicitly defined based on the distance between sequences, on the number of SLVs, DLVs, TLVs, on the occurrence frequency of sequences, and on the assigned sequence identifier. With this total order, the construction of the tree consists of building a minimum spanning forest in a graph [START_REF] Kruskal | On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF], where each sequence is a node and the link weights are defined by the total order. By construction, the pairs at distance δ will be joined before the pairs at distance δ + 1.

Use case 2: Querying typing databases

A related problem is querying typing databases for similar typing profiles. Given a set P of d profiles of length m each, a profile u not necessarily in P but with the same length m as those in P , and k such that 0 < k < m, the problem is to find all profiles v ∈ P such that H(u, v) ≤ k. One may be also interested on local phylogenetic patterns, but those can be inferred from found profiles using for instance the goeBURST algorithm.

Once we define the value for k, we can address this problem as follows. We index all d profiles in the database as before in linear time O(md), and given a query profile u, we enumerate all candidate profiles v. We then verify as before all candidate pairs and we return only those satisfying H(u, v) ≤ k.

Since u may not be in P , we rely neither on LCP clustering nor on RMQ. By using the SA we find candidate matches through binary search, identifying lower Algorithm 2: Algorithm for querying typing databases.

1 Input: An SA S for a set P of d profiles of length m each, concatenated as a string s of length md; a profile u of length m; an integer threshold 0 < k < m.

2 Output: The set X of distinct profiles that are at Hamming distance at most k from u, i.e., X = {v ∈ P | H(u, v) ≤ k}.

3 Initialization: Initialize a hash table HT to track verified profiles v and let L = ⌊ m k+1 ⌋. tion on the profiles, the expected value for the number of profiles matching u is no more than B σ L : we have B blocks in total and each block can only match at most one other block in u (since they must be aligned; line 9). Moreover, since we are not relying on the LCP array in this case, the verification step (line 13) takes O(m) time. Hence, the algorithm requires on average the following running time

O(m log md + m • B σ L ).
Let us analyze this further to obtain the relevant condition on k. We have the following:

m • B σ L = m • ⌊ m ⌊m/(k+1)⌋ ⌋ • d σ ⌊ m k+1 ⌋ ≤ m • ( m ⌊m/(k+1)⌋ ) • d σ m k+1 -1 ≤ m 2 d σ m k+1 -1 .
By some simple rearrangements we have that:

m 2 d σ m k+1 -1 = m 2 d (md) log σ log md ( m k+1 -1) = m(md) 1-(m-k-1) log σ (k+1) log md .
Consequently, in the case when the algorithm requires O(m log md) time on average. This algorithm was implemented and integrated in INNUENDO Platform, which is publicly available [START_REF]INNUENDO Platform[END_REF]. The INNUENDO Platform is an infrastructure that provides the required framework for data analyses from bacterial raw reads sequencing data quality insurance to the integration of epidemiological data and visualization.

k < (m -k -1) • log σ + (k + 1
As such, rapid methods for classification and search for closely related strains are a necessity for quick navigation through the platform database entries. More information about the project can be found at its website [START_REF] Innuendo | A novel cross-sectorial platform for the integration of genomics in surveillance of foodborne pathogens[END_REF].

As a starting point and for the purpose of this study, a subset of 2312 wgMLST profiles of Escherichia coli retrieved from Enterobase [START_REF] Sergean | EnteroBase[END_REF] were included in the IN-NUENDO database as well as their ancillary data and predefined core-genome cluster classification. Two tab-separated files containing the wgMLST and cgMLST profiles for the Escherichia coli strains were also created to allow storing information on the currently available profiles and for updating with profiles that will become available upon the platform analyses.

One of two index files are used depending on the type of search we want to perform: classification or search for k-closest. The cgMLST index file is used for strain classification, which relies on a nomenclature designed for the cgMLST profiles. As such, and since a pre-classification was performed on the database of Escherichia coli strains, we continued using it for comparison purposes. However, when searching for the k-closest profiles, we take into consideration all targets available in the wgMLST profiles using the wgMLST index file for a higher discriminatory power.

Each time a new profile is generated from the platform, it requires classification.

The INNUENDO Platform performs the classification step based on the approach described in our use case 2 with a given maximum of k differences over core genes. In the case of the search for the k-closest, it is useful to define the input data for visualization methods according to a defined number of differences on close strains. For each profile used as input for the search, the method searches for the k-closest strains considering at most k differences among all wgMLST loci. Since duplicate matches can occur between the profiles used for each search, the final file used as input for the visualization methods is the intersection of the results of the k-closest profiles between each input strain. The set of strain identifiers are then used to query the INNUENDO database to get the profiles and ancillary data to be sent to PHYLOViZ Online [START_REF] Ribeiro-Gonçalves | Phyloviz online: web-based tool for visualization, phylogenetic inference, analysis and sharing of minimum spanning trees[END_REF] for further analysis, namely with the goeBURST algorithm.

The drawback of using this method for classification and search is the need for rebuilding the index each time there is a new profile, which will depend on the number of profile entries on the database. Nevertheless, the number of updates is rather smaller compared to the number of queries and the index can be build in the background, with search functionalities still using the old index during the process. In our implementation, the index and related data structures are serialized in secondary memory and they are accessed by mapping them into memory. The implementation of the underlying tool is made publicly available [START_REF]Fast MLST searching and querying[END_REF].

The above described approaches in combination with the features offered by the INNUENDO Platform allow microbiologists to quickly and efficiently search for strains close to their strain of interest, allowing a more targeted, focused and simple visualization of results.

Experimental evaluation

We evaluated the proposed approach to compute the pairs of profiles at distance at most k using both real and synthetic datasets. We used real datasets obtained through different typing schemas, namely whole-genome multi-locus sequence typing (wgMLST) data, core-genome multi-locus sequence typing (cgMLST) data, and single-nucleotide polymorphism (SNP) data. Table 2 summarizes the real datasets used. We should note that wgMLST and cgMLST datasets contain sequences of integers, where each column corresponds to a locus and different values in the same column denote different alleles. Synthetic datasets comprise sets of binary sequences of variable length, uniformly sampled, allowing us to validate our theoretical findings.

We implemented both versions described above in the C programming language: one based on binary search over the SA; and another one based on finding clusters in the LCP array. Since allelic profiles can be either string of letters or sequences of integers, we relied on libdivsufsort library [START_REF] Mori | A lightweight suffix-sorting library[END_REF] and qsufsort code [START_REF] Larsson | Suffix sorting implementation to accompany the paper Faster Suffix Sorting[END_REF][START_REF] Larsson | Faster suffix sorting[END_REF], respectively. For RMQ over the LCP array, we implemented a fast well-known solution that uses constant time per query and linearithmic space for pre-processing [START_REF] Bender | The LCA problem revisited[END_REF].

All tests were conducted on a machine running Linux, with an Intel(R) Xeon(R) CPU E5-2630 v3 @ 2.40GHz (8 cores, cache 32KB/4096KB) and with 32GB of RAM. All binaries where produced using GCC 5.3 with full optimization enabled.

Synthetic datasets

We first present results with synthetic data for different values of d, m and k. All synthetic sequences are binary sequences uniformly sampled. Results presented in this section were averaged over ten runs and for five different sets of synthetic data.

The bound proved in Theorem 1 was verified in practice. For k satisfying the conditions in Theorem 1, the running time of our implementation grows almost linearly with n, the size of the input. We can observe in Fig. 1 a growth slightly above linear. Since we included the time for constructing the SA, the LCP array and the RMQ data structure, with the last one in linearithmic time, that was expected.

We also tested our method for values of k exceeding the bound shown in Theorem 1. For d = m = 4096 and a binary alphabet, the bound for k given in Theorem 1 is no more than ⌊m/(2 log m)⌋ = 170. For k above this bound we expect that proposed approaches are no longer competitive with the naïve approach. As shown in Fig. 2, for k > 250 and k > 270 respectively, both limits above the predicted bound, the running time for both computing pairwise distances by finding lower and higher bounds in the SA, and by processing LCP based clusters, becomes slower than the running time of the naïve approach.

In Fig. 3 we have the running time as a function of the number d of profiles, for different values of m and for k satisfying the bound given in Theorem 1. The running time for the naïve approach grows quadratically with d, while it grows linearly for both computing pairwise distances by finding lower and higher bounds in the SA, and by processing LCP based clusters. Hence, for synthetic data, as described by Theorem 1, the result holds.

Real datasets

For each dataset in Table 2, we ranged the threshold k accordingly and compared the approaches discussed in methods section with the naïve approach that computes the distance for all sequence pairs. Results are provided in Table 3.

In most cases, the approach based on the LCP clusters is the fastest up to two orders of magnitude compared to the naïve approach. As expected, in the case when data are not uniformly random, our method works reasonably well for smaller values of k than the ones implied by the bound in Theorem 1. As an example, the upper bound on k for C. jejuni would be around 200, but the running time for the naïve approach is already better for k = 64. We should note however that the number of candidate profile pairs at Hamming distance at most k is much higher than the expected number when data are uniformly random. This tells us that we can design a simple hybrid scheme that chooses a strategy (naïve or the proposed method)

depending on the nature of the input data. It seems also to point out clustering effects on profile dissimilarities, which we may exploit to improve our results. We leave both tasks as future work for the full version of this article.

We incorporated the approach based on finding lower and higher bounds in the SA in the implementation of goeBURST algorithm, discussed in methods section.

We did not incorporate the approach based on the LCP clusters as the running time did not improve much as observed above. Since running times are similar to those reported in Table 3, we discuss only the running time for C. jejuni. We need only to index the input once. We can then use the index in the different stages of the algorithm and for different values of k. In the particular case of goeBURST, we use the index twice: once for computing the number of neighbors at a given distance, used for untying links according to the total order discussed in the description of goeBURST algorithm in methods section, and a second time for enumerating pairs at distance below a given threshold. Note that the goeBURST algorithm does not aim to link all nodes, but to identify clonal complexes (or connected components)

for a given threshold on the distance among profiles [START_REF] Francisco | Global optimal eBURST analysis of multilocus typing data using a graphic matroid approach[END_REF]. In the case of C. jejuni dataset, and for k = 52, the running time is around 36 seconds, while the naïve approach takes around 115 seconds, yielding a three-fold speedup. In this case we get several connected components, i.e., several trees, connecting the most similar profiles. We provide the tree for the largest component in Fig. 4, where each node represents a profile. The nodes are colored according to one of the loci for which profiles in this cluster differ. Note that this tree is optimal with respect to the criterion used by the goeBURST algorithm, not being affected by the threshold on the distance. In fact, since this problem is a graphic matroid, the trees found for a given threshold will be always subtrees of the trees found for larger thresholds [START_REF] Papadimitriou | Combinatorial Optimization: Algorithms and Complexity[END_REF].

Comparing this tree with other inference methods is beyond the scope of this article;

the focus here was on the faster computation of an optimal tree under this model.

In many studies, the computation of trees based on pairwise distances below a given threshold, usually small compared with the total number of loci, combined with ancillary data, such as antibiotic resistance and host information, allows microbiologists to uncover evolution patterns and study the mechanisms underlying the transmission of infectious diseases [START_REF] Francisco | PHYLOViZ: phylogenetic inference and data visualization for sequence based typing methods[END_REF].

Conclusions

Most distance-based phylogenetic inference methods rely directly or indirectly on Hamming distance computations. The computation of a distance matrix is a common first step for such methods, taking Θ(md 2 ) time in general, with d being the number of sequences or profiles and m the length of each sequence or profile. For large-scale datasets this running time may be problematic; however, for some methods, we can avoid to compute all-pairs distances [START_REF] Pardi | Distance-based methods in phylogenetics[END_REF].

We addressed this problem when only a truncated distance matrix is needed, i.e., one needs to know only which pairs are at Hamming distance at most k.

This problem was motivated by the goeBURST algorithm [START_REF] Francisco | Global optimal eBURST analysis of multilocus typing data using a graphic matroid approach[END_REF], which relies on a truncated distance matrix by construction. Both the problem and techniques discussed here are related to average-case approximate string matching [START_REF] Fredriksson | Average-optimal single and multiple approximate string matching[END_REF][START_REF] Barton | Fast algorithms for approximate circular string matching[END_REF]. We proposed here an average-case linear-time and linear-space algorithm to compute the pairs of sequences or profiles that are at Hamming distance at most k, when

k < (m-k-1)•log σ log md
, where σ is the size of the alphabet. We integrated our solution in goeBURST demonstrating its effectiveness using both real and synthetic datasets.

We must note however that our analysis holds for uniformly random sequences and, hence, as observed with real data, the presented bound may be optimistic. It is thus interesting to investigate how to address this problem taking into account local conserved regions within sequences. Moreover, it might be interesting to consider in the analysis null models such as those used to evaluate the accuracy of distancebased phylogenetic inference methods [START_REF] Saitou | Introduction to Evolutionary Genomics[END_REF].

The proposed approach is particularly useful when one is interested in local phylogenies, i.e., local patterns of evolution, such as searching for similar sequences or profiles in large typing databases, as in our use case 2. In this case we do not need to construct full phylogenetic trees, with tens of thousands of taxa. We can focus our search on the most similar sequences or profiles, within a given threshold k. There are however some issues to be solved in this scenario, namely, dynamic updating of the data structures used in our algorithm. Note that after querying a database, if new sequences or profiles are identified, then we should be able to add them while keeping our data structures updated. Although more complex and dynamic data structures are known, a technique proposed recently for adding dynamism to otherwise static data structures can be useful to address this issue [START_REF] Munro | Dynamic data structures for document collections and graphs[END_REF]. This and other challenges raised above are left as future work. The current ability to rapidly sequence whole microbial genomes, has the promise to revolutionize these fields by allowing the identification of thousands of potentially clinically relevant targets in the genome. NGS data can be used to detect outbreaks in hospital settings or in the food industry, e.g., by monitoring the spread of antimicrobial resistance, an ever-growing concern. It can help also in the development of vaccines by helping, for instance, to determine targets conserved in the entire bacterial population.

However, it is becoming clear that the bottleneck shifted from the production of sequence data to its analysis. One of the major challenges is on how phylogenetic inference methods can be scaled up to analyze thousands of genetic loci in thousands of isolates. Usually, computing genetic evolutionary distances among a set of typing profiles or taxa dominates the running time of many of these methods.

It is important also to note that most of genetic evolution distance definitions rely, even if indirectly, on computing the pairwise Hamming distance among sequences or profiles.

In this work, we propose an average-case linear-time algorithm to compute pairwise Hamming distances among a set of taxa under a given Hamming distance threshold. This article includes both a theoretical analysis and extensive experimental results concerning the proposed algorithm. We further show how this algorithm can be successfully integrated into a well-known phylogenetic inference method, and how it can be used to speedup querying local phylogenetic patterns over large typing databases.

This work is an extension of the article "Towards distance-based phylogenetic inference in averagecase linear-time", presented at WABI 2017. It includes several revisions, including those raised by reviews of the workshop version of the paper, and it presents another real application of the developed methods, showing how it can be successfully used to speedup querying local phylogenetic patterns over large typing databases. The algorithm was integrated in INNUENDO platform, that is being developed under the INNUENDO project (http://www.innuendoweb.org/).

We wish to confirm that there are no known conflicts of interest associated with this publication.

The manuscript has been revised taking into account reviewers comments, and it has been read and approved by all named authors. Answers to reviewers' comments are also provided with our resubmission.

With my best regards,

Alexandre Francisco

Personal Cover Click here to download Personal Cover amob_cover_letter_r1.pdf

Theorem 1

 1 Given d profiles of length m each over an integer alphabet Σ of size σ > 1 with the letters of the profiles being independent and identically distributed random variables uniformly distributed over Σ, and the maximum Hamming distance 0 < k < m, Algorithm 1 runs in O(md) average-case time and space if k < (m -k -1) • log σ log md . Proof Let us denote by s the string of length md obtained after concatenating the d profiles. The time and space required for constructing the SA and the LCP arrays for s and the RMQ data structure over the LCP array is O(md).

18 ℓp := ℓ 19 else

 1819 if ℓp = -1 then 20 Pairs enumeration: 21 foreach Ct, with 0 ≤ t ≤ k do 22 foreach (p, q) ∈ Ct × Ct : p < q do 23 if (p, q) / ∈ HT then 24 HT := HT ∪ {(p, q)} 25

3

  and both indices correspond to the starting position of the ith block in their profiles.

Theorem 2

 2 and higher bounds in the SA, as discussed before. Hence, given the k + 1 nonoverlapping blocks of length L = ⌊ m k+1 ⌋ for u, we search for each one of them in O(L log md) time. Since we have k + 1 blocks, it takes O(kL log md) = O(m log md) time to search for all k + 1 blocks in u. Finally, we can then verify and report all candidate profiles v ∈ P as detailed in Algorithm 2. Although, in the worst case, Algorithm 2 runs in time O(md + m log md), as we may have d matches at most, we can prove a similar average case as in Theorem 1. Given a profile u and a set of d profiles of length m each, all over an integer alphabet Σ of size σ > 1, with the letters of the profiles being independent and identically distributed random variables uniformly distributed over Σ, the SA for the string s of length md obtained after concatenating the d profiles, and the maximum Hamming distance 0 < k < m, Algorithm 2 runs in O(m log md) average-case time if k < (m -k -1) • log σ + (k + 1) • log log md log md . Proof Let us denote by B the total number of blocks over s and by L the block length. We set L = ⌊ m k+1 ⌋ and thus we have that B = d⌊ m L ⌋. By the stated assump-

4 X 7 h 8 foreach ℓ ≤ j ≤ h do 9 if 14 if δ ≤ k then 15 X

 47891415 := ∅ 5 foreach 0 ≤ i ≤ k do 6 ℓ := LowerBinSearch(S, s, u[iL..(i + 1)L -1]) := HigherBinSearch(S, s, u[iL..(i + 1)L -1]) Aligned(j) = i then 10 v := ⌊S[j]/m⌋ 11 if v / ∈ HT then 12 HT := HT ∪ {v} 13 δ := H(v, u) := X ∪ {v} 16 Finalize: Return the set X.

  It uses the cgMLST index file for the search since the classification is constructed based on those number of loci. If the method returns at least one match, it classifies the new profile with the classification of the closest. If not, a new classification is assigned. A new entry is then added to the INNUENDO database as well as to the cgMLST and wgMLST profiles files and the index files are updated.

FiguresFigure 1

 1 Figures

Figure 2

 2 Figure 2 Synthetic datasets, with σ = 2 and m = 4096. Running time for computing pairwise distances by finding lower and higher bounds in the SA, and by processing LCP based clusters, as function of the number d of profiles and for different values of k.

Figure 3

 3 Figure 3 Synthetic datasets, with σ = 2 and k = ⌊m/(2 log m)⌋ according to Theorem 1.Running time for computing pairwise distances naïvely, by finding lower and higher bounds in the SA, and by processing LCP based clusters, as a function of the number d of profiles.

Figure 4

 4 Figure 4 The tree inferred for the largest connected component found with k = 52 for the C. jejuni dataset. Image produced by PHYLOViZ[START_REF] Nascimento | PHYLOViZ 2.0: providing scalable data integration and visualization for multiple phylogenetic inference methods[END_REF].

  Figure1Click here to download Figurertest.pdf

Table 1

 1 Data structures used in our approach for each step.

	1

Table 2

 2 Real datasets used in the experimental evaluation. (*)Dataset provided by the Molecular Microbiology and Infection Unit, IMM.

	Dataset	Typing method	Profile length	Number of distinct elements	Reference
	Campylobacter jejuni	wgMLST	5446	5669	(*)
	Salmonella enterica	wgMLST	3002	6861	[13]
	Salmonella typhi	SNP	22143	1534	[36]
	Streptococcus pneumoniae	cgMLST	235	1968	[37, 38, 39]

Table 3

 3 Time and percentage of pairs processed for each method and dataset.

	Dataset	k	t (s)	Naïve pairs (%)	Binary search t (s) pairs (%)	LCP clusters t (s) pairs (%)
		8	108.59		100	0.22	0.06	0.17	0.06
	C. jejuni	16 32	109.30 108.60		100 100	0.48 3.52	0.32 5.45	0.34 2.67	0.32 5.45
		64	108.60	100	231.05	99.98	162.36	99.98
		8	89.85		100	1.04	2.37	0.95	2.37
	S. enterica	16 32	87.26 85.36		100 100	7.16 36.29	12.69 33.22	6.73 30.76	12.69 33.22
		64	84.63		100	254.45	82.44	187.15	82.44
		89	28.83		100	16.63	91.48	12.02	91.48
	S. typhi	178	28.32		100	46.98	99.91	32.03	99.91
		890	30.04		100	113.57	100	129.14	100
		8	0.56		100	0.02	0.93	0.02	0.93
	S. pneumoniae	16 32	0.57 0.56		100 100	0.05 0.20	1.71 4.42	0.04 0.15	1.71 4.42
		64	0.58		100	5.63	73.36	5.01	73.36

Acknowledgements

This work was partly supported by the Royal Society International Exchanges Scheme, and by the following projects: BacGenTrack (TUBITAK/0004/2014) funded by FCT (Fundação para a Ciência e a Tecnologia) / Scientific and Technological Research Council of Turkey (Türkiye Bilimsel ve Teknolojik Araşrrma Kurumu, T ÜB ˙ITAK), PRECISE (LISBOA-01-0145-FEDER-016394) and ONEIDA (LISBOA-01-0145-FEDER-016417) projects co-funded by FEEI (Fundos Europeus Estruturais e de Investimento) from "Programa Operacional Regional Lisboa 2020" and by national funds from FCT, UID/CEC/500021/2013 funded by national funds from FCT, and INNUENDO project [25] co-funded by the European Food Safety Authority (EFSA), grant agreement GP/EFSA/AFSCO/2015/01/CT2 ("New approaches in identifying and characterizing microbial and chemical hazards"). The conclusions, findings, and opinions expressed in this review paper reflect only the view of the authors and not the official position of the European Food Safety Authority (EFSA).

Competing interests

The authors declare that they have no competing interests.

Author's contributions MC, APF, SPP and CV conceived the study and contributed for the design and analysis of the methods and experimental evaluation. APF, SPP and CV implemented Algorithm 1 and run the experiments. JAC conceived the case study 2 and contributed with the biological background. APF and BRG implemented Algorithm 2 and integrated it in INNUENDO Platform. All authors contributed to the writing of the manuscript. All authors have read and approved the final manuscript.

Author details

1 Faculdade de Medicina, Instituto de Microbiologia and Instituto de Medicina Molecular, Universidade de Lisboa, Lisboa, PT. 2 Department of Informatics, King's College London, London, UK. 3 INESC-ID Lisboa, Rua Alves Redol 9, 1000-029 Lisboa, PT. 4 Instituto Superior Técnico, Universidade de Lisboa, Lisboa, PT. 5 Instituto Superior de Engenharia de Lisboa, Instituto Politécnico de Lisboa, Lisboa, PT.