Method for temporal calibration of a phylogenetic tree

Abstract

The invention relates to a method for temporal calibration of a phylogenetic tree by digital processing, and comprises: - receiving in digital format a set of events dated isolations or transmissions of a pathogenic agent and a set of events to be dated isolations or transmissions of the pathogenic agent, organized in the form of a phylogenetic tree according to an order of traversal from the root to the leaves. Several events have a date to be determined; - forming a Hessian matrix of a logarithmic objective function F of a molecular clock model of the pathogenic agent, defining the relationship between the genetic distances, the dates and the durations between the events, the Hessian matrix H comprising a so-called local block T defining the branches between said events of the phylogenetic tree. The block T has a size proportional to the number of events of the phylogenetic tree whose date is to be determined.;-identify the dates to be determined by a Newton optimization method. Figure to be published with the abstract: Fig. 1.

Description

Title of the invention: Method for temporal calibration of a phylogenetic tree

[0001] The invention relates to the processing of epidemic data, and in particular the processing of genetic sequences to extrapolate information concerning epidemic transmissions.

[0002] An infectious epidemic involves pathogen transmission events. The dates of these transmission events are a valuable tool to aid decision-making for epidemic control. The use of these dates makes it possible, in particular, to determine the entry of a pathogen into a geographical unit (hospital department, city, country) or to identify peaks in the frequency of transmission events. The dates of transmission events also make it possible to retrospectively estimate the reproduction rate of the epidemic, called R0.

[0003] The dates of transmission events are not known exactly. These dates can be estimated through phylodynamic analysis of pathogen samples. This phylodynamic analysis involves recovering samples from infected or colonized hosts, obtaining heritable characteristics of the pathogens, typically DNA or RNA sequences, reconstructing the phylogeny (genealogy) of the pathogens from these sequences, and estimating transmission dates from the phylogeny.

[0004] Phylodynamics is based on the creation of a phylogenetic tree whose leaves are samples of pathogens and whose nodes coincide with pathogen transmission events between hosts when these nodes connect leaves from distinct hosts. Each leaf represents a directly observed pathogen and each node represents a divergence event, not directly observed. Since epidemics can involve a large number of cases, determining the dates of events must be feasible by numerical methods potentially involving significant computing resources.

[0005] Known dating methods are based on models of evolution of RNA or DNA sequences, exploiting the known dates of the leaves of the phylogenetic tree. It is generally considered that the sampling date corresponds to each leaf. The dating of the nodes of the tree is designated by the term temporal calibration. Temporal calibration is a crucial step in any phylodynamic analysis. Thus, its precision and speed of execution are preponderant parameters of an epidemic surveillance system based on the analysis of genetic sequences of the pathogen. The aim is therefore to fix as precisely as possible transmission events, for example by quantifying the durations between two parent nodes.

[0006] Many temporal calibration methods have been developed. To achieve accurate and robust results, known temporal calibration methods have a common problem of computation time which increases faster than the number of samples in the phylogenetic tree to be analyzed (typically exponentially). For a large epidemic, it is currently not possible to perform a temporal calibration in a reasonable time. Conversely, in a reduced time, it is currently not possible to obtain sufficient precision and robustness for data processing.

[0007] The invention aims to resolve one or more of these drawbacks. The invention thus relates to a method for temporal calibration of a phylogenetic tree by digital processing, comprising: -receiving in digital format a set of events dated isolations or transmissions of a pathogenic agent and a set of events to be dated isolations or transmissions of the pathogenic agent, the events being organized in the form of a phylogenetic tree comprising a root, leaves, nodes ordered according to an order of traversal from the root to the leaves in the phylogenetic tree, the root, the leaves and the nodes each being associated with one of said events, the root, the leaves and the nodes being connected by branches having a length corresponding to a genetic distance between events, several of said events having a date to be determined, each of said branches being associated with an ascending event and a descending event of the phylogenetic tree according to said order; - formation of a Hessian matrix of a logarithmic objective function F of a molecular clock model of the pathogenic agent, defining the relationship between the genetic distances, the dates and the durations between the events, the Hessian matrix H comprising a so-called local block T defining the branches between said events of the phylogenetic tree, the block T having a size proportional to the number of events of the phylogenetic tree whose date is to be determined; -identify the dates to be determined in a solution vector w, by a Newton optimization method of solving the equation H(F(v)) s = g(F(v)), with v a current vector, with g a gradient function, and s an unknown or Newton direction, by fixing the solution vector w to the value of the current vector v when the variation of the current vector v between two iterations of the Newton method is less than a threshold or when the variation of the objective function F between two iterations of the Newton method is less than a threshold.

[0008] The invention also relates to the following variants. Those skilled in the art will understand that each of the characteristics of the following variants can be combined independently with the above characteristics, without constituting an intermediate generalization.

[0009] According to a variant, the Hessian matrix H is of the form:

[0010] [Math.l]

[0011] with P a so-called global block corresponding to parameters to be determined of the molecular clock model of the pathogenic agent;

[0012] with C and C blocks corresponding to the partial derivatives each involving an event and a global parameter of the molecular clock model of the pathogen.

[0013] According to another variant, at least one global parameter of the molecular clock model of the pathogenic agent is to be determined.

[0014] According to yet another variant, the resolution of the linear system in Newton optimization comprises a sparse LU decomposition or a sparse QR decomposition implemented on the block T.

[0015] According to a variant, the events associated with the leaves are dated and the events associated with the nodes are to be dated.

[0016] According to another variant, the length of each branch is defined in units of genetic distance or in number of genetic substitutions.

[0017] According to yet another variant, the molecular clock model of the pathogenic agent is a relaxed clock model with additive variance, comprising parameters defining a rate of evolution in units of genetic distance per unit of time and a dispersion of the rate of evolution representative of the variability of the speed of evolution between different branches of the phylogenetic tree.

[0018] According to another variant, said function F is of the form

[0019] a * log [3 - log T(a) +(a-1) * log d -[3* d with T a Gamma distribution having the parameters a=p *1 / (1+ œ) which corresponds to the shape of the distribution and [3= 1 / (1+ œ) which corresponds to the rate of the distribution, with co a dispersion parameter of the relaxed clock model with additive variance.

[0020] According to yet another variant, in said block T, a row index i corresponding to a single node and the same column index value j corresponding to this single node, the row indices i and column indices j being assigned in the reverse order of the traversal of the phylogenetic tree.

[0021] According to a variant, each node is associated with a single ascending event and at most two descending events.

[0022] According to another variant, the amplitude of displacement in the direction of the vector s during an iteration of Newton's resolution is defined to respect the rule that a descendant event is always later than its ascendant event in the phylogenetic tree.

[0023] According to another variant, the method comprises a calculation of a level of uncertainty of said solution by a Wald method, including a calculation of the diagonal of the inverse of the Hessian matrix H of the likelihood function applied to said vector w.

[0024] According to yet another variant, a Schurr complement method is used to obtain the diagonal of the inverse of the Hessian matrix H when calculating the uncertainty level of the solution.

[0025] The invention also relates to a digital system, comprising processing means configured to implement a time calibration method as defined above.

[0026] Other characteristics and advantages of the invention will emerge clearly from the description given below, for information purposes only and in no way limiting, with reference to the appended drawings, in which:

[0027] [Fig-1] is a schematic illustration of an example of a phylogenetic tree;

[0028] [Fig.2] illustrates the structure of a block of a Hessian matrix obtained for the tree phylogenetics of [Fig.l];

[0029] [Fig.3] illustrates an example of an infrastructure for implementing a calibration according to the invention;

[0030] [Fig.4] illustrates the estimated rate of evolution for a method according to the invention and a reference to the state of the art, with a first sample example;

[0031] [Fig.5] illustrates the percentage error tMRCA for the method according to the invention and the reference to the state of the art, with the first sample example;

[0032] [Fig.6] illustrates comparatively the percentage of error on the dates for the reference to the state of the art and for the invention, with the first sample example;

[0033] [Fig.7] illustrates comparatively the percentage of error on the dates for the process according to the invention and the reference of the state of the art, with the first example sample;

[0034] [Fig.8] illustrates the estimated rate of evolution for a method according to the invention and a reference to the state of the art, with a second sample example;

[0035] [Fig.9] illustrates the percentage error tMRCA for the method according to the invention and the reference to the state of the art, with the second sample example;

[0036] [Fig. 10] illustrates comparatively the percentage of error on the dates for the method according to the invention and the reference of the state of the art, with the second example of sample;

[0037] [Fig. 11] illustrates comparatively the percentage of error on the dates for the method according to the invention and the reference of the state of the art, with the second sample example.

[0038] [Fig. 12] illustrates the comparative execution time between the method of the invention and the reference method of the state of the art, as a function of the sample size, for a given digital processing system;

[0039] [Fig. 13] illustrates the evolution of the execution time of the method of the invention for large samples.

[0040] The invention relates to a method for temporal calibration of a phylogenetic tree by digital processing. The method comprises receiving digital input data. The digital input data corresponds to a phylogenetic tree to be calibrated.

[0041] The digital input data thus includes: -a set of dated events of isolation or transmission of a pathogen. The heritable characteristics of isolated pathogens can be obtained through DNA or RNA sequencing of the pathogen, or through other data sources such as Raman spectroscopy, the study of methylation profiles. Analyses to obtain the heritable characteristics of isolated pathogens can be carried out at a time lag compared to the isolations themselves. And - a set of events dating from isolation or transmission of the pathogenic agent.

[0042] The digital input data of the events are organized in the form of a phylogenetic tree comprising: a root, leaves, nodes ordered according to an order of traversal from the root to the leaves in the phylogenetic tree.

[0043] The purpose of the temporal calibration of the phylogenetic tree is to determine dates of isolation or transmission for events, these dates being initially undetermined.

[0044] The root, leaves and nodes are each associated with one of the events. The root, leaves and nodes are connected by links or branches. The branches have a length corresponding to a genetic distance between the events they connect. Each of the branches is associated with an ascending event and a descending event in the phylogenetic tree. The duration of the branches can be deduced from the dates of the events with which these branches are associated.

[0045] For the sake of simplification, the dates of the dated events correspond to the leaves, that is to say isolation events, the dates of the events to be dated correspond to transmission nodes. It is also possible to envisage that leaves have unknown dates or, on the contrary, that nodes are associated with dates of known events.

[0046] The branches of the phylogenetic tree can be represented by a set of triplets (identification of an ascending event, identification of a descending event, length of the branch), each triplet representing a branch between an ascending event and a descending event. The known length of this branch is for example defined in units of genetic distance, for example in number of substitutions or mutations.

[0047] The leaves of the tree typically represent the observed DNA or RNA sequences. A leaf has no descendants and therefore does not appear as an ancestor in the list of triplets. The tree also has a root (oldest node). The root has no ancestors and therefore does not appear as a descendant in the list of triplets.

[0048] In an earlier step of creating the phylogenetic tree, known per se from the prior art, the root is typically determined by an outgroup method or by a systematic analysis of all possible roots to select the root that maximizes the correlation between, on the one hand, the genetic distance between the root and the leaves and, on the other hand, the date of the leaves. This method is generally referred to as a root-to-leaf regression.

[0049] The relationship between genetic distances and durations between events or time intervals (also called evolutionary durations) between ancestors and descendants is described by a molecular clock model. Many molecular clock models are known per se. A known and robust molecular clock model is the relaxed clock model with additive variance. Such a molecular clock model will be advantageously implemented within the scope of the invention.

[0050] A molecular clock model typically includes one or more parameters that describe the relationship between genetic distances and time intervals at the scale of the phylogenetic tree as a whole. These parameters are referred to as global parameters.

[0051] The local parameters of the molecular clock model are the dates of the events (in particular the nodes) to be estimated. Each date is local to a particular event or node.

[0052] In the relaxed clock model with additive variance, two global parameters are defined: the evolutionary rate and the dispersion of the evolutionary rate. The evolutionary rate is expressed in units of genetic distance per unit of time and describes the overall rate of evolution in the phylogenetic tree. The dispersion of the evolutionary rate is a quantity without biological meaning and describes the variability of the evolutionary rate between branches of the tree. A dispersion parameter equal to zero defines a so-called strict clock model. In which the evolutionary rate is constant throughout the phylogenetic tree.

[0053] The time calibration method aims to estimate the most probable values ​​of the set of unknown dates and, jointly (if applicable), of the unknown global parameters of the molecular clock model.

[0054] The most probable values ​​of the dates of the events are obtained by maximizing the value of a logarithmic objective function associated with the molecular clock model of the pathogen dependent on the unknown dates of the events.

[0055] For the calibration of the phylogenetic tree, we form a Hessian matrix of this logarithmic objective function. As a reminder, a Hessian matrix H of a real-valued function F, applied to parameters (xl,...,xn), is defined as: Hij(f)=ô2f / ôxiôxj for i and j between 1 and n.

[0056] This Hessian matrix includes a so-called local block T. This block T includes the second partial derivatives of the objective function with respect to the dates of the events to be determined. The block T thus has a size proportional to the number of events in the phylogenetic tree whose date is to be determined.

[0057] In the context of the invention, the block T has the property of being sparse, that is to say that it necessarily contains zero values. Indeed, the second partial derivative of the objective function with respect to the dates of two nodes is zero, since these nodes are separated in the tree by at least one other node. In practice, each row or column of the block T contains at most 4 non-zero values ​​because a node has at most three neighbors connected in a binary tree. The number of non-zero values ​​in the block T therefore increases linearly with the number of dates to be determined. There is also a permutation of the block T such that the indices of its rows and columns (each corresponding to a node) respect the order of traversal of the tree from the leaves to the root. Such a permutation makes it possible to further reduce the computation time to implement the Newton optimization detailed later.

[0058] Newton's optimization methods involve the estimation and inversion of the Hessian matrix. The inventor has also found that Newton's resolution methods theoretically present a cubic computation time with the dimension of the matrix for usual matrices but present a linear computation time with the dimension of this matrix if this matrix is ​​sparse. The inventor having found that the Hessian matrix is ​​sparse in its application to a phylogenetic tree for a relaxed clock model with additive variance, it appears that the time of time calibration is here linear with the number of dates to be determined. Thus, Newton's optimization methods, which at first sight seem incompatible with the search for a reduced computation time, have proven to be particularly appropriate for time calibration.

[0059] The invention thus makes it possible to carry out a temporal calibration of a phylogenetic tree in a time which evolves linearly with the number of samples processed. Thus, it becomes possible to carry out a temporal calibration in a reasonable time with high precision, even for a large-scale epidemic.

[0060] The method of the invention is probabilistic and is based on a maximum likelihood estimation of the parameters of the evolution model. The optimization of the method according to the invention is carried out without approximation since it is based on the analytical expression of the set of partial derivatives of the relaxed molecular clock model with additive variance.

[0061] Figures 1 and 2 respectively illustrate an example of a simplified phylogenetic tree and a corresponding leaf-to-root permutation of an example of a Hessian matrix for this phylogenetic tree. The phylogenetic tree thus comprises nodes 111 to 119 and leaves 201 to 210. In the block T illustrated in [Fig.2], the possibly non-zero values ​​are represented by a cross and the always zero values ​​by a zero. The row and column numbers of the block T reflect the position of the interior nodes of the phylogenetic tree.

[0062] Newton-Raphson optimization involves iteratively searching for the value of a parameter vector that maximizes an objective function or, more formally, that cancels the derivative of that objective function.

[0063] At each iteration, the parameter vector is moved in a direction given by the product of the inverse of the Hessian matrix and the gradient of the objective function. This movement corresponds in practice to the exact solution of the maximum under the simplifying assumption that the objective function is a quadratic form. When this assumption is true, the Newton-Raphson method directly maximizes the objective in a single step. In the general case where this simplifying assumption is an approximation (for example with a time calibration), a finite number of iterations is necessary to achieve the objective.

[0064] At each iteration, the objective value is compared to the previous value and the algorithm stops when the objective value is virtually equal to the previous one, which signals that the process has reached a plateau. The objective function here is the log-likelihood of the molecular clock model and the parameter vector is composed of the dates to be estimated and the global parameters (if any). When parameters must obey constraints (positivity, inequality for example), they are typically replaced by a parameter projected onto the real numbers (from negative infinity to positive infinity) using a link function such as the logarithm function, which projects a non-negative value onto the real numbers. The dates are here expressed as seniority relative to the present, the present being typically defined as the sampling date of the last leaf.The dates are non-negative, the date of an ancestor cannot be later than the date of its oldest descendant (or leaf). Similarly, the global parameters of the molecular clock model relaxed to . additive variance are both non-negative, a negative rate of change or dispersion having no biological meaning. In this context, all parameters to be estimated are projected onto real numbers by a logarithmic link function.

[0065] The following notations may be used: -di is the genetic distance of a branch ending with an event of index i; -li is the duration of the branch ending with the event with index i; -f is the age of a leaf corresponding to the event with index i; noting p(i) the parent of i and c(i) the child of i, we therefore have 1; = tp^-fou l=tp-t to simplify the notation; -the evolution rate p of the relaxed clock model with additive variance and a parameter co a dispersion parameter (equal to 1 for a strict model in the following implementation, or equal to 0 for a strict model in the Didelot publication) according to the principle taught by the Didelot publication (Didelot, X., Siveroni, I. & Volz, EM 'Additive Uncorrelated Relaxed Clock Models for the Dating of Genomic Epi-demiology Phylogenies. Mol. Biol. Evol. 38, 307-317 (2021)').

[0066] The genetic distance d, interpreted as the number of changes in a branch in the relaxed clock model with additive variance, follows a Gamma T distribution with the parameters a=p *1 / (1+ co) which corresponds to the shape of the distribution; [3= 1 / (1+ co) which corresponds to the distribution rate.

[0067] This distribution has a mean of p *1, a variance of p *1 / (1+ co) and a coefficient of variation of V((l+ co) / p *1

[0068] The corresponding density (simplified by removing the index i) can be defined by: P(dll, p , co)= (|3“ / E(a)) * d“ 1 * e = d(p *l / (l+œ»-l * ed / (l+ œ) * (1+ (0)-( M *1 / (1+ «m / *! / (1+

[0069] We can then define an objective function or branch logarithmic probability, associated with the molecular clock model: L = a * log [3 - log T(a) +( a-1) * log d -[3* d =a * log (|3*d) - log T(a) - log d -[3* d = (p *1 / (1+ co)) * log (d / (1+ (o))-log T(p *1 / (1+ (o))- log d -d / (l+ a»)

[0070] We can then define an objective or logarithmic probability function of the tree by:

[0071] [Math.2] *■“*4 ■ p( z J 'J i+c*. - 1+ü) A -J&V -sÿ> # r[(( p( ,^

[0072] By a Newton-Raphson optimization method, we jointly maximize the ages of the nodes tb ..., Li as well as the parameters of the clock model. The optimization method must guarantee the order constraint of the tree, namely tP(i)>ti>tc(i) for the set of events i.

[0073] All parameters including node dates are strictly positive and their error is proportional, so we optimize their logarithm.

[0074] Newton-Raphson optimization proves advantageous on large phylogenetic trees, due to its rapid convergence in number of iterations, even if each iteration is in theory more computationally intensive.

[0075] Newton-Raphson optimization is constrained. The optimization algorithm must guarantee that the constraints are respected during the iterations, the logarithmic optimization function being undefined outside the feasibility interval.

[0076] We denote by xs the vector of parameters at step s (including the model parameters, here after the dates of the n-1 nodes). Let Axs= VA(xs)=d A / dxs be the gradient of the logarithmic objective function and A2xs= V2 A (xs)=d2 A / dxs2 the Hessian of the logarithmic objective function for xs. Each node has a date f = exi and the order constraint of the phylogenetic tree imposes that tp(i)>ti.

[0077] We can define the amplitude of an iteration 7, and the iteration rule xs+i = xs - ys( A2xJ 1 Ax„ with ( A2xJ 1 Ax, the Newton direction.

[0078] The iteration is feasible (whose result has a finite objective function) if and only if, for all values ​​of i, tp(i)jS+i>tijS+i, that is to say that exP®'s+1> exis+1 noting that the exponential function increases monotonically,

[0079] For all i>l: 0< xp(i),s+1- xi>s+1 = xp(i)>s+1 + ys( Axp(i),s)- xi>s - ys Axi>s

[0080] Thus, an optimal iteration amplitude ys respects the following rule:

[0081] 0 <ys*< min(xp(i)jS- xi> s) / ( Axi>s - Axp(i)>s )= ys+ for i>l. In this interval, the optimal iteration amplitude is the one that maximizes the objective function. We advantageously use Brent's method to determine the optimal amplitude that maximizes the objective function in this interval.

[0082] When the distance between two successive values ​​of the parameter vector xs is less than a certain threshold (for example 1 / 1,000,000), the iterations of Newton's optimization can be stopped.

[0083] Newton optimization involves the inversion of the Hessian matrix. The inversion of a Hessian can be unstable. To get around this problem, we can implement a Tikhonov regularization for the Newton step, with p= (0+1) * (A2x+O* tr(A2x). I )1 Ax I being an identity matrix. P tends towards the Newton iteration before regularization when O tends towards 0 and towards the gradient of this iteration when O tends towards infinity. In case of premature termination of the execution of the Newton optimization due to numerical instability, the value of O is increased.

[0084] The Hessian matrix can help to accelerate the convergence of the optimization and provides approximation information regarding the reliability index. The implemented Newton optimization has a linear complexity with the variable n due to the sparse nature of the Hessian matrix, which allows for rapid convergence at a fairly low processing cost at each iteration.

[0085] We can partition the gradient of the parameters into a vector t of the logarithms of the date gradients of the nodes, and into a vector of the gradient of the parameters linked to the model, i.e. Ax=[t,0]T

[0086] Which can be translated by the following form of the Hessian:

[0087] [Math.3] TC (n -1) *( n-1) H = ! CP , k*k

[0088] The P block appears when parameters of the molecular clock model are unknown prior to the time calibration. The P block is said to be global because it involves parameters common to all nodes of the tree. The size of the P block is defined by the chosen molecular clock model, and is independent of the size of the tree.

[0089] Blocks C and C' include partial derivatives involving a node and a global parameter.

[0090] Rather than performing the complete inversion of the Hessian matrix, one can also perform a block inversion of the matrix called the Schur complement method, by sparse LU decomposition, when the sparse matrix H respects the root to leaf permutation.

[0091] This block inversion can be expressed as follows in the more generic case illustrated above: [Math.4] H l = r^+T^c^e-cr^c'T 1 -[8-CT l C) CT

[0092]

[0093]

[0094]

[0095]

[0096]

[0097]

[0098] [Math.5] = +B AB . -AB -BA HAS . with ( B - CT l C y 1 = and B = CT 1 Matrix A is a k by k matrix whose inversion requires little computation time. The Newton-Raphson direction can be expressed as follows: [Math.6] _ 'Th-B At e-Bt)' A(e-Bt) Moreover, the convergence speed of the time calibration is quadratic: each iteration reduces the error to the square. Once the optimal values ​​(and possibly global parameters) have been obtained, an estimate of their uncertainty (e.g., in the form of a 95% confidence interval) can be obtained by the Wald method. This method relies on the asymptotically Gaussian nature of the distribution of the parameters estimated at maximum likelihood, the variance of the Gaussian distribution of each parameter being the diagonal entry of the inverse of the negative of the Hessian matrix. As with the Newton direction solution, the Wald method relies on the inversion of the Hessian matrix. Such an inversion is generally infeasible in practice for large problems. However, this matrix is ​​typically defined negative at maximum likelihood.The combination of negative-definiteness, sparseness, and tree structure allows us to obtain the Wald variances, the diagonal of the inverse of the negative of the Hessian matrix, very efficiently using the sparse Cholesky decomposition. Indeed, the inverse of a sparse matrix is ​​usually dense. But with a permutation from the root to the leaves, we can benefit from the fact that the inverse of the Cholesky decomposition of the T block is sparse. The diagonal entries (and only these) of the inverse of the Hessian matrix are obtained without wasting time on the off-diagonal entries, by multiplying only the necessary vectors of the inverse of the Cholesky decomposition. This entire procedure allows obtaining an uncertainty estimate for all parameters (dates and global parameters) with linear computational complexity, where a standard implementation with cubic complexity would be infeasible.

[0099] The Schurr complement method can be implemented to obtain the diagonal of the inverse of the Hessian matrix H when calculating the uncertainty level of the solution.

[0100] The performance of a time calibration method can be compared to the methods TREETIME (frequentist, approximate, sequence-based approach, described in Sagulenko, P., Puller, V. & Neher, RA TreeThne: Maximum-likelihood phylodynamic analysis. Virus Evol. 4, (2018)), and TREEDATER (frequentist, approximate, distance-based approach, described in Volz, EM & Frost, SDW Scalable relaxed clock phylogenetic dating. Virus Evol. 3, (2017).

[0101] As illustrated in the diagrams of Figures 4 to 7, the estimation error of a method according to the invention is illustrated with the so-called Treedater method, assumed to be the most accurate known method. The simulations were carried out on 1024 simulated phylogenies of 50 isolates presenting a good quality temporal signal (evolution rate q=10, dispersion co=1, sampling of the 1% of the most recent isolates). The estimation error of the method according to the invention is generally lower than that of Treedater, for all the parameters including the evolution rate q, the most recent common ancestor (tMRCA) and the dates of the nodes. The error on the dates rises to more than 60% for Treedater and does not exceed 30% for the method of the invention. [Fig.4] illustrates the estimated rate of evolution (on the left the invention, on the right Treedater), [Fig.5] illustrates the percentage of error tMRCA (on the left the invention, on the right Treedater), [Fig.6] illustrates comparatively the percentage of error on the dates for treedater (on the ordinate) and for the invention (on the abscissa), and [Fig.7] illustrates according to another presentation the percentage of error on the dates for the process of the invention and for treedater (on the left the invention, on the right Treedater).

[0102] As illustrated in the diagrams of Figures 8 to 11, the estimation error of a method according to the invention is illustrated with the so-called Treedater method. The simulations were carried out here on 1024 simulated phylogenies of 50 isolates presenting a temporal signal of poor quality (evolution rate q=10, dispersion co=10, sampling of the 1% of the most recent isolates). The estimation error of the method according to the invention is generally lower than that of Treedater, for all parameters including the evolution rate q, the most recent common ancestor (tMRCA) and the dates of the nodes. The maximum date error rises to 91% for the invention and 7410% for Treedater. The maximum tMRCA error is 179% for the invention and 25,380% for Treedater. [Fig.8] illustrates the estimated rate of evolution (on the left the invention, on the right Treedater), [Fig.9] illustrates the percentage of error tMRCA (on the left the invention, on the right Treedater), [Fig.10] illustrates comparatively the percentage of error on the dates for treedater (on the ordinate) and for the invention (on the . abscissa), and [Fig.l 1] also illustrates comparatively the percentage of error on the dates (on the left the invention, on the right Treedater).

[0103] [Fig. 12] illustrates (with a logarithmic scale) the comparative execution time between the method of the invention and treedater, as a function of the sample size. The results for Treedater correspond to the upper envelope, the results for the invention correspond to the lower envelope. It can be seen that the difference in execution time between the method of the invention and Treedater increases cubically with the sample size, as anticipated.

[0104] [Fig. 13] illustrates the evolution of the calculation time of the time calibration method according to the invention as a function of the sample size. As anticipated by the theory, the calculation time evolves substantially linearly with the sample size and therefore makes it possible to envisage a time calibration for very large samples.

[0105] [Fig. 3] schematically illustrates an example of infrastructure 1 for implementing the invention. The infrastructure is based on terminals 301 to 303 configured to retrieve dates and information on isolation events, typically corresponding to the leaves of the phylogenetic tree. The terminals 301 to 303 retrieve information on the isolation events or determine this information from samples. The terminals can 301 to 303 receive or implement DNA or RNA sequencing of these isolation events.

[0106] The reconstruction of the phylogeny can be carried out by means of dedicated processing means or on a digital processing system 320 also intended for the temporal calibration of the phylogenetic tree. For this purpose, the information coming from the terminals 301 to 303 can be collected by the appropriate processing system by means of a communication network 310. The processing system creates in a manner known per se a phylogenetic tree whose leaves are the events of isolation of the pathogen and whose nodes coincide with the events of transmission of the pathogen between hosts when these nodes connect leaves from distinct hosts. Each leaf represents a directly observed pathogen and each node represents a divergence event, generally not directly observed.

[0107] Known dating methods are based on models of evolution of RNA or DNA sequences, exploiting the known dates of the leaves of the phylogenetic tree. It is generally considered that the date of isolation corresponds to each leaf. An order of traversal from the root to the leaves is defined in the phylogenetic tree. The branches of the phylogenetic tree have a length corresponding to a genetic distance calculated by the appropriate processing system.

[0108] The processing system 320 recovers the generated phylogenetic tree. The system processing 320 implements the temporal calibration method as described previously, on the basis of this phylogenetic tree. The solution determined by this temporal calibration method is then stored in a database 330.

Claims

Claims

1. Method for temporal calibration of a phylogenetic tree by digital processing, comprising: - receiving in digital format a set of events dated isolations or transmissions of a pathogenic agent and a set of events to be dated isolations or transmissions of the pathogenic agent, the events being organized in the form of a phylogenetic tree comprising a root, leaves, nodes ordered according to an order of traversal from the root to the leaves in the phylogenetic tree, the root, the leaves and the nodes each being associated with one of said events, the root, the leaves and the nodes being connected by branches having a length corresponding to a genetic distance between events, several of said events having a date to be determined, each of said branches being associated with an ascending event and a descending event of the phylogenetic tree according to said order;- formation of a Hessian matrix of a logarithmic objective function F of a molecular clock model of the pathogen, defining the relationship between the genetic distances, the dates and the durations between the events, the Hessian matrix H comprising a so-called local block T defining the branches between said events of the phylogenetic tree, the block T having a size proportional to the number of events of the phylogenetic tree whose date is to be determined;- identify the dates to be determined in a solution vector w, by a Newton optimization method of solving the equation H(F(v)) s = g(F(v)), with v a current vector, with g a gradient function, and s an unknown or Newton direction, by fixing the solution vector w to the value of the current vector v when the variation of the current vector v between two iterations of the Newton method is less than a threshold or when the variation of the objective function F between two iterations of the Newton method is less than a threshold.;

2. Time calibration method according to claim 1, in which the Hessian matrix H is of the form: [Math. 7] [T Cl with P a so-called global block corresponding to parameters £7 = LC PJ to be determined from the molecular clock model of the pathogenic agent; with C and C' blocks corresponding to the partial derivatives each involving an event and a global parameter of the molecular clock model of the pathogen.

3. Method for temporal calibration of a phylogenetic tree according to claim 2, in which at least one global parameter of the molecular clock model of the pathogenic agent is to be determined.

4. A method for temporal calibration of a phylogenetic tree according to any preceding claim, wherein the solving of the linear system in the Newton optimization comprises a sparse LU decomposition or a sparse QR decomposition implemented on the block T.

5. A time calibration method according to any preceding claim, wherein the events associated with the leaves are dated and the events associated with the nodes are to be dated.

6. A method of temporal calibration of a phylogenetic tree by digital processing according to any one of the preceding claims, wherein the length of each branch is defined in units of genetic distance or in number of genetic substitutions.

7. A method for temporal calibration of a phylogenetic tree according to any one of the preceding claims, wherein the molecular clock model of the pathogen is a relaxed clock model with additive variance, comprising parameters defining a rate of evolution in units of genetic distance per unit of time and a dispersion of the rate of evolution representative of the variability of the rate of evolution between different branches of the phylogenetic tree.

8. A time calibration method according to claim 5, wherein said function F is of the form a * log [3 - log T(a) +( a-1) * log d -[3* d with T a Gamma distribution having the parameters a=p *1 / (1+ œ) which corresponds to the shape of the distribution and [3= 1 / (1+ œ) which corresponds to the rate of the distribution, with co a dispersion parameter of the relaxed clock model with additive variance.

9. A method for temporal calibration of a phylogenetic tree according to any one of the preceding claims, wherein in said block T, a row index i corresponding to a single node and the same column index value j corresponding to this single node, the row indices i and column indices j being assigned in the reverse order of traversal of the phylogenetic tree.

10. A method of temporal calibration of a phylogenetic tree according to any one of the preceding claims, wherein each node is associated with a single ascending event and at most two descending events.

11. A method for temporal calibration of a phylogenetic tree according to any preceding claim, wherein the displacement amplitude in the direction of the vector s during an iteration of the Newton resolution is defined to respect the rule that a descendant event is always later than its ascendant event in the phylogenetic tree.

12. Method for temporal calibration of a phylogenetic tree according to any one of the preceding claims, comprising a calculation of a level of uncertainty of said solution by a Wald method, including a calculation of the diagonal of the inverse of the Hessian matrix H of the likelihood function applied to said vector w.

13. A method for temporal calibration of a phylogenetic tree according to claim 12, wherein a Schurr complement method is used to obtain the diagonal of the inverse of the Hessian matrix H when calculating the uncertainty level of the solution.

14. Digital system, characterized in that it comprises processing means configured to implement a time calibration method according to any one of the preceding claims.

Images