17.dating_the_selection.Rmd

---
title: "Dating the selection"
output: github_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, fig.retina = 2)
library(tidyverse)
library(ggpubr)
library(cowplot)
library(rehh)
library(coda)
```

The software [startmrca](https://github.com/jhavsmith/startmrca) [@Smith2018-xt] is designed to estimate the timing of sweeps based on patterns of haplotype diversity and length that arise as a result of mutation and recombination after a sweep has occurred. 

`startmrca` relies on the concept of a "focal allele" which in-theory should represent the allele that is under selection, however, in-practise it need only represent an allele on the same haplotype as the selected allele in strong linkage disequilibrium with it. 

In order to estimate the timing of selection at a locus we need the following information;

  - The position of the focal allele
  - The identity of the focal allele at this position (0 or 1 for a biallelic SNP)
  - A list of individuals with haplotypes containing the focal allele
  - A list of individuals without the focal allele


#### Visualisation of haplotypes at the Haem Peroxidase Locus BLFC01000154_250000_300001

Our strongest signature of a selective sweep was identified at locus `BLFC01000154_250000_300001` in the inshore population.

```{r, results='hide'}
hap_names2pops <- function(hn){
  pops <- rep("IN",length(hn))
  sos <- grepl("RS[123]",hn)
  nos <- grepl("^AR",hn)
  ins <- (!sos & !nos)
  pops[sos]="SO"
  pops[nos]="NO"
  pops
}



vcf_file <- "data/hpc/startmrca/regions/BLFC01000154_250000_300001_aaref.vcf"

hh <- data2haplohh(hap_file = vcf_file,
                   polarize_vcf = FALSE,
                   vcf_reader = "data.table")

position <- 281245
marker_index <- which(positions(hh)==position)
```


Firstly we plot the haplotypes themselves in a region +/- 200 markers around the focal locus. We see that the highly homozygous haplotypes at the bottom all have the derived allele whereas few of those in the top have this allele. 


```{r, results='hide'}
hh_subset <- subset(hh, select.mrk = (marker_index-200):(marker_index+200))


plot(
  hh_subset,
  mrk = 201,
  group_by_allele = TRUE,
  ignore.distance = TRUE,
  col = c(NA, "red"),
  linecol = c("lightblue", "lightpink"),
  mrk.col = "black",
  cex = 0.1,
  cex.lab.hap = 0.2,
  hap.names = hap_names2pops(hap.names(hh))
)
```

A furcation plot shows there is one haplotype with very clear EHH.  This is clearly the candidate for selection that was identified using EHH statistics and the haplotypes under selection seem to segregate effectively based on the focal allele. 


```{r, results='hide'}
furcation <- calc_furcation(hh,mrk = marker_index)

hn <- hap.names(hh)

plot(furcation,hap.names = hap_names2pops(hn),cex=0.5,cex.lab=0.3)
```

The furcation plot is designed to capture relationships between haplotypes as a function of distance away from a focal point.  Another way to view the sequences is to consider the relationships between their core haplotype sequences.  This plot serves as a useful complement to the previous plot.  Here we can not only see that the focal allele accurately identifies the selected haplotype, but also .. that it completely identifies them (ie there are no selected haplotypes that do not have the focal allele).  We can see this because haplotypes with this allele form a single monophyletic clade.

```{r}
library(phytools)
library(ggtree)

hh_gg <- hh_subset@haplo %>% 
  as.data.frame() %>% 
  rownames_to_column("haplotype") %>% 
  pivot_longer(-haplotype,names_to = "position", values_to = "allele") %>% 
  extract(position,into = "pos",regex = "V([0-9]+)",convert = TRUE) 

njtree <- hh_subset@haplo %>% dist.gene() %>% nj %>% midpoint.root()

tree_data <- data.frame(tiplab=hn) %>% 
   mutate(pop=hap_names2pops(hn), has_marker_allele = (hh_subset@haplo[,201]==1))

gnjt <- ggtree(njtree,ladderize = FALSE) %<+% tree_data + 
#  geom_tiplab(aes(label=pop,color=pop),align = TRUE,size=2) +
  geom_tippoint(aes(color=pop),alpha=0.9,size=1) + 
  ggsci::scale_color_startrek() + 
  theme(legend.position = "none")

ggsave(gnjt,filename = "figures/hap_tree_peroxidase.png",width = 3, height = 5)
```


```{r}
gnjt + geom_facet(panel="Haplotypes",data = hh_gg %>% filter(allele==1), 
                geom= geom_point,mapping = aes(x=pos),color="blue",shape="|",size=1) 
```

```{r}
# Identify haplotypes in inshore that are not the selected hap.  This is used later for GEVA 

pop="inshore"

haps_alleles <- data.frame(allele=haplo(hh)[,marker_index]) %>% 
  rownames_to_column("hap") %>% 
  extract(col = hap,into = "loc_code","([^_]+)", remove = FALSE) %>% 
  mutate(location = case_when(
    loc_code %in% c("BR","AI") ~ "inshore",
    loc_code %in% c("AR") ~ "northoffshore",
    loc_code %in% c("RS1","RS2","RS3") ~ "southoffshore"
  ))

haps_alleles %>% filter(allele==0) %>% filter(location=="inshore")
```


## Running STARTMRCA

Our goal was to estimate the timing of selection on the gene `s0150.g24` which contained 75 sites whose allele frequencies and make them potential candidates for the beneficial allele.  We chose 3 of these spanning the width of the gene and then ran startmrca on each as follows;

- We used `vcftools` to extract a phased vcf file across a 1Mb region centered on `s0150.g24`. This extraction was performed based on the vcf file where the ancestral allele has been assigned to the reference.
- We then ran startmrca using the script [03_run_tmrca.sh](data/hpc/startmrca/03_run_tmrca.sh) which itself calls [run_startmrca.R](data/hpc/startmrca/run_startmrca.R). 

For each of our three loci we check convergence using the Gelman diagnostic.  Then report the median estimate for the TMRCA.  Dates shown are based on a mutation rate of 1.2e-8 and a recombination rate of 3.2e-8.  This mutation rate is the same as used in other analyses (SMC++, fastSimCoal) and the recombination rate is based on a linkage map for *Acropora millepora* published by Dixon et al [@Dixon2015-x] which had a length of 1358 centimorgans (cM). Assuming a genome size of 430mb for *A. millepora* this corresponds to a mean recombination rate of 3.16e-8 per base per generation, which we rounded to 3.2e-8. 


```{r, include=FALSE}
read_chain <- function(path){
  tmp <- read_rds(path) 
  tmp$t.chain %>% 
    as.data.frame() %>% 
    mutate(p = path) %>% 
    extract(p,into=c("locus","rep"),regex="(B.*)_([0-9]+).rds") %>% 
    filter(V4==1) %>% 
    slice_tail(n=5000) %>% 
    slice_sample(n=1000)
}
if ( !file.exists("data/r_essentials/17_tmrca_data.rds")){
  tmrca_data <- list.files("data/hpc/startmrca/results/","*.rds", full.names = TRUE) %>% 
    map_dfr(read_chain) 
  write_rds(tmrca_data,"data/r_essentials/17_tmrca_data.rds")
} else {
  tmrca_data <- read_rds("data/r_essentials/17_tmrca_data.rds")
}



to_mcmc_list <- function(locus_data){
  locus_data %>% dplyr::select(rep,time=V1) %>% group_by(rep) %>% group_map( ~ mcmc(.x)) %>% mcmc.list()
}



locus_groups <- split(tmrca_data,tmrca_data$locus)
locus_mcmc_lists <- locus_groups %>% lapply(to_mcmc_list)
plot(locus_mcmc_lists[[1]])
gelman.diag(locus_mcmc_lists[[1]])

gelman_pe <- function(x){
  if( length(x)>1){
    gd <- gelman.diag(x)
    return(gd$psrf)
  }else{
    return(c(10,10))
  }
}

gdr <- locus_mcmc_lists %>% sapply(gelman_pe) %>% t() %>% as.data.frame() %>% rownames_to_column("locus")

converged <- gdr %>% filter(V2<2) %>% pull(locus)
```


```{r}
tmrca_data %>% 
  filter(locus %in% converged) %>% 
  group_by(locus) %>% 
  dplyr::select(time=V1,rep) %>% 
  summarise(median(time)*5, nchains=length(unique(rep))) %>% knitr::kable()
```