-
Notifications
You must be signed in to change notification settings - Fork 4
/
coala-examples.Rmd
77 lines (63 loc) · 1.98 KB
/
coala-examples.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
---
title: "Examples"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Examples}
%\VignetteEngine{knitr::rmarkdown}
%\usepackage[utf8]{inputenc}
---
```{r echo=FALSE}
library(knitr)
opts_chunk$set(comment = NA,
fig.width = 7,
fig.height = 5,
fig.align = "center")
```
## 1. Obtain the neutral distribution of Tajima's D
For 20 samples from a panmictic population:
```{r}
library(coala)
model <- coal_model(20, 2000) +
feat_mutation(2) +
feat_recombination(1) +
sumstat_tajimas_d()
stats <- simulate(model, seed = 15)
plot(density(stats$tajimas_d, na.rm = TRUE),
main = "Neutral Distribution of Tajiam's D")
```
## 2. Generate the site frequency spectrum
For a neutral model with two populations and migration:
```{r}
model2a <- coal_model(c(10, 10), 100) +
feat_mutation(10) +
feat_recombination(5) +
feat_migration(0.5, symmetric = TRUE) +
sumstat_sfs(population = "all")
stats <- simulate(model2a, seed = 20)
barplot(stats$sfs / sum(stats$sfs),
names.arg = seq_along(stats$sfs),
col = 3)
```
And again, but now with a bottleneck in one population:
```{r}
model2b <- model2a +
feat_size_change(0.1, population = 2, time = 0.25) +
feat_size_change(1, population = 2, time = 0.5)
stats <- simulate(model2b, seed = 25)
barplot(stats$sfs / sum(stats$sfs),
names.arg = seq_along(stats$sfs),
col = 4)
```
## 3. Plot the nucleotide diversity against the mutation rate:
```{r}
model3 <- coal_model(10, 50) +
feat_mutation(par_prior("theta", sample.int(100, 1))) +
sumstat_nucleotide_div()
stats <- simulate(model3, nsim = 40)
mean_pi <- sapply(stats, function(x) mean(x$pi))
theta <- sapply(stats, function(x) x$pars[["theta"]])
plot(theta, mean_pi, pch = 19, col = "orange")
abline(lm(mean_pi ~ theta), col = "red2", lty = 3)
```
If you have a nice example for using coala, feel free to extend this
vignette via a pull request [on GitHub](https://github.com/statgenlmu/coala)!