Method for the temporal calibration of a phylogenetic tree

Abstract

The invention relates to a method for the temporal calibration of a phylogenetic tree by means of digital processing, and comprises: - receiving, in digital format, a set of dated events involving the isolation or transmission of a pathogen, and a set of events to be dated involving the isolation or transmission of the pathogen, organised in the form of a phylogenetic tree following a traversal order from the root to the leaves. A plurality of events have a date to be determined; - forming 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, wherein the Hessian matrix H comprises what is referred to as a local block T defining the branches between the events on the phylogenetic tree. The block T has a size proportional to the number of events on the phylogenetic tree for which the date is to be determined; - identifying the dates to be determined by means of Newton's optimisation method.

Description

METHOD FOR TEMPORAL CALIBRATION OF A PHYLOGENETIC TREE

[0001] [The invention relates to the processing of epidemic data, and in particular to 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. Using these dates makes it possible, in particular, to determine the entry of a pathogen into a geographical unit (hospital ward, 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 epidemic's reproduction rate, known as 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 phylogenetic tree (genealogy) of the pathogens from these sequences, and estimating transmission dates from the phylogenetic tree.

[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. Because epidemics can involve a large number of cases, determining the dates of events must be feasible by numerical methods that potentially involve significant computing resources.

[0005] Known dating methods are based on models of evolution of RNA or DNA sequences, exploiting known dates from leaves of the phylogenetic tree. The sampling date is generally considered to correspond to each leaf. The dating of the nodes of the tree is referred to as temporal calibration. Temporal calibration is a crucial step in any phylodynamic analysis. Thus, its precision and speed of execution are key parameters of an epidemic surveillance system based on the analysis of the pathogen's genetic sequences. The goal is therefore to fix transmission events as precisely as possible, 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: computational time 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 short time, it is currently not possible to obtain sufficient accuracy and robustness for data processing.

[0007] The invention aims to solve one or more of these drawbacks. The invention thus relates to a method for temporal calibration of a phylogenetic tree by digital processing, as defined in the appended claim 1.

[0008] The invention also relates to variants of the dependent claims. Those skilled in the art will understand that each of the features of the dependent claims and the description can be independently combined with the features of an independent claim, without constituting an intermediate generalization.

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

[0010] 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:

[0011] Fig. 1 is a schematic illustration of an example of a phylogenetic tree;

[0012] Fig.2 illustrates the structure of a block of a Hessian matrix obtained for the phylogenetic tree of Figure 1;

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

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

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

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

[0017] Fig.7 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 first example sample;

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

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

[0020] 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 sample;

[0021] 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.

[0022] 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;

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

[0024] A method for temporally calibrating a phylogenetic tree by digital processing is provided. The method comprises receiving digital input data. The digital input data corresponds to a phylogenetic tree to be calibrated.

[0025] 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.

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

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

[0028] 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 branch is associated with an ascending event and a descending event in the phylogenetic tree. The length of the branches can be deduced from the dates of the events with which these branches are associated.

[0029] 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. We can also consider that leaves have unknown dates or on the contrary that nodes are associated with dates of known events.

[0030] 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.

[0031] 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.

[0032] 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.

[0033] 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.

[0034] 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.

[0035] 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.

[0036] 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 time and describes the overall rate of evolution in the phylogenetic tree. The dispersion of the evolutionary rate is a biologically meaningless quantity 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.

[0037] The time calibration process 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.

[0038] 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 depending on the unknown dates of the events. The logarithmic objective function is a likelihood function associated with the molecular clock model. It will be noted that in the context of the invention, an exact logarithmic objective function directly representing the molecular clock model is advantageously implemented, without introducing any approximation.

[0039] 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 (x1,...,xn), is defined as: Hij(f)=ô 2 f / ôxiôxj for i and j between 1 and n.

[0040] It should be noted that in the context of the invention, an analytical Hessian function is implemented. By not using an approximation of the Hessian function, the invention benefits from a theoretically quadratic convergence speed, which reduces the squared error at each iteration of a Newton optimization.

[0041] 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.

[0042] 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, as long as 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 connected neighbors 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, called reverse order. Such a permutation makes it possible to further reduce the computation time to implement the Newton optimization detailed later.

[0043] Newton's optimization methods involve the estimation and inversion of the Hessian matrix. The inventor also found that Newton's resolution methods theoretically presented a cubic computation time with the dimension of the matrix for usual matrices but presented a linear computation time with the dimension of this matrix if this matrix was sparse. The inventor having found that the Hessian matrix was 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 glance seem incompatible with the search for a reduced computation time, proved to be particularly appropriate for time calibration.A Newton optimization is an iterative method, each iteration adding a Newton direction s to a current vector v, with s and v such that H(F(v)) s = g(F(v)), with H a Hessian function.

[0044] 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, carrying out a temporal calibration becomes possible in a reasonable time with high accuracy, even for a large-scale epidemic. Moreover, this calibration is carried out on the basis of an exact method and not on the basis of an approximation.

[0045] 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.

[0046] Figures 1 and 2 illustrate, respectively, 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 T block illustrated in Figure 2, possibly non-zero values ​​are represented by a cross and always-zero values ​​by a zero. The row and column numbers in the T block reflect the position of the interior nodes of the phylogenetic tree.

[0047] 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.

[0048] At each iteration, the parameter vector is shifted in a direction given by the product of the inverse of the Hessian matrix and the gradient of the objective function. This shift 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 (e.g., with time calibration), a finite number of iterations is necessary to reach the objective.

[0049] 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 The parameter vector is composed of the dates to be estimated and the global parameters (if any). When parameters must obey constraints (e.g., positivity, inequality), 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 antiquity relative to the present, the present being typically defined as the sampling date of the last leaf. The dates are non-negative, as the date of an ancestor cannot be later than the date of its oldest descendant (or leaf). Similarly, the global parameters of the relaxed molecular clock model with additive variance are both non-negative, as a negative evolutionary rate or dispersion has no biological meaning.In this context, all parameters to be estimated are projected onto real numbers by a logarithmic link function.

[0050] The following notations can be used: -di is the genetic distance of a branch ending with an event of index i; -h is the duration of the branch ending with the event with index i; -ti is the age of the node or leaf corresponding to the event with index i; denoting p(i) the parent of i and c(i) the child of i, we therefore have = tp^-ti or l=t p -t to simplify the notation; -the rate of evolution p of the relaxed clock model with additive variance and a dispersion parameter w, equal to 0 for a strict model, 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 Epidemiology Phylogenies. Mol. Biol. Evol. 38, 307-317 (2021 )').

[0051] The function F can be a sum, over all branches of the tree, of the logarithms of the probability density of the genetic distance of the branch under a Gamma distribution. The genetic distance d, interpreted as the number of changes in a branch in the relaxed clock model with additive variance, follows a Gamma distribution with parameters a=p *l / (1 + co) which corresponds to the shape of the distribution; [3= 1 / (1 + co) which corresponds to the distribution rate, with I the time elapsed between the ascending event and the descending event of the branch, with p a rate of evolution.

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

[0053] The corresponding density (simplified by removing the index i) can be defined by: where T is the gamma function.

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

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

[0056] [Math. 2]

[0057] F = S =2 -1 log node or leaf with index i, excluding the root node.

[0058] Using a Newton-Raphson optimization process, we jointly maximize the ages of the nodes ti , ... , t n -i as well as the clock model parameters. The optimization method must guarantee the tree order constraint, namely t P (i)> tj>t C(i) for all events i.

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

[0060] 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.

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

[0062] We denote by x s the parameter vector at step s (including the model parameters, here after the dates of the n-1 nodes). Let Ax s = VA(x s)=d A / dxs the gradient of the logarithmic objective function and A 2 x s = V 2 A (x s )=d 2 A / dxs 2 the Hessian of the logarithmic objective function for x s . Each node has a date tj= e xi and the order constraint of the phylogenetic tree imposes that t P (i)>ti.

[0063] We can define the amplitude of an iteration y s and the iteration rule Xs+i = x s - y s ( HAS 2 Xs) -1 Ax s , with ( A 2 x s ) _1 Ax s Newton's direction.

[0064] The iteration is feasible (whose result has a finite objective function) if and only if, for all values ​​of i, t P (j), s +i >tj, s +i , that is to say that e xp(i) s+1 > e xi s+1 noting that the exponential function increases monotonically,

[0065] Vi>1 : 0< x P (i),s+i - Xi,s+1 = x P (i),s+i + Ys( Ax P (i),s)- Xi,s - Ys xi,s

[0066] Thus, an optimal iteration amplitude Ys ​​respects the following rule:

[0067] 0 <Ys*< min(x P (i),s- x iiS ) / ( Axi, s - Ax P (i), s )= Ys + for i>1. 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.

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

[0069] Newton optimization involves inverting the Hessian matrix. Inverting a Hessian can be unstable. To circumvent this problem, a Tikhonov regularization can be implemented for the Newton step, I being an identity matrix. P tends towards Newton's iteration before regularization when <t>tends towards 0 and towards the gradient of this iteration when <t>tends towards infinity. In case of premature termination of the execution of Newton's optimization due to numerical instability, the value of <t>.

[0070] The Hessian matrix can help 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 relatively low processing cost at each iteration.

[0071] 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

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

[0073] [Math. 3]

[0075] The P block appears when parameters of the molecular clock model are unknown prior to time calibration. The P block is called 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.

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

[0077] 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.

[0078] This block inversion can be expressed as follows in the more generic case illustrated above:

[0079] [Math. 4]

[0083] Matrix A is a k by k matrix whose inversion requires little computation time.

[0084] The Newton-Raphson direction can be expressed as follows:

[0085] [Math. 6]

[0087] Moreover, the convergence speed of the time calibration is quadratic: each iteration reduces the error to the square.

[0088] 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 character 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 negative definite at maximum likelihood.The combination of negative-definiteness, sparseness, and tree structure allows us to obtain the Wald variances, i.e., 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 ​​generally 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 block T is sparse.

[0089] 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.

[0090] The Schur 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.

[0091] 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 TreeTime: 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).

[0092] 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 p=10, dispersion œ=1, sampling of the most recent 1% of 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 p, 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.Figure 4 illustrates the estimated rate of evolution (on the left the invention, on the right Treedater), Figure 5 illustrates the percentage error tMRCA (on the left the invention, on the right Treedater), Figure 6 illustrates comparatively the percentage error on the dates for treedater (on the ordinate) and for the invention (on the abscissa), and Figure 7 illustrates according to another presentation the percentage error on the dates for the invention process and for treedater (on the left the invention, on the right Treedater).

[0093] 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 p=10, dispersion co=1O, sampling of the 1% of 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 p, 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.Figure 8 illustrates the estimated rate of evolution (on the left the invention, on the right Treedater), Figure 9 illustrates the percentage error tMRCA (on the left the invention, on the right Treedater), Figure 10 illustrates comparatively the percentage error on the dates for treedater (on the ordinate) and for the invention (on the abscissa), and Figure 11 also illustrates comparatively the percentage error on the dates (on the left the invention, on the right Treedater).

[0094] Figure 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 gap in execution time between the method of the invention and Treedater increases exponentially with the sample size.

[0095] Figure 13 illustrates the evolution of the calculation time of the temporal 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 temporal calibration for very large samples.

[0096] Figure 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 isolation events or determine this information from samples. The terminals may 301 to 303 receive or implement DNA or RNA sequencing of these isolation events.

[0097] The reconstruction of the phylogenetic tree 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.

[0098] Known dating methods are based on models of RNA or DNA sequence evolution, exploiting the known dates of the leaves of the phylogenetic tree. The isolation date is generally considered to correspond to each leaf. A traversal order 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.

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

[0100] The calibrated phylogenetic tree can then be used by epidemiologists to anticipate, manage, predict, or treat an outbreak of the pathogen. The dates of these transmission events usefully inform decision-making for epidemic control, for example, by identifying the date of entry of a pathogen into a geographical unit (hospital ward, city, country), for example, by identifying whether the estimated transmission date between two hosts coincides with a period of contact between these hosts, or by identifying peaks in the frequency of transmission events. The dates of transmission events can also be followed by a determination of the reproduction rate R o of the epidemic.< / t> < / t> < / t>

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, H a Hessian 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. A method of time calibration according to claim 1 , in which the Hessian matrix H is of the form: [Math. 7] ” - [c T- pl with P a so-called global block corresponding to parameters to be determined of 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 pathogenic agent.

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 a Gamma distribution having the parameters a= *1 / (1 + co) which corresponds to the shape of the distribution and [3= 1 / (1 + co) 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 of temporal calibration of a phylogenetic tree according to any one of the preceding claims, 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 Schur 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