Phylogenetics IV: Zika virus outbreak investigation

Phylodynamics

This practical is a continuation of the previous page. We are interested in the following questions:

• At what speed was Zika virus spreading through South America?
• Could we have predicted the total number of cases?

To address these questions we're going to transition from a principally phylogenetics outlook into a more population genetics focused approach. We're also going to shift from classical statistical methods such as maximum likelihood into a Bayesian approach.

What are the differences?

• Phylogenetics is principally concerned with the reconstruction of phylogenetic trees to understand how individuals are related.
• Population genetics views the phylogeny as secondary in importance, a means to the broader end of understanding population dynamics.

In the present context, the broader popgen aim is to quantify the growth rate of the Zika virus population during the outbreak and leveraging that understanding to predict the course of the epidemic.

• Classical statistics only uses the observed data to make inference about unknown parameters.
• Bayesian statistics combines the observed data with prior information to help parameter inference.

When you have reliable or well-informed prior beliefs then clearly using that information provides Bayesians with an advantage. There are practical and philosophical differences too that often make Bayesian approaches more tractable, particularly for complex problems involving many or high-dimensional parameters. However, when one does not have reliable or well-informed prior beliefs, the advantages of Bayesian statistics come at the cost of needing to find a sensible 'objective', 'reference' or 'non-informative' prior.

What is phylodynamics? One view is it's nothing more than a buzzword (a common affliction) meaning simply the application of population genetics to situations in which the underlying relatedness structure can be represented by a phylogeny (remember often this is not true, think of pedigrees).

Coalescent theory

Underlying modern population genetics is coalescent theory which provides a starting point for understanding how population processes leave signals in molecular sequence data by influencing the shape and size of the phylogenies or pedigrees relating the individuals in a population. We will only consider phylogenies here.

To understand how population processes perturb phylogenies, first it is necessary to know what to expect in the absence of interesting dynamics. The standard neutral model makes the following key assumptions:

• The population is well-mixed and stable in size
• There are no fitness differences between individuals
• Individuals can be randomly sampled from the population

You can think of these assumptions as describing an idealized population. In its simplest form, the coalescent provides a statistical probability distribution for phylogenies of individuals sampled at random from this idealized population.

In R, we are going to simulate phylogenies under the standard neutral model using the coalescent. The following code will simulate a new phylogeny (also called a coalescent tree) for 10 individuals every time you press return until you stop it by pressing Ctrl C

source("~/R/myutils/coalescent.R")
while(TRUE) {
  tree = sim.coalescent.tree(n=10)
  draw.tree(tree)
  readline()
}

☛ You can open a new R plotting window in Linux with the x11() command.

Try plotting 12 coalescent trees side-by-side on the same scale. What patterns do you notice in the trees?

par(mfrow=c(3,4),mar=c(0,2,1,0)+1)
for(i in 1:12) draw.tree(sim.coalescent.tree(10),ylim=c(0,5),lwd=2)

The y-axis on these simulated phylogenies shows time in coalescent time units. This is not calendar time and is different again to what you might call phylogenetic time units, in which phylogenies are usually scaled in terms of the expected number of substitutions per site. Instead, it is the number of generations divided by PN_e, where P is the ploidy (Zika virus is haploid so P=1) and N_e is the effective population size. For pathogen populations, the question of what does the effective population size correspond to biologically is an involved one. See e.g. the methods of Dearlove and Wilson (Philosophical Transactions of the Royal Society B 368: 20120314) if you are interested.

If you need to know, the conversion factor is

1 phylogenetic time unit = θ/2 coalescent time units

where θ/2 is the expected number of substitutions per site per coalescent time unit and θ can be estimated directly as the mean pairwise genetic distance.

The point of coalescent time units, confusing though they are, is that populations all behave the same under the standard neutral model once you've applied a linear scaling for the size of the gene pool.

For comparison, we're going to simulate trees in an exponentially growing population. The effect of exponential growth on phylogenies is one of the simplest signatures of population dynamics that is captured by coalescent theory.

par(mfrow=c(3,4),mar=c(0,2,1,0)+1)

for(i in 1:12) draw.tree(expgrowth(sim.coalescent.tree(10),5),ylim=c(0,1),lwd=2)

Notice that the scale of the y-axis is much smaller in these simulations. Since the population size is changing, it is necessary to specify that the time units are scaled by the size of the population at the present (i.e. at the time of sampling). The scale is smaller because the depth of a coalescent tree depends on the size of the population, so a population that was much smaller in the recent past can only have a shallow tree. What else do you notice comparing these simulations to the previous simulations under no growth?

The key points you should have noticed are:

• Coalescent trees are symmetric, meaning every topology is equally likely
• In the standard neutral model,
• There is a lot of variability in branch lengths
• On average, branches nearer the root tend to be the longest and tips tend to be the shortest
• In exponentially growing populations,
• The whole tree is compressed
• But the deeper branches are distorted more than the tips because further back in the phylogeny, the population was smaller
• This distortion of the relative length of branches near the root versus the tips compared to that expected under no growth is the signature of a growing population
• There is overlap in the distributions of coalescent trees with different growth rates, which means that estimates of growth rates are necessarily subject to substantial statistical uncertainty.

It is only because population growth leaves a signature in phylogenies, and through the phylogenies, on patterns of sequence variation, that we are able to estimate growth rates, and population dynamics more generally, from genetic data, subject to the model assumptions.

Often the phylogenies of an exponentially growing population are referred to as star-shaped. This is shorthand and shouldn't be taken too literally. Only under exceptionally strong growth would the distortion of branch lengths be so strong as to generate phylogenies appearing to have no internal branching structure, in which case the overall compression of the phylogeny would be extreme, likely resulting in very low diversity. Don't forget that estimation of the phylogeny itself is also subject to statistical uncertainty.

BEAUti and the BEAST

Only because of the exceptional service rendered to the scientific community by the Bayesian phylogenetics package BEAST are the irredeemably feeble acronyms used to name these programs tolerated.

BEAUti (Bayesian Evolutionary Analysis Utility) is a graphical interface that allows a very rich family of evolutionary models to be fitted to molecular sequence data in a Bayesian approach. BEAST (Bayesian Evolutionary Analysis Sampling Trees) implements the statistical analysis set up in BEAUti using a technique known as Markov chain Monte Carlo (MCMC). The results of BEAST can be interpreted using the packaged software Tracer or through standard software such as R.

We are going to set up an analysis of the Asian and South American Zika virus genomes to estimate the rate of global expansion of the outbreak from the earliest Malaysian sequence sampled in 1996 through to the 2015/2016 South American genomes. We're going to go back to the original sequence data to do this, so we need to extract just those sequences in R:

subtree.seq = match(sapply(subtree$tip.label,function(s) unlist(strsplit(s,"_"))[1]),rownames(a))
b = a[subtree.seq,]
rownames(b) = subtree$tip.label
write.fasta(b,"genomes.subtree.fasta")
Next, start BEAUti by opening the Terminal Emulator and typing beauti &

From the File menu, click Import Data and choose the genomes.subtree.fasta file just created.

You will see that beauti is organized by panels named Partitions, Taxa, and so on. We're going to go through these panels to set up the model and analysis.

Taxa. Use this panel to define groups of individuals for whom you would like to estimate the time of their most recent common ancestor. Previously we used TempEst to get an idea of the timing of the origin of the South American outbreak. To do the same analysis in BEAST, click the + in the bottom left corner, name the new taxon crown and use the panels to include all South American genomes, i.e. exclude only HQ234499.1, EU545988.1, JN860885.1, KF993678.1 and KJ776791.2

Repeat, but call the new taxon stem and select the tick box under the Stem? column. The stem and crown will date the beginning and end of the branch separating the South American sequences from the Asian sequences.

Tips. Similar to TempEst, use this panel to specify the sampling years. Tick the Use tip dates box and click the Parse dates button, selecting Defined by a prefix and its order, with order second and prefix underscore (_). The sampling years should appear in the Date column.

Sites. Use this panel to specify the substitution model. As previously I would recommend the HKY85 model with Empirical base frequencies and no rate heterogeneity, but you are free to experiment.

Clocks. I recommend the strict clock model for this analysis, which assumes a constant substitution rate throughout the phylogeny. Relaxed clocks may be useful but there is a trade-off, introducing more statistical uncertainty into estimates of dates.

Trees. Here you can specify the Exponential Growth tree prior. Advanced users may wish to experiment with other growth models including the nonparametric Bayesian skyline. I recommend you start with a distance-based (UPGMA) tree rather than a random tree because this will reduce the MCMC burn-in.

Priors. At this point it is necessary to state your subjective prior beliefs regarding the parameters you intend to estimate. If you are not well-informed, you can either

• Conduct prior sensitivity analysis, in which you run the whole analysis with different priors and test their effects on the final results.
• Aim to specify some sort of objective/reference/non-informative prior.

For tmrca(crown), tmrca(stem) and treeModel.rootHeight do not change the default Using Tree Prior. This option leads to non-sensical inference because it puts multiple priors on the same parameters.

For kappa, the default prior is log-normal with a peak around 3. This looks like a reasonable informative prior since transition:transversion ratios are known to generally exceed 2.

For exponential.popSize, the default prior is a non-informative one over x prior. This is equivalent to placing equal prior weight on any value of the logged parameter, and is a reasonable-looking 'reference' prior. An alternative would be Infinite Uniform, equivalent to placing equal prior weight on any positive parameter value. Both priors suffer the major disadvantage of failing to guarantee that the MCMC will converge on a proper posterior distribution, a problem that may come to light if the data are not sufficiently informative as to the parameters.

For exponential.growthRate, the default prior is a double exponential (Laplace) distribution with mean absolute value (the scale) of 1.0. This implies equal prior weight to population contraction as to growth. An alternative would be an exponential distribution that only allows positive values, corresponding to population growth. However, this would make it impossible to know whether the data are able to rule out population contraction without help from the prior. The choice of scale is subjective, and an alternative 'reference' prior would be the Infinite Uniform. If unsure (you should be), try conducting a prior sensitivity analysis.

MCMC. Specify the File name stem as beast.mcmc1 and click Generate BEAST File, saving it as beast.mcmc1.xml in your Downloads directory. Now repeat with File name stem beast.mcmc2 and this time save as beast.mcmc2.xml. This is essential because multiple runs are required to establish whether the MCMC has successfully converged to a proper posterior distribution.

Running BEAST. The next step is to run the two chains. Open two Terminal windows and navigate in both to your Downloads folder. In one type

beast beast.mcmc1.xml

and in the other type

beast beast.mcmc2.xml

This will run both cores flat out at 100%, which may make other programs run slowly. You can lower the priority of the BEAST runs by first identifying their process numbers using top or ps (the processes are called java) then typing

renice +5 PROCESSNUMBER

It took a little over 10 minutes to run the two chains in parallel on my virtual machine.

Understanding the output of BEAST is the next page of this practical.