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") 

☛ When working with graphics in R on xfce you may find it useful to move the Terminal and the R Graphics Device to a separate Workspace by clicking the top left hand corner of the window and selecting Move to Another Workspace. Shift between workspaces using Ctrl Alt Left and Ctrl Alt Right.

#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 the command line by clicking Applications, Terminal Emulator and change directory to your Downloads folders. The following command will launch PhyML

cd ~/Downloads
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 35 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

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.12, indicating that on average 12 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 ~/GitHub/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.