Phylogenetics II: Zika virus outbreak investigation

Phylogenetic reconstruction

This practical is a continuation of the previous page

Distance-based methods

The simplest methods for phylogenetic reconstruction base inference directly on the pairwise genetic distance of the sequences. This throws away most of the information in the alignment by summarizing it through a simple summary statistic. For these simple distance-based methods, it is sufficient to analyse the data in R.

#Inspect the pairwise distance matrix 
PD = matrix(0,nrow(a),nrow(a)) 
PD[lower.tri(PD)] = DIST 
PD[upper.tri(PD)] = t(PD)[upper.tri(PD)] 
rownames(PD) = rownames(a); colnames(PD) = rownames(PD) 
#Only a limited number of rows and columns will fit on the screen 
PD[1:5,1:5] 

#Visualize as an image: because we used the --reorder option in mafft, this immediately reveals relatedness structure in the data 
image(1:nrow(a),1:nrow(a),PD,xlab="Individual",ylab="Individual",main="Pairwise distance") 

#One of the simplest ways to infer a phylogenetic tree is with UPGMA (Unweighted Pair Group Method with Arithmetic Mean). This is a clustering algorithm that iteratively joins the branches ancestral to the most similar sequences, computing similarity iteratively as the mean between each pair of groups that are merged when their branches are joined. 
upgma.tree = as.dendrogram(hclust(as.dist(PD),method="average")) 
plot(upgma.tree,xlab="",ylab="Distance") 

#The phangorn package allows us to use neighbour joining, another distance-based method, instead.
require(phangorn)
nj.tree = NJ(PD)
plot(nj.tree)

#The resulting tree needs to be rooted. A simple method is midpoint rooting.
nj.tree = midpoint(NJ(PD))
plot(nj.tree)

☛ You can save your current R session to the working directory with the save.image() command.

Distance-based methods are quick and often referred to as dirty. Among the limitations of distance-based methods

• The methods are ad hoc algorithms that do something sensible at each step, but do not add up to a principled approach to inference
• In particular, there is no explicit model of evolution, making the assumptions difficult to scrutinize
• With no probabilistic model of sequence evolution, no theoretical guarantees of performance can be applied to these methods
• For the same reason, it is difficult to quantify statistical uncertainty associated, although bootstrapping is commonly employed
• Since information is thrown away when summarizing the sequences by pairwise distance alone, distance methods cannot be expected to perform as well as alternatives based on the full data
• In practice, these concerns are born out with worse empirical performance than more sophisticated alternatives

Distance-based methods are certainly quick, but in practice their accuracy is often not horribly worse than sophisticated alternatives. While they are useful data visualization tools, there is no good reason not to apply more accurate methods wherever possible.

Maximum likelihood methods

Maximum likelihood (ML) methods apply an explicit, probabilistic, model of sequence evolution to the full sequence data, and estimate the tree and branch lengths through a principled statistical approach. ML methods are

• Consistent, meaning that as more data becomes available, the estimate gets better
• Unbiased for large datasets
• Amenable to formal hypothesis tests and construction of confidence intervals

ML methods attempt to find the tree and branch lengths that are most compatible with the data. They consider how probable the observed data would have been for different trees and branch lengths. No matter how improbable the observed data are, they optimize the tree and branch lengths so that no alternative tree and branch lengths could make the observed data more probable.

In practice, it is usually impossible to guarantee that a computer algorithm has found the global maximum likelihood estimates because the number of tree topologies is too large to explore exhaustively: for n sequences there are

(2n-5)(2n-7)(2n-9) ... (5)(3)(1)

possible unrooted bifurcating labelled trees. Instead, various search strategies are employed that aim to find the best possible estimate in a reasonable amount of time. These algorithms often use distance-based methods as a starting point.

Since maximum likelihood methods must be highly optimized, most widely-used methods are available as stand-alone command line programs. We will use PhyML which is simple to use and fast up around 200 sequences. Beyond that there are more specialized tools such as RAxML.

Open an existing docker session or execute a new one if necessary. The following command will launch PhyML

cd ~
phyml -i genomes.relab.phylip -b 0 -v 0 -c 1 -s BEST

Typing phyml -h will list all the options. The above command line implies the following defaults:

• -d nt specifies nucleotide not amino acid sequences
• -p instructs PhyML to initialize the search with a minimum parsimony (distance-based) tree
• -m HKY85 specifies the Hasegawa, Kishino and Yano (1985) model of nucleotide sequence evolution
• -f e indicates that nucleotide frequencies should be estimated from the empirical frequencies in the alignment
• -t e implies that the transition:transversion ratio parameter should be estimated by maximum likelihood along with the tree and branch lengths

The explicitly selected options do the following:

• -i specifies the Phylip format sequence alignment
• -b 0 instructs PhyML not to perform bootstrap replicates
• -v 0 -c 1 indicates that a single substitution rate should be applied to all columns in the alignment
• -s BEST combines NNI and SPR search strategies

When I ran through the practical, this command took 40 seconds to run. Note that bootstrapping is considered good practice as a means of quantifying the uncertainty in the estimated tree. Typically, 100 or more bootstrap replicates are conducted and the tree plotted with bootstrap support (0-100) beside each split. PhyML can perform regular and approximate bootstrap, and you are free to explore these capabilities.

The tree can be visualized, for example in figtree, by typing

figtree genomes.relab.phylip_phyml_tree.txt &

☛ You can open a new window or tab in Terminal from the File menu

Figtree has a flexible and intuitive interactive graphical interface that can be used to fine-tune how the tree is displayed. Alternatively, switching back to R, the tree can be read in and displayed

require(phangorn)
ml.tree = ladderize(midpoint(read.tree("genomes.relab.phylip_phyml_tree.txt")))
plot(ml.tree)
#Output the midpoint rooted version for later use
write.tree(ml.tree,"genomes.midpoint.tree.txt")

For users with the time to learn, R provides even more fine-grained control over the way the tree is displayed. See the ape package documentation for details.

☛ In R access help by preceding any command with a question mark, e.g. ?ladderize

Among the benefits of ML methods is the interpretability of the estimated tree. For instance, the branch lengths report the expected number of substitutions per site along the branch, given the data and the assumed model of sequence evolution. This provides an alternative to the 'raw' pairwise genetic distance by calculating the sum of branch lengths between each pair of sequences, a quantity known as the phylogenetic or cophenetic distance. In R, it can be obtained by

# Pairwise phylogenetic distance matrix
TPD = cophenetic(ml.tree)
summary(TPD[lower.tri(TPD)])
# Mean branch length 
mean(ml.tree$edge.length)

When I ran the code, I obtained a mean phylogenetic distance of 0.14, indicating that on average 14 substitutions per 100 bases separate each pair of genomes, somewhat higher than the mean pairwise genetic distance of 0.087. Usually phylogenetic distance is higher than genetic distance because

• ML methods can account for the saturating effect of repeat mutation on 'raw' genetic distance
• Missing bases in the alignment caused by ambiguity or indels can be imputed by ML methods, and may increase the distance
• Recombination can artificially inflate phylogenetic distances by giving the appearance of repeat mutation

The problem of recombination and detecting it in phylogenies is the subject of the next section.

Rooting out recombination

A key assumption of phylogenetics methods is that recombination, also known as horizontal gene transfer, has not played an important role in the ancestral history of the sample. In other words, a single tree is assumed to accurately represent the relatedness of the individuals at all positions in the sequence.

This contrasts strongly with the assumption of a pedigree - e.g. your family tree - in which recombination is expected to have occurred every generation. This is why you get half of each paternally-inherited autosome (chromosomes 1-22) from your paternal grandfather and grandmother respectively, and half of each maternally-inherited autosome from your maternal grandfather and grandmother respectively. So there is a change in ancestry along every chromosome for every generation in your pedigree. This means that relatedness in any population changes constantly throughout the human autosome sequences.

Trying to reconstruct a single tree from highly recombining sequences is a case of model misspecification and risks misleading inference. This is why it is important to

• Ask whether it is sensible to reconstruct a phylogeny in the first place
• Test for evidence of recombination in the estimated phylogeny

Unlike humans, most infectious organisms are not obligately sexual, and therefore recombination might be rare or absent. In this practical, we will use ClonalFrameML to test for evidence of recombination in the sequences.

	A	C	G	T
A	-	fC	κfG	fT
C	fA	-	fG	κfT
G	κfA	fC	-	fT
T	fA	κfC	fG	-

The HKY85 substitution model
The rate of substitution from one base to another is assumed to be proportional to the overall frequency of that base (f_A,f_C,f_G,f_T), adjusted by a factor κ representing the relative rate of transitions (purine-to-purine or pyrimidine-to-pyrimidine) versus transversions (purine-to-pyrimidine or vice versa).

ClonalFrameML assumes an HKY85 substitution model and takes the ML phylogeny as input. To extract the ML estimate of the transition:transversion ratio from the PhyML output, type at the command line

grep 'Transition/transversion ratio' genomes.relab.phylip_phyml_stats.txt

Inserting this value into the -kappa argument (I got 4.6), run ClonalFrameML by typing

ClonalFrameML genomes.midpoint.tree.txt genomes.relab.fasta zika -kappa 4.6 -prior_mean "0.1 0.001 0.1 0.015" -prior_sd "0.1 0.001 0.1 0.015" > zika.out.txt
Rscript /tmp/ClonalFrameML/src/cfml_results.R zika
firefox ~/zika.cfml.pdf &

There are a large number of options which are summarized by running ClonalFrameML with no trailing arguments.  We have modified only two arguments compared to their default settings, the transition:transversion ratio as mentioned above, and the mean and standard deviation (sd) of the prior distribution of the branch lengths of the tree. ClonalFrameML is a Bayesian method (more on this later) and uses prior knowledge to help inform its parameter estimates. The prior knowledge on four parameters is represented by gamma distributions with means and standard deviations specified on the command line. The fourth parameter is the branch length, and I have set this to have a prior mean and standard deviation roughly equal to the mean branch length in the ML tree.

This analysis took 45 seconds to run for me. ClonalFrameML produces a pdf file showing the recombination-corrected phylogeny against a plot of the positions of substitutions (small vertical bars) and recombination events (long horizontal dark blue lines) in the ancestry history of the sequences.

The presence of multiple long horizontal dark blue lines indicates there have been numerous recombination events in the history of the sequenced genomes. You can see that in any particular row, these recombination events pick out heterogeneity in the substitution density across the genome (left to right) on the corresponding branch of the tree.

The deep branches of the phylogeny in particular seem to have been the focus of recombination. This suggests recombination may have been an important force early in the evolutionary history of these genomes. In contrast, there is very little recombination near the tips of the tree, indicating that recombination has been rare or absent during recent evolution.

Most of the recombination appears to have occurred on the branches leading to the Senegal sequences KF383120 (ArD 142623), KF383117 (ArD 128000) and KF383118 (ArD 157995), which were sampled between 1997 and 2001, and the branch separating KF383120 and DQ859064 (the Spondweni virus outgroup). This will be important to bear in mind when interpreting the phylogeny, particularly because recombination can distort the tree and cause misleading inference of relationships between sequences. Where there has been no recombination between sequences, we can reconstruct the phylogenetic subtree with more confidence.

This practical continues on the next page.

Phylogenetics I: Zika virus outbreak investigation

This tutorial now runs in Docker. Instructions for setting up under MacOS are here and setting up under Ubuntu are here.

Before beginning the tutorial, create a working directory to store your results. Open a Terminal and type

WDIR=~/Downloads/phylogenetics

mkdir $WDIR

To start the tutorial, type in a Terminal (slightly different for Ubuntu: see here)

docker run -v $WDIR:/home -e DISPLAY=$(ifconfig en0 | grep 'inet ' | cut -f 2 -d' '):0 -it --name practical dannywilson/teaching-phylogenetics:march2020

During the tutorial, you will need to open multiple sessions. To start a new session open a new Terminal window and type

docker exec -it practical /bin/bash

To check the graphics are working correctly, type in your docker session

xclock &

which should open a new window showing a clock face with hour and minute hands.

Introduction

At the beginning of 2015, several cases of patients presenting symptoms of mild fever, rash, conjunctivitis and joint pain were reported by clinicians in South America. Investigation revealed infections with Zika virus, an obscure relative of dengue and chikungunya, that until the late 2000s had caused sporadic benign infections in Africa and Asia. Outbreaks emerged in Micronesia in 2007 and French Polynesia in 2013-14, during which an unusual increase in the incidence of severe neurological complications was first noticed. Over the course of 2015 and into 2016, clinicians in Brazil reported an upswing in the numbers of children born with microcephaly. The finger of blame was pointed towards Zika virus, a flavivirus, members of which group had been known to exhibit neurotropism and thereby carry the risk of causing neurological defects. As of March 2016, the epidemiological history and expectations for the ongoing spread of Zika virus were very uncertain.

The purpose of this practical is to demonstrate the tools available for phylogenetic analysis during early outbreak investigation. The aims are to address the following questions:

• What was the geographic origin of the 2015/2016 Zika outbreak?
• How was the virus related to previous Zika outbreaks?
• When did the 2015/2016 Zika virus strain diverge from other variants?
• Was there anything special about the 2015/2016 strain?
• At what speed was Zika virus spreading through South America?
• Could we have predicted the total number of cases?

Genome sequences

As of January 2016, only a modest number of Zika virus whole, or near-whole, genome sequences were available:

Click here for the Google Doc. The first step in the analysis is to obtain the whole genome sequences from Genbank using the accession numbers provided in publications or found through a Genbank search query.

In your web browser, navigate to https://www.ncbi.nlm.nih.gov/genbank/

☛ You can open Firefox in your docker session by typing firefox &

Copy the accession numbers from the spreadsheet and paste them into the query box, and click search. The 31 sequences will appear, split across two pages.

To show the sequences, click Summary, FASTA, and to download click Send to: File. (Leave Complete Record checked). By default, it will save the genomes to a file called sequence.fasta.txt in your working directory.

☛ Your working directory has different names from inside and outside the docker container. From inside docker, its location is /home. From outside docker it is wherever you set WDIR to reference.

Multiple sequence alignment

Phylogenetic reconstruction relies on the identification of character traits whose patterns of presence and absence is informative about shared evolutionary ancestry. In the past, many traits were used, but modern phylogenetics utilizes, almost exclusively, molecular sequence traits: DNA, RNA or protein sequences.

To compare any character across individuals, first homology must be established. Only traits that are homologous, i.e. directly co-inherited, can be used as the basis for valid comparison. Therefore the first challenge is to identify sequence homology between the virus genomes. By visualizing the Zika genomes in jalview, you will see they are all different lengths, partly because of differences in primers used for sequencing, partly because of mutations.

Open jalview by opening a docker session in Terminal and typing jalview &. (Ignore the annoying messages it opens.)

In jalview, open the genomes you downloaded via the File, Input Alignment, From File menu, and navigating to sequence.fasta.txt in your working directory, which is /home from inside docker.

You can assist visualization by selecting, from the menu, Colour, Nucleotide. Scanning through the sequences you will notice that there are similarities between some of the genomes, but it is very messy and inconsistent.

We're going to use mafft to align the sequences by identifying regions of sequence similarity. Sequences will be padded with indels (- characters, representing insertions or deletions) to improve the similarity within columns of the sequence alignment.

Open a docker session in the Terminal. Typing mafft -h will produce a bewildering list of all the mafft options. Multiple sequence alignment is a very difficult problem in bioinformatics and even the most popular aligners can only be considered heuristic at best. In this practical, I will give you a recipe for aligning the sequences. If you want to learn more, Google for mafft and read the online documentation as a starting point.

In docker, navigate to your working folder and run mafft as follows (increase the --thread option if you want to use more than two processors)
cd ~
mafft --thread 2 --anysymbol --bl 62 --op 1.53 --ep 0.123 --reorder --retree 1 --treeout --maxiterate 0 --localpair sequence.fasta.txt > genomes.mafft.fasta

☛ Some users have found they need to replace sequence.fasta.txt with sequence.fasta in the command above. Type cat sequence.fasta.txt  to check the file exists and is not empty.

This will take 10-20 minutes with two processors depending on your machine. Once the alignment is complete, you can open the genomes.mafft.fasta file in jalview. If the process says it was killed, you may need to increase the maximum system memory docker is allowed to use.

A good alignment between sequences that are genuinely homologous should

• Maximize the percent identity in each column of the alignment
• Produce a sequence alignment which is not spuriously longer than the original sequences
• Not contain indels for the majority of sequences in most parts of the alignment

Between closely related individuals, the sequences should be identical for most of their length. The quality of the alignment is critical to all downstream analyses because they are based on the differences that remain between the sequences after alignment. These differences are phylogenetically informative, and they can be better visualized by selecting from the menu Colour, Percent Identity.

Basic alignment summaries

Before moving on to phylogenetic reconstruction, we're going to quantify the alignment quality in R. In a new docker session, type R and follow the instructions below.

#To load the functions for reading the alignment, type
source("https://raw.githubusercontent.com/danny-wilson/danny-wilson.github.io/main/files/fasta.R")
#Read in the alignment
setwd("~")
a = toupper(read.fasta("genomes.mafft.fasta",as.char=TRUE))

#Number of sequences 
nrow(a)

#Alignment length
ncol(a)

#Maximum sequence length
max(rowSums(a!='-'))

#Number of sites with one or more indel: for closely related sequences, most sites would not be expected to be indels, although differences in primers may exacerbate indel rates near the ends of the sequences
sum(colSums(a=='-')>0)

#Total frequency of indels in the alignment
mean(a=='-')

When I ran through the practical, the total alignment length was 10835 bases, only 0.2% longer than the longest sequence (10808 bases), indicating that mafft had found the sequences to be largely collinear, and had introduced indels only sparingly. Although 997 columns had at least one indel, only 3.1% of all bases were indels. These summary statistics indicated no major problems with the alignment process.

Tidying up the alignments

Nevertheless, sequence alignments can contain a variety of idiosyncrasies that are helpful to address prior to downstream analysis so they do not cause analysis artifacts. It is also helpful to tidy up the sequences before quantifying the sequence similarity of aligned bases.

#Nucleotide composition: should be in keeping with known base frequencies in the organism.
table(a)

#Because of sequencing quality limitations, there are often ambiguous base calls such as R (purine, A/G), Y (pyrimidine, C/T) and N (any), in addition to the indel symbol (-). Downstream analyses differ in their treatment of such ambiguity codes. However, they are all roughly equivalent in the sense that they convey little or no information as to what the base would have been. Therefore to simplify matters, it is reasonable to replace them all with 'N'. 
gd = a=='A' | a=='C' | a=='G' | a=='T' | a=='U' | a=='N'
a[!gd] = 'N'

#For RNA sequences, there can be inconsistency in whether uracil or its DNA counterpart thymine is reported. When sequences come from different sources, they might use different conventions. To simplify matters, replace all uracil codes with thymines.
a[a=='U'] = 'T'

#In downstream analyses, computational complexity strongly depends on the number of different phylogenetic patterns, i.e. the occurrence of As/Cs/Gs/Ts/Ns across individuals in a single column of the alignment. Ambiguous base calls that are scattered throughout the alignment fairly randomly can massively increase the number of different phylogenetic patterns. Particularly for columns that are otherwise invariant, this added cost comes with no benefit because only variable bases are phylogenetically informative. Therefore, it is reasonable at invariant sites to replace Ns with the consensus base.
cons = apply(a,2,function(x) names(sort(table(ifelse(x=='N',NA,x)),decreasing=TRUE))[1])
is.variable = apply(a,2,function(x) nlevels(factor(ifelse(x=='N',NA,x)))>1)
a[,!is.variable] = sapply(which(!is.variable),function(j) ifelse(a[,j]=='N',cons[j],a[,j]))

#The sequence labels in Genbank are often long and contain unusual characters that can upset downstream analysis programs. Truncate them at the first space
lab = sapply(rownames(a),function(s) unlist(strsplit(s,' '))[1])

#Write the modified sequences in FASTA format
rownames(a) = lab
write.fasta(a,"genomes.relab.fasta")

#One of the bugbears of bioinformatics is file format conversion. Many programs use different file formats. In the phylogenetics setting, PHYLIP format is common. For later use, we're also going to write the sequences in PHYLIP format.
outfile = "genomes.relab.phylip"
nc = max(nchar(lab))
lab = sapply(lab,function(s) paste0(s,paste0(rep(" ",nc-nchar(s)),collapse="")))
cat(nrow(a),ncol(a),file=outfile)
cat("\n",file=outfile,append=TRUE)
for(i in 1:nrow(a)) {
  cat(lab[i],paste0(a[i,],collapse=""),file=outfile,append=TRUE)
  cat("\n",file=outfile,append=TRUE)
}

Basic diversity measures

As a way of familiarizing yourself with the sample, it is helpful to quantify the relatedness of the individuals in terms of nucleotide similarity and total number of variable sites. Basic measures of pairwise genetic diversity between individuals also provide the starting point for distance-based phylogenetic tree building.

#Number of variable sites: for closely related sequences, most sites would not be expected to vary
(S = sum(colSums(t(t(a)!=cons))>0))

#Mean pairwise diversity: the mean, across all pairs of sequences, of the proportion of the sequence at which they vary. For closely related sequences, should be much less than 0.1
a2i = c('A'=1,'C'=2,'G'=3,'T'=4)
i = matrix(a2i[a],nrow(a),ncol(a))
DIST = rowMeans(apply(i,2,dist)>0,na.rm=TRUE)
(PI = mean(DIST))

#Distribution of pairwise diversity across sequences
summary(DIST)

The number of variable sites (4362/10835) and the mean pairwise diversity (0.087) are perhaps on the high side. Since the alignment quality metrics gave no cause for concern, this could indicate the presence of one or more outlier genomes in the sample, or it might simply reflect the high inherent diversity of the viruses.

This practical continues on the next page

Create a simple MLST pipeline

The idea is to obtain the multi locus sequence type (MLST) for an assembled genome by building a simple bioinformatics pipeline. The starting material needed is

1. One or more assembled genomes, in FASTA format, to multi locus sequence type
2. A multilocus sequence typing database comprising
1. For every locus, the allele sequences in FASTA format labelled by their allele numbers
2. A list of sequence types (STs) defined by the multi-locus allele numbers (allelic profile)
3. A linux server with blastn installed

I will illustrate by making a simple pipeline for Staphylococcus aureus.

Phylogenetics in Practice Lecture Downloads

Powerpoint slides of the lectures can be downloaded from here:

phylogenetics_I.pptx
phylogenetics_II.pptx

Instructions for launching the Azure Virtual Machine from Ubuntu Linux

1. Once you have accepted your invitation to the Oxford DTC group, navigate to portal.azure.com
2. In the search box at the top of the dashboard type Phylogenetics-WILSON-image
3. Click the button at the top called +Create VM
4. Under Basics, make the following settings
   Name: Phylogenetics-YOURSURNAME
   User name: Choose a username (you won't use it)
   Password: Choose a password (you won't use it)
   Resource group: Use existing, 3_Phylogenetics_Wilson
5. Under Size, choose D2S_V3 Standard (£64.32/month)
6. Under Settings, use the defaults except
   Network security group (firewall): select the existing Phylogenetics-nsg
   Auto-shutdown: On
7. Click the Bell icon at the top of the screen. It will open a panel that says Deployment in progress. Once the deployment has succeeded, click Go to resource
8. Click Connect and copy the IP address after the @ sign
9. Open Remmina on your local computer
10. From the menu selection Connection, New and enter the following settings
    Connection name: Phylogenetics
    Server: paste the IP address you copied
    User name: wilson
    Password: wilson17
    Colour depth: True colour (24 bpp)
    Then click the Save button
11. Double-click the new Phylogenetics connection you set up
12. If the connection fails first time, try username wilson password wilson17 again in the box that appears.

Remember to stop your virtual machine after you are done, otherwise the DTC will continue to be charged. Also make sure any aborted attempts to create VMs are stopped.

Instructions for launching the Azure Virtual Machine from Windows or Mac

1. Once you have accepted your invitation to the Oxford DTC group, navigate to portal.azure.com
2. In the search box at the top of the dashboard type Phylogenetics-WILSON-image
3. Click the button at the top called +Create VM
4. Under Basics, make the following settings
   Name: Phylogenetics-YOURSURNAME
   User name: Choose a username (you won't use it)
   Password: Choose a password (you won't use it)
   Resource group: Use existing, 3_Phylogenetics_Wilson
5. Under Size, choose D2S_V3 Standard (£64.32/month)
6. Under Settings, use the defaults except
   Network security group (firewall): select the existing Phylogenetics-nsg
   Auto-shutdown: On
7. Click the Bell icon at the top of the screen. It will open a panel that says Deployment in progress. Once the deployment has succeeded, click Go to resource
8. Click Connect and copy the IP address after the @ sign
9. Open Microsoft Remote Desktop on your local computer
10. Click New and enter the following settings
    Connection name: Phylogenetics
    PC name: paste the IP address you copied
    User name: wilson
    Password: wilson17
    Click the red close button in the top left of the Edit Remote Desktops window
11. Double-click the new Phylogenetics connection you set up
12. Select Connect Always when prompted

Remember to stop your virtual machine after you are done, otherwise the DTC will continue to be charged. Also make sure any aborted attempts to create VMs are stopped.

Update. If Remmina fails to connect open a Terminal and type

1. ssh -X wilson@IPaddress
2. Enter password wilson17

You can run programs remotely by typing, e.g. firefox &. You can open as many terminal windows or tabs as you like to help organize your remote session.

Phylogenetics V: Zika virus outbreak investigation

This is a continuation of the previous page of the practical.

The aim of this part of the practical is to understand the output of the Bayesian phylogenetics/population genetics software BEAST. The output can be read in by standard software such as R, but BEAST has its own bespoke graphical interface called Tracer that we will explore.

Tracer. To visualize the output of the BEAST analysis, we will use Tracer. To launch it, at the Terminal type tracer &

From the File menu, select Import Trace File, navigate to your Downloads folder and select beast.mcmc1.log. Repeat to read in beast.mcmc2.log. Notice that the two files are listed in the top left of the screen, and below that is a Combined output of the two runs.

Assessing MCMC mixing and convergence

Click on clock.rate in the list of parameters on the left side bar of tracer. Select the Trace panel

Sometimes the trace is likened to a hairy caterpillar. What it shows is the way that the parameter clock.rate was sampled (i.e. simulated) from the posterior distribution by the MCMC over the iterations (labelled State on the x-axis).

You can see that the MCMC explores the posterior distribution by taking a random walk through the parameter values. This causes the evident auto-correlation in parameter values from iteration to iteration.

The first million iterations are greyed out. This phase of the run is called the burn-in before which the MCMC has reliably converged on the posterior distribution from its initial parameter values. Identifying the length of the burn-in is an inexact science, and we just used the BEAST default of one million iterations.

Following the burn-in, the hairy caterpillar should look fairly flat with relatively modest auto-correlation. If the caterpillar is obviously skewiff or the auto-correlation shows a long lag (evident as waves) this means that the MCMC is not mixing well, i.e. the random walk is not efficiently exploring the parameter values. A poorly-mixing MCMC needs running for longer, or the proposals used in the random walk need optimizing.

A useful heuristic for how well the MCMC is mixing is the ESS (effective sample size), shown in the left side bar for each parameter. The ESS, which is calculated for all iterations after the burn-in, takes account of the auto-correlation to calculate the approximate number of independent samples from the posterior distribution. An ESS below 100 suggests serious problems with the MCMC. Every parameter listed in the left side bar needs an adequate ESS.

The next step is to evaluate convergence of the MCMC. This was the reason for running multiple chains. In the top left panel labelled Trace Files, click first on beast.mcmc1.log and then hold Ctrl and click beast.mcmc2.log. As you Ctrl-click the second run, you should see its trace super-imposed on top of the trace of the first run. If the two runs have both converged, they should lie on top of each other. You can improve the graph by dragging the right frame of the Tracer window to widen the window until the Colour by drop-down list is visible. Select Colour by: Trace File. This will make the second trace purple. The two traces should both be well-mixed and should not show systematic differences such as a different mean value. Failure to converge is often accompanied by bad mixing, which will make the two traces clearly non-overlapping during some periods of the MCMC.

If you are satisfied that for every parameter the burn-in is long enough, the ESS is large enough and the MCMC shows good mixing and convergence, it is time to move on.

Interpreting the Bayesian posterior distribution

The molecular clock rate is a fundamental quantity for the interpretation of phylogenies because it provides the real-time substitution rate. Only with this information is it possible to convert phylogenetics time units or coalescent time units to calendar time. Previously we estimated the clock rate using TempEst. Now we will use BEAST to estimate this parameter.

Click the Combined run from the Trace Files list in the top left corner, select clock.rate from the left side bar and choose the Marginal Prob Distribution panel on the right.

This shows you the posterior distribution of the clock rate approximated by the MCMC: it is essentially a histogram (technically, Tracer shows a smoothed kernel density estimate) of the parameter values sampled by the MCMC runs that you visualized in the trace plots.
The posterior distribution is a Bayesian statement about the probable values of the parameter. It takes into account the prior distribution you specified in BEAUti and the observed data as interpreted through the evolutionary model you specified in BEAUti.

How do you use the posterior distribution?

Point estimates. If you want to quote a single number to represent your estimate of the parameter, you can summarize the posterior by taking some form of average, such as the:

• Posterior mean. This represents the expected value of the parameter, averaging over the uncertainty in the posterior.
• Posterior median. There is 50% posterior probability that the true parameter lies below this value and 50% that it lies above it. Unlike the mean, it has the desirable property that the median of any transformation (e.g. logarithm) of the parameter equals the transformation of the median. This property is known as invariance to transformation and is shared with maximum likelihood estimates.
• Posterior mode. This is the value with the highest posterior density, so in a sense represents one's 'best guess'.

If you want to know more about how to motivate your choice of point estimate, look up Bayesian decision theory, which will introduce the concept of risk/utility functions.

To get Tracer to provide the point estimate of your choosing, click on the Estimates panel. I obtained a posterior mean of 7.6×10^-4. Do not blithely copy out 7.6E-4, because this is not scientific notation. Do not copy out 7.6413×10^-4 unless (a) you are confident that the MCMC provides you with 5 significant figures of precision and (b) this level of precision is needed. BEAST helpfully provides a standard error of the mean. I got 3.6×10^-6 which indicates that the true posterior mean, after accounting for the limited number of iterations of the MCMC, lies somewhere between the sample mean plus and minus two standard errors, i.e. 7.6413×10^-4 ± 2×3.6×10^-6 = (7.569×10^-4, 7.713×10^-4). This indicates that only one or two significant figures are appropriate. To obtain greater precision combine more chains or run each chain for longer. As a guide, the precision increases only with the square root of the number of iterations, so 100 times as many iterations are needed for 10 times the precision.

I obtained a posterior median of 7.6×10^-4 as well. For symmetrical distributions, the mean and median will match, as here. For some reason, Tracer could not provide the posterior mode but it would probably have been similar.

Credibility intervals. The point estimate does not convey the statistical uncertainty associated with your inference. (This is different to the uncertainty represented by the standard error of the mean.) Common statistics used to summarize the uncertainty in the posterior distribution include the:

• 95% highest posterior density (HPD) interval. This is the narrowest interval within which 95% of the posterior density lies.
• 95% equal-tailed interval. This is the interval defined by the 2.5% and 97.5% percentiles. There is 2.5% probability the true parameter value lies below the interval and 2.5% probability it lies above it. Like the median, this interval is invariant to transformations of the parameter.

Of course you can choose 90% instead of 95% credibility intervals, or other values, but you must state what they are. I obtained a 95% HPD interval of (5.3×10^-4, 9.9×10^-4). Tracer does not produce the equal-tailed interval, even though it is easier to calculate.

Unlike the semantic contortions of classical confidence intervals, a Bayesian can say "I believe there is a 95% probability that the clock rate is between 5.3×10^-4 and 9.9×10^-4."

This is a respectably narrow credibility interval, in the sense that there is only around a two-fold difference between the upper and lower bound.

• How do BEAST's estimates of the clock rate compare to TempEst's?

The age of the South American Zika virus outbreak

Now you're equipped to interpret the other evolutionary parameters. While you could manually take the clock rate and apply it to one of your previously estimated phylogenies to try to date the emergence of the South American outbreak, this is unnecessary because BEAST co-estimates the phylogeny along with the evolutionary parameters and, because we set it up in BEAUti, it reports the age of the most recent common ancestor (MRCA) of the South American sequences in tmrca(crown) and the age of the MRCA of the South American sequences and their next most closely related non-South American sequence in tmrca(stem). It is important to estimate both quantities because, even if the South American outbreak is monophyletic, these MRCAs can only provide bounds on (i.e. book-end) the date the outbreak began. Note that because we provided the sampling dates of the tips in years, the MRCA date estimates are also in years.

• Can you calculate the credibility interval within which there is at least a 95% posterior probability that the true date of the origin of the South American outbreak began?

How rapidly did the Asian/South American Zika virus pandemic spread?

The exponential.growthRate parameter is the rate of population increase (or decrease, if negative) per year. What does the posterior distribution look like? How quickly did the outbreak spread? Do you obtain a different rate if you analyse the South American sequences alone?

Advanced users might want to explore BEAUti for other methods of estimating population size changes in the Asian/South American Zika genomes or the South American genomes alone. What does the non-parametric Bayesian Skyline method infer?

Bayesian inference of the phylogeny

So far we have seen how BEAST estimates evolutionary parameters and tree statistics in a Bayesian fashion, but we have not seen how it estimates the phylogeny, also known as a genealogy in a population genetics context. A major difference between the way PhyML and BEAST estimate phylogenies is that the former estimates an unrooted tree in which the branch lengths are in phylogenetic time units (i.e. expected numbers of substitutions per site). In contrast, BEAST is able, when provided with a clock rate specified through the prior or estimated from informative data, to estimate a rooted tree in which the branch lengths are in calendar time. This rooted tree is constrained so that the tips must occur at the specified sampling times, whereas the unrooted tree estimated by PhyML is unconstrained.

Another major difference is that PhyML estimates a point estimate (the maximum likelihood tree) whereas BEAST samples trees from the posterior distribution. Therefore it is necessary to summarize this posterior distribution in some way to provide a point estimate. A consensus tree can be constructed which incorporates nodes that have high posterior probability, as evidenced by their appearance in the majority of sampled trees. The branch lengths might be specified, for example, by the posterior median of each branch taken as an average over the sampled trees in which it appears.

To build a maximum clade credibility (MCC) tree in BEAST:

• Launch TreeAnnotator (type treeannotator at the command line)
• Set the Burnin to 1000 samples
• Set the Input Tree File to beast.mcmc1.trees
• Select a name for the Output File, e.g. beast.consensus.tree
• Click Run

To view the consensus tree, launch figtree. As for the individual parameters, it is important to quantify the uncertainty associated with the phylogeny point estimate. Explore the options to add labels to the internal nodes of the tree, among which one can choose the posterior probability of the split. If you obtained bootstrap values previously for the maximum likelihood tree, how do the Bayesian posterior probabilities compare?

Hackathon ideas

• There are now over 600 whole or near-whole Zika virus genomes. Can you apply what you've learned about reference-based mapping earlier to extend your phylogenetic analysis to leverage the information containined in hundreds of genomes? 
• Can you apply what you've learned to understand the early dynamics of other outbreaks such as Foot and Mouth, SARS, Ebola and Coronavirus?