---
title: "Molecular Clock Dating of Influenza H3N2"
author: "Erik Volz"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{H3N2}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# Introduction
`treedater` fits a strict or relaxed molecular clock to a phylogenetic tree and estimates evolutionary rates and times of common ancestry. The calendar time of each sample must be specified (possibly with bounds of uncertainty) and the length of the sequences used to estimate the tree. 

`treedater` uses heuristic search to optimise the TMRCAs of a phylogeny and the substitution rate. 
An uncorrelated relaxed molecular clock accounts for rate variation between lineages of the phylogeny which is parameterised using a Gamma-Poisson mixture model.

To cite:

* E.M. Volz and Frost, S.D.W. (2017) [Scalable relaxed clock phylogenetic dating](https://doi.org/10.1093/ve/vex025). Virus Evolution.

The most basic usage is 
```r
dater( tre, sts, s)
```

where 

* `tre` is an `ape::phylo` phylogeny, 
* `sts` is a named vector of sample times for each tip in `tre`
* `s` is the length of the genetic sequences used to estimate `tre`

## Invoking treedater from the command line

You can also use treedater from the command line without starting R using the `tdcl` script: 
```
./tdcl -h
Usage: ./tdcl [-[-help|h] [<logical>]] [-[-treefn|t] <character>] [-[-samplefn|s] <character>] [-[-sequenceLength|l] <double>] [-[-output|o] [<character>]]

-t <file> : file name of tree in newick format  
-s  <file> : should be a comma-separated-value file with sample times in format <taxon-id,sample-time> and no header
-l <length> :  the integer length of sequences in alignment used to construct the tree 
-o <file>: name of file for saving output 
```
Note that you may need to modify the first line of the `tdcl` script with the correct path to `Rscript` or `littler`.

# Influenza H3N2 HA 
This data set comprises 177 HA sequences collected over 35 years worldwide with known date of sampling.
We estimated a maximum likelihood tree using [iqtree](http://www.iqtree.org/). 
We will use the sample dates and ML tree to fit a molecular clock and estimate a dated phylogeny. 
First, load the tree (any method can be used to load a phylogeny into [ape::phylo format](http://ape-package.ird.fr/)): 

```{r}
require(treedater)
(tre <- ape::read.tree( system.file( 'extdata', 'flu_h3n2_final_small.treefile', package='treedater') ))
```
Note that this tree does not have a root, and in the process of fitting a molecular clock, we will estimate the best root location. 

```{r}
seqlen <- 1698 # the length of the HA sequences used to reconstruct the phylogeny
```

To fit the molecular clock, we will need the sample time for each lineage. 
Note that the date of sampling is incorporated into the name of each lineage, which is common in viral phylogenetics studies. 
The package includes a convenient function for extracting these dates:
```{r}
sts <- sampleYearsFromLabels( tre$tip.label, delimiter='_' )
head(sts)
```

How are samples distributed through time? 
```{r}
hist( sts , main = 'Time of sequence sampling') 
```

The basic usage of the treedater algorithm is as follows:
```{r}
dtr <- dater( tre , sts, seqlen, clock = 'strict' )
dtr
```
This produces a *rooted* tree with branches in calendar time. Note that if we invoked `dater` with a rooted input tree, it would _not_ estimate the root position. In this way, you can also set the root location in other ways, such as by using an outgroup. 
It is also good practice to provide at least one initial guess of the substitution rate using the `omega0` parameter, but if we omit that value as we have done here, treedater will attempt to guess good starting values. 
You can also specify an uncorrelated relaxed clock using `clock='uncorrelated'`. 

Lets see how long it takes to run treedater:

```{r}
rt0 <- Sys.time()
dtr <- dater( tre , sts, seqlen, clock = 'strict' )
rt1 <- Sys.time()
rt1 - rt0 
```
You can speed up treedater by providing a rooted tree, or by providing an educated guess of the substitution rate, or by using parallel computing with the `ncpu` option.

Note the returned value includes estimated substition rates and TMRCAs. 
The `dtr` object extends `ape::phylo`, so most of the methods that you can use in other R packages that use that format can also be used with a `dater` object. Lets plot the tree. 

```{r}
plot( dtr , no.mar=T, cex = .2 )
```

It looks like there are a couple of recent lineages that dont seem to fit well with the ladder-like topology. 
We can further examine this by doing a root-to-tip regression using the fitted tree and estimated node times which also shows a couple of outliers:
```{r}
rootToTipRegressionPlot( dtr )
```
It is always a good idea to visualize these distances to ensure that there is enough 'clock signal' in the data to reliably estimate rates and dates.  We will examine these outliers in the next section. 



## Detecting and removing outliers / testing for relaxed clock

To find lineages that dont fit the molecular clock model very well, run

```{r}
outliers <- outlierTips( dtr , alpha = 0.20) 
```

This returns a table in ascending order showing the quality of the molecular clock model fit for each lineage.
Now lineages could be selected for removal in various ways. 
Lets remove all tips that dont have a very high q-value :

```{r}
tre2 <- ape::drop.tip( tre, rownames(outliers[outliers$q < 0.20,]) )
```

Now we can rerun `dater` with the reduced tree:

```{r}
dtr2 <- dater(tre2, sts, seqlen, clock='uncorrelated', ncpu = 1)  # increase ncpu to use parallel computing
dtr2
```

After removing the outliers, the coefficient of variation of rates is much lower, suggesting that a strict clock model may be appropriate for the reduced tree. We can test the suitability of the strict clock with this test: 


```{r}
rct <- relaxedClockTest( tre2, sts, seqlen, ncpu = 1 ) # increase ncpu to use parallel computing
```

Note that the `ncpu` option enabled parallel computing to speed up this test. 

This test indicates a relaxed clock. 
Nevertheless, lets re-fit the model to the reduced tree using a strict clock for comparison:


```{r}
dtr3 <- dater( tre2, sts, seqlen, clock='strict' )
dtr3
```

```{r}
plot( dtr3 , no.mar=T, cex = .2 ) 
```

The rate is higher than the initial estimate with the relaxed clock and the recently-sampled outlying lineages have been removed. 

## Parametric bootstrap

Estimating confidence intervals for rates and dates is straightforward using a parametric bootstrap: 

```{r}
rt2 <- Sys.time()
(pb <- parboot( dtr3, ncpu = 1) )# increase ncpu to use parallel computing
rt3 <- Sys.time()
```

How fast was it? Note that the `ncpu` option would enable parallel computing. 

```{r}
rt3 - rt2 
```

We can also plot the estimated number of lineages through time with confidence intervals: 

```{r}
plot( pb ) 
```

If the *ggplot2* package is installed, we can use that instead:

```{r}
if ( suppressPackageStartupMessages( require(ggplot2)) )
  (pb.pl <- plot( pb , ggplot=TRUE) )
```

Note repeated bottlenecks and seasonal peaks of LTT corresponding to when samples are taken during seasonal epidemics.

The package also includes methods for nonparametric bootstrapping if you have already computed a bootstrap distribution of phylogenies. 

## Missing sample times
Suppose we only know some of the sample times to the nearest month, a common occurance in viral phylogenetic studies. 
To simulate this, we will put uncertainty bounds on some sample times equal to a +/- 2-week window. 
We create the following data frame with columns `lower` and `upper`: 

```{r}
sts.df <- data.frame( lower = sts[1:50] - 15/365, upper = sts[1:50] + 15/365 )
head(sts.df )
```

In this case, we constructed the data frame with bounds for the first 50 samples in the tree, but we could also manually construct a data frame for a few selected samples where times of sampling are uncertain, or for all of the samples. 

Now re-run treedater with the uncertain sample times. The vector `sts` provided here gives an initial guess of the unknown sample times. 

```{r}
(dtr4 <- dater( tre2, sts, seqlen, clock='strict', estimateSampleTimes = sts.df ) )
```

Note that the estimated rates and dates didnt change very much due to uncertain sample dates in this case.