00.Simulation.Rmd

---
######################################
# Click on "Run Document" in RStudio #
######################################
title: "Quick Introduction to Simulation"
author: "Miguel de Navascués"
output: learnr::tutorial
runtime: shiny_prerendered
bibliography: src/references.bib
---

These learning resources are available at [GitHub](https://github.com/mnavascues/ABCRFtutorial) & [Zenodo](http://doi.org/10.5281/zenodo.1435503).

## Setup (tl;dr)

This unit uses R packages `learnr` [@Schloerke2020] and `phyclust` [@Chen2011], function to read *ms* output from `readms.output.R` script distributed with `ms` [@Hudson2002], and functions to calculate summary statistics from *ms* simulations in file `ms_sumstats.R` are modified from [@Stadler2009a; @Stadler2009b].


```{r setup, include=TRUE}
library(learnr)
# function ms() (Hudson's coalescent simulator)
library(phyclust, quietly=T)
# function read.ms.output()
source("src/readms.output.R", local =T)
# functions to calculate summary statistics from ms output
source("src/ms_sumstats.R", local =T)
```

## What is a simulation?

Discuss:

- What is a computer simulation?
- What are the sources of randomness in a simulation?
- How are random numbers produced in a computer?

## Simulation of a coin toss

### Exercise 1

Simulate a ten coin tosses ($n=10$) using R function `sample(_____)` and probability of getting head $p$. Use a seed for random number generator so your simulation is reproducible:


```{r simulation_coin_toss, exercise=TRUE, exercise.lines=5}
n <- 10
p <- 0.5
sample(x=c("H","T"),size=_____, replace=_____, prob=_____)
```

```{r simulation_coin_toss-hint}
n <- 10
p <- 0.5
sample(x=c("H","T"),size=n, replace=T, prob=c(p,1-p))
```

```{r simulation_coin_toss-solution}
set.seed(123456)
n <- 10
p <- 0.5
sample(x=c("H","T"),size=n, replace=T, prob=c(p,1-p))
```

### Exercise 2

Simulate directly the number of heads ($h'$) using `rbinom(_____)`:

```{r simulation_coin_toss_rbinom, exercise=TRUE}
n <- 10
p <- 0.5
num_of_sims = 1
rbinom(n=_____,size=_____,prob=_____)
```

```{r simulation_coin_toss_rbinom-solution}
n <- 10
p <- 0.5
num_of_sims = 1
rbinom(n=num_of_sims,size=n,prob=p)
```



## Simulation of population genetic data

Simulation via the coalescent. For convenience, we will use the R package `phyclust` [@Chen2011] that includes a version of coalescent simulator `ms` in this tutorial. In practice, use of other simulators (e.g. msprime [@Kelleher2016], SLiM [@Haller2017], fastsimcoal [@Excoffier2011] or the original implementation of ms [@Hudson2002]) is recommended.

### Exercise 3

Simulate with `ms()` a sample of $n=10$ individuals sequenced for a mitochondrial DNA locus of lenght $l=10000$bp in a population of size $N=10000$ with mutation rate (per bp) $\mu=10^{-8}$:

```{r coalescent_sim, exercise=TRUE}
n<-10
l<-10000
N<-10000
u<-1e-8
num_of_sims<-1
_____
ms(nsam=_____,nreps=_____,opt=_____)
```

```{r coalescent_sim-hint}
n<-10
l<-10000
N<-10000
u<-1e-8
num_of_sims<-1
theta <- _____
ms(nsam=_____,nreps=_____,opt=paste("-t",theta))
```


```{r coalescent_sim-solution}
n<-10
l<-10000
N<-10000
u<-1e-8
num_of_sims<-1
theta <- 2*N*u*l 
ms(nsam=n,nreps=num_of_sims,opt=paste("-t",theta))
```

We can put this text into a more usable form (i.e. matrix of polymorphisms, etc) using a function distributed with the original software `ms` and that we have loaded in the setup with:

```{r eval=FALSE, echo=TRUE}
source("src/readms.output.R")
```

### Exercise 4

Use `read.ms.output()` function to read ms output and print only the genotype matrix (i.e. `gametes`).

```{r coalescent_sim_segsites, exercise=TRUE, exercise.lines=10}
n<-10
l<-10000
N<-10000
u<-1e-8
num_of_sims<-1
theta <- 2*N*u*l 
msout <- ms(nsam=n,nreps=num_of_sims,opt=paste("-t",theta))
sim_data <- read.ms.output(____)
_____
```

```{r coalescent_sim_segsites-hint}
n<-10
l<-10000
N<-10000
u<-1e-8
num_of_sims<-1
theta <- 2*N*u*l 
msout <- ms(nsam=n,nreps=num_of_sims,opt=paste("-t",theta))
sim_data <- read.ms.output(txt=msout)
sim_data$_____
```



```{r coalescent_sim_segsites-solution}
n<-10
l<-10000
N<-10000
u<-1e-8
num_of_sims<-1
theta <- 2*N*u*l
msout <- ms(nsam=n,nreps=num_of_sims,opt=paste("-t",theta))
sim_data <- read.ms.output(txt=msout)
sim_data$gametes
```


### Exercise 5

Calculate summary statistics (number of segregating sites, number of haplotypes, $\pi$ and Tajima's $D$) for the simulations. Functions `S()`, `NH()`, `PI()`, `TajimaD()` and `FuLiD()` have been loaded using:

```{r eval=FALSE, echo=TRUE}
source("src/ms_sumstats.R")
```


```{r summary_stats, exercise=TRUE, exercise.lines=15}
n<-10
l<-10000
N<-10000
u<-1e-8
num_of_sims <- 10
theta <- 2*N*u*l
msout <- ms(nsam=n,nreps=num_of_sims,opt=paste("-t",theta))
sim_data <- read.ms.output(txt=msout)

```

```{r summary_stats-solution}
n<-10
l<-10000
N<-10000
u<-1e-8
num_of_sims <- 10
theta <- 2*N*u*l
msout <- ms(nsam=n,nreps=num_of_sims,opt=paste("-t",theta))
sim_data <- read.ms.output(txt=msout)

S(sim_data)
NH(sim_data)
PI(sim_data)
TajimaD(sim_data)
FuLiD(sim_data)
```



#   

#### References