the evolutionary and demographic history of a measurable evolving viral population1 can be inferred by phylodynamic analysis of longitudinally sampled sequences.2 in particular , the shape of a phylogenetic tree often reflects the characteristics and the interactions among evolutionary and ecological processes .
for example , phylogenies of viruses under continuous positive selection , such as inter - host influenza or intrahost human immunodeficiency virus ( hiv ) , usually exhibit a marked staircase - like topology.3,4 pathogens under weak or absent positive selection show a more balanced tree shape , as it happens within the measles serotypes.5 trees displaying a star - like topology or increasing root to tip distance with sampling time are typical of an exponentially growing population , whilst the opposite pattern can usually be associated to constant or decreasing population sizes , as for example in dengue virus inter - host phylogenies.6 phylogenetic tree shapes can also be coupled with phenotypic traits ( either numeric or categorical ) or geographic correlates ( such as geographic origin of sampled strains ) , which can be analyzed via comparative methods in terms of evolutionary or spatiotemporal dynamics.7 several statistical approaches exist for determining if and how serially sampled sequences evolve under a strict or relaxed molecular clock.8,9 however , such methods do not give indications about the phylogenetic tree shape ( eg , staircase- or star - like ) and the related temporal structure ( eg , if sequences sampled at the same time point tend to cluster together and to be direct ancestors of sequences sampled at later time points ) .
some of these measures consider both topology and branch lengths of a tree , as well as phenotypic tip traits,1012 while others evaluate only the tree topology in relation to geographic or phenotypic characters associated with the sampled strains.1315 finally , there are purely topological measures based on tree symmetry / balance.1619 the temporal clustering ( tc ) is a recently developed statistic , which takes into account phylogenetic tree topology and sampling time of the tips.20 the tc statistic assesses the temporal structure of a phylogenetic tree , by evaluating the order of time changes from internal nodes to tips .
it is based on the maximum parsimony reconstruction of ancestral characters implemented in phylogeography,13 but it has been modified to prevent the estimation of temporally impossible state changes in tip dated trees ( ie , an earlier time point emerging from an ancestor assigned to a later time point )
. it also allows the comparison of phylogenies inferred from data sets including different number of times points and/or sampled sequences per time point . currently , there is no available software implementation of the tc statistic , although it can be calculated with a series of manual steps using the macclade program ( http://macclade.org/ ) . in this work ,
phylotempo is written in r , a free software environment for statistical computing and graphics ( http://www.r-project.org/ ) .
along with the tc implementation , several other tree topological measures were integrated in phylotempo using pre - existing r libraries , in a user - friendly graphical framework .
the program was tested on several longitudinally sampled intra - host hiv and simian immunodeficiency virus ( siv ) population data sets .
the results showed how the comparison of the tc statistic with other topological measures can provide multifaceted insights on the dynamic processes shaping the evolution of pathogenic viruses .
the original formulation of the tc statistic20 requires a phylogenetic tree with n taxa sampled at t different discrete time points .
a state ( ie , time ) transition matrix is then defined , where the cost of going from later to earlier time points is infinite ( ie , time - irreversibility ) .
the other costs are usually defined as integer linearly increasing with the time points ordering .
ancestral tree states are inferred using fitch s parsimony algorithm.21 a non - normalized tc score is then calculated by summing all the state changes across the tree branches according to the cost matrix weights . a tree with a perfect temporal structure , ie ,
a tree in which all tips sampled at time point ti are monophyletic and directly emerge from time point ti1 , would have a non - normalized tc equal to t1 .
conversely , a tree with the least temporal structure would have a maximal non - normalized tc , equal to the sum of the n number of taxa multiplied by the corresponding w weights of the cost matrix , ie , 
 max=j=1tnjw1j .
the normalized tc rescales the non - normalized tc value in the interval , by considering a background distribution of tc statistics obtained by shuffling the time points associated to the tree tips ( keeping the topology fixed ) and reestimating the ancestral characters .
specifically , the normalized tc statistic is where the savg and smax are -respectively- the average and the maximum non - normalized tc values observed in the randomized trees , while min is the minimum theoretical non - normalized tc , equal to t1 .
the numerator represents the deviation from the null hypothesis ( ie , no temporal clustering ) , while the denominator represents the range of possible values for the given number of taxa and time points .
coupled with the tc statistic , phylotempo includes also the following tree topology measures and tests of hypothesis : aldous graphical test and lilkelihood ratio test to decide if tree fit the yule or the uniform models;22 colless and sackin s shape statistics , both under the yule or uniform hypotheses;23 cherry count;24 pybus gamma.25 in addition , a simple tree statistic called  staircase - ness  is introduced , counting the proportion of sub - trees that are imbalanced ( ie , sub - trees where the left child contains more leaves than the right child , or vice - versa ) compared against the distribution of such proportions obtained from random trees . see the supplementary material for the properties of this measure .
all the code has been written in the r language . besides the standard core library set of r , the following r libraries have been used ( including their dependencies ) :  ape  ,  ade4  ,  phybase  ,  phylobase  ,  phangorn  ,  doby  ,  infotheo  ,  aptreeshape  ,  diversitree  ( http://www.r-phylo.org ) . the required input of phylotempo is a phylogenetic tree file in  newick  format and a two - column text file in which each tip name present in the phylogenetic tree is associated with its corresponding time of sampling ( a numeric value such as days or years ) .
if present , negative branch lengths are set to zero and then all branch lengths are added a 10 value .
the tree is rooted on the tip that gives the highest linear correlation between the root - to - tip distance and the sampling time of the tip , and finally it is ladderized .
the vector of sampling times is then discretized into time intervals by using an equal - frequency binning , where the optimal number of bins is the square root of the vector size .
the maximum allowed number of discrete time intervals is nine , and each time bin needs to contain at least two tips .
the tc statistic calculation is made upon the previous theoretical description . however , in this new implementation the ancestral characters are estimated using maximum likelihood26 rather than parsimony .
a major advantage of maximum likelihood is that it also allows for an optimized estimate of the weights of the transition cost matrix .
the number of tip randomizations is set to 300 by default , but the value can be modified by the user .
both graphical and text output are produced , where figures are plotted in multiple windows , text is printed in the r command - line window and results are saved in a tabulated file .
the graphical plots include : the phylogenetic tree with ancestral character state probabilities drawn with pie charts at internal nodes ; the tc statistic compared versus the randomized background distribution ; a linear correlation plot between the sampling times of tips and root - to - tip distances ; a kruskal - wallis test comparing distribution of root - to - tip distances with the discretized time points ; the staircase - ness , aldous , sackin s and cherry count statistics with the corresponding background randomizations .
the text output reports the aforementioned results as well as the p - values from the statistical tests . also ,
a script that allows the analysis of multiple trees in  nexus  format , from an a posteriori , eg , trees as output from mrbayes,27 or bootstrap analysis has been made available .
phylotempo is distributed under the gnu general public license and is available at https://sourceforge.net/projects/phylotempo/ for download .
the first data set included intra - host hiv-1 phylogenetic trees , inferred from serially sampled envelope ( env ) c2-v5 sequences , from nine untreated subjects with fast disease progression,4 named as the  shankarappa  data set after the first author of the paper .
the second data set included intra - host siv trees , inferred from env gp120 sequences sampled longitudinally from four experimentally infected rhesus macaques that were cd8-depleted before infection and progressed to aids within 75118 days post infection.28 the third data set included intrahost hiv-1 gag p24 trees from six untreated subjects enrolled in the options cohort29 all carrying the hla - b*5701 allele strongly associated with slower disease progression , that were followed longitudinally from early infection up to seven years.30 in table 1 the text output of phylotempo is reported , after running the program on each proof - of - concept data set ( note that for simplicity not all indicators output by phylotempo are shown ) . when comparing the former calculation of the tc statistic based on parsimony with the new one based on maximum likelihood , in general we found a high degree of linear  correlation (   0.8 , combining the three data sets , data not shown ) .
on average , the tc exhibited a weak linear  correlation with any of the other tree topology measures implemented in phylotempo ( average  = 0.10 , standard deviation 0.17 ) , including also dn / ds values ( estimated via the nei - gojobori method averaging across all positions ) .
the maximum value obtained was  = 0.42 , found with respect to the root - to - tip - distance vs. sampling time correlation .
the average tc statistic for the shankarappa data set was 0.29 ( st.dev 0.10 ) , for the siv data set was 0.11 ( st.dev 0.09 ) .
the options data set , based on the highly conserved hiv-1 gag p24 gene , allowed us to evaluate the effect of including or excluding identical sequences and resulted in a tc value of 0.27 ( st.dev 0.11 ) and 0.36 ( st.dev 0.19 ) when analyzing all sequences or only unique sequences , respectively .
figure 1 shows two of the trees analyzed ( options and the siv data sets , respectively ) , and includes the maximum likelihood estimate of the ancestral time states , with state probabilities reported as pie charts at each internal node .
figure 2 reports the placement of the tc statistic , as well as all the other tests , with respect to background random distributions or null hypotheses .
in addition , the correlation plot between the root - to - tip distances and the sampling time is shown , along with the boxplots of the root - to - tip distances stratified by the time intervals .
on average , the running time of phylotempo on an input phylogenetic tree with 100150 leaves ( 34 time points and 300 randomizations ) takes less than 5 minutes using a standard desktop computer . running times for trees of 300 or 400 tips increase to half or one hour .
in this paper we presented phylotempo , a set of scripts in the r language that calculates the tc clustering statistics and other measures of phylogenetic tree shape , with a comprehensive text and graphical output .
the choice of the r software environment gives to the tool the advantage to be available for many platforms ( microsoft windows , mac , or linux ) and , since r features a plethora of libraries both for phylogenetic analysis and graphics , to be ready for the inclusion of other functions related to the analysis of phylogenetic tree shape and comparative statistics .
although other programs that calculate tree shape statistics are available , such as the java application treestat ( http://tree.bio.ed.ac.uk/software/treestat/ ) , and path - o - gen ( http://tree.bio.ed.ac.uk/software/pathogen/ ) , as well as several command - line functions in r , this is the first that implements the tc statistic merging in a user - friendly interface both graphical and text outputs . in addition , phylotempo is capable of generating an a posteriori tc statistic , reading a tree ensemble in  nexus  format , such as the output by mrbayes ( http://mrbayes.sourceforge.net/ ) . as a future perspective in the context of bayesian analysis , a theoretical approach to derive an analytical formula for the tc statistic is advisable , allowing the avoidance of the time - consuming tree randomization for each tree .
first , tc does not correlate with any previously described topological tree measure implying that the new statistic evaluates aspects of the evolutionary process not captured by other methods .
second , tc does not correlate with estimated dn / ds ratios in different data sets .
several studies have interpreted temporally structured phylogenies as evidence of sequential viral population bottlenecks driven by continuous selection pressure.2,4,31,32 the trees inferred from the options data sets include hiv-1 sequences from patients with the hla - b*5701 allele that has been associated with slower disease progression , possibly due to strong positive or purifying selection driving viral escape from cytotoxic t lymphocyte recognition.33 interestingly , the tc calculated for the options data sets are not significantly different ( p = 0.21 from a t - test ) from those calculated for the shankarappa data sets .
the finding suggests that temporally structured genealogies may reflect intra - host evolutionary processes that are similar in two groups of patients characterized by different rates of disease progression and that may not be related to selection pressure .
however , it is important to point out that the subjects in the shankarappa data set were followed for a shorter period of time than the options subjects , and that the intervals between longitudinal samples were overall shorter in the first data set ( data not shown ) .
the low tc values for the shankarappa data set may simply reflect an incomplete turnover of the viral quasispecies , which can require up to 22 months,34 causing an intermix of sequences sampled at different time points .
moreover , archival viral strains expressed in cellular reservoirs would decrease the temporal structure of serially sampled genealogies , because sequences from later time points may share a most recent common ancestor with sequences collected much earlier in infection.35 therefore , the tc statistic could be a powerful tool to investigate the extent and impact of latent viral reservoirs in intra - host hiv-1 evolution .
finally , it is interesting to note that the siv data sets show the lowest tc values .
this may be the result of the relatively short time of infection in these animals , as well as a consequence of the depletion of cd8 t cells right after infection.28 in conclusion , the present work describes a practical and user - friendly implementation of a novel statistic to evaluate the shape of phylogenetic trees , inferred from longitudinal samples of measurably evolving viral populations , which can provide significant insights on underlying evolutionary processes linked to infection dynamics and pathogenesis .
the staircase - ness measure counts the ( i ) proportion of sub - trees that are imbalanced ( ie , sub - trees where the left child contains more leaves than the right child , or vice - versa ) . an alternative formulation ( ii )
is to make the average of all the min(l , r)/max(l , r ) values of each sub - tree , where l and r are the number of leaves in the left and right children of a subtree . in this work
we compared the staircase - ness values against the distribution of such proportions obtained from random trees .
however , there are also a few properties of this measure that are worth to be analyzed analytically .
first of all , the staircase - ness of perfectly balanced binary trees is always zero , whichever formulation is used . on the other hand ,
the staircase - ness of perfectly imbalanced trees ( ie , ladder - like trees ) is always one when counting the proportions ( ie , formulation i ) , whilst depends on the number of leaves when performing the average ( ie , formulation ii ) .
specifically , the staircase - ness values for perfectly imbalanced trees using formulation ( ii ) is 
 sn=1ni=1n1i , where n is the number of subtrees .
this formula tends to zero as n increases since the limit of the series 
 limnsn=1ni=1n1i converges to zero , as previously demonstrated by e. cesro ( http://en.wikipedia.org/wiki/ces%c3%a0ro_mean .
the distribution of the staircase - ness values obtained by simulating random trees ( function rtree(number_of_tips , branch_length = runif(1 ) ) of the r library  ape  ) does not pass the shapiro - wilk normality test ( p<<0.0001 , even by considering only trees with a number of tips > 300 ) , neither resembles a gamma distribution , whose parameters had been fit on the actual data ( p<<0.0001 , using a kruskal / wallis test on simulations ) .
however , the average values of both definitions values look stable across all the tree sizes ( fig .
s1 ) , while the standard deviation seems to decrease by increasing the tree size .
the limits of the average staircase - ness values for formulation ( i ) and ( ii ) are close to 0.61 and 0.64 , respectively .
distribution of staircase - ness values from random trees ( 5,000 simulations ) , by varying tree size ( from 3 to 2,000 leaves ) .
notes : upper panels show results for formulation ( i ) , whilst lower panels for formulation ( ii ) . left panels represent the scatterplot of all staircase - ness values depending on the tree size , with a global average and standard deviation indicated in red .
central panels show the boxplots of staircase - ness values by stratifying for tree sizes ( 5 equal - width intervals spanning tree sizes between 3 and 2,000 ) , with the corresponding stratified average and standard deviation .
the right panels show the histograms for all staircase - ness values , and compares them with simulated distributions whose parameters have been fit on the empirical data ( gaussian and gamma functions ) .
the staircase - ness measure counts the ( i ) proportion of sub - trees that are imbalanced ( ie , sub - trees where the left child contains more leaves than the right child , or vice - versa ) . an alternative formulation ( ii )
is to make the average of all the min(l , r)/max(l , r ) values of each sub - tree , where l and r are the number of leaves in the left and right children of a subtree . in this work
we compared the staircase - ness values against the distribution of such proportions obtained from random trees .
however , there are also a few properties of this measure that are worth to be analyzed analytically .
first of all , the staircase - ness of perfectly balanced binary trees is always zero , whichever formulation is used . on the other hand ,
the staircase - ness of perfectly imbalanced trees ( ie , ladder - like trees ) is always one when counting the proportions ( ie , formulation i ) , whilst depends on the number of leaves when performing the average ( ie , formulation ii ) .
specifically , the staircase - ness values for perfectly imbalanced trees using formulation ( ii ) is 
 sn=1ni=1n1i , where n is the number of subtrees .
this formula tends to zero as n increases since the limit of the series 
 limnsn=1ni=1n1i converges to zero , as previously demonstrated by e. cesro ( http://en.wikipedia.org/wiki/ces%c3%a0ro_mean .
the distribution of the staircase - ness values obtained by simulating random trees ( function rtree(number_of_tips , branch_length = runif(1 ) ) of the r library  ape  ) does not pass the shapiro - wilk normality test ( p<<0.0001 , even by considering only trees with a number of tips > 300 ) , neither resembles a gamma distribution , whose parameters had been fit on the actual data ( p<<0.0001 , using a kruskal / wallis test on simulations ) .
however , the average values of both definitions values look stable across all the tree sizes ( fig .
s1 ) , while the standard deviation seems to decrease by increasing the tree size .
the limits of the average staircase - ness values for formulation ( i ) and ( ii ) are close to 0.61 and 0.64 , respectively .
distribution of staircase - ness values from random trees ( 5,000 simulations ) , by varying tree size ( from 3 to 2,000 leaves ) .
notes : upper panels show results for formulation ( i ) , whilst lower panels for formulation ( ii ) . left panels represent the scatterplot of all staircase - ness values depending on the tree size , with a global average and standard deviation indicated in red .
central panels show the boxplots of staircase - ness values by stratifying for tree sizes ( 5 equal - width intervals spanning tree sizes between 3 and 2,000 ) , with the corresponding stratified average and standard deviation .
the right panels show the histograms for all staircase - ness values , and compares them with simulated distributions whose parameters have been fit on the empirical data ( gaussian and gamma functions ) .
 serially - sampled nucleotide sequences can be used to infer demographic history of evolving viral populations . 
 the shape of a phylogenetic tree often reflects the interplay between evolutionary and ecological processes . 
 several approaches exist to analyze the topology and traits of a phylogenetic tree , by means of tree balance , branching patterns and comparative properties . the temporal clustering ( tc ) statistic is a new topological measure , based on ancestral character reconstruction , which characterizes the temporal structure of a phylogeny . here 
 , phylotempo is the first implementation of the tc in the r language , integrating several other topological measures in a user - friendly graphical framework . 
 the comparison of the tc statistic with other measures provides multifaceted insights on the dynamic processes shaping the evolution of pathogenic viruses . 
 the features and applicability of phylotempo were tested on serially - sampled intra - host human and simian immunodeficiency virus population data sets . 
 phylotempo is distributed under the gnu general public license at https://sourceforge.net / projects / phylotempo/. 
Introduction
Methods
Implementation
Results
Discussion
Supplementary Material
On the properties of the staircase-ness measure