-
Notifications
You must be signed in to change notification settings - Fork 0
/
ex2_simple.R
65 lines (48 loc) · 2.13 KB
/
ex2_simple.R
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
library(slendr)
init_env()
# exercise 2 -- part a ----------------------------------------------------
simulate_afs <- function(Ne) {
# create a slendr model with a given N size
pop <- population("pop", N = Ne, time = 1)
model <- compile_model(pop, generation_time = 1, simulation_length = 100000)
# simulate a tree sequence
ts <-
msprime(model, sequence_length = 10e6, recombination_rate = 1e-8) %>%
ts_mutate(mutation_rate = 1e-8)
# get a random sample of names of 10 individuals
samples <- ts_names(ts) %>% sample(10)
# compute the AFS vector (dropping the 0-th element added by tskit)
afs <- ts_afs(ts, list(samples))[-1]
afs
}
afs_1 <- simulate_afs(1000)
afs_10 <- simulate_afs(10000)
afs_30 <- simulate_afs(30000)
afs_1 %>% plot(type = "o", main = "AFS, Ne = 1000")
afs_10 %>% plot(type = "o", main = "AFS, Ne = 10000")
afs_30 %>% plot(type = "o", main = "AFS, Ne = 30000")
# exercise 2 -- part b ----------------------------------------------------
afs_observed <- c(2520, 1449, 855, 622, 530, 446, 365, 334, 349, 244,
264, 218, 133, 173, 159, 142, 167, 129, 125, 143)
# we know that the Ne is between 1000 and 30000, so let's simulate
# a bunch of AFS vectors for different Ne values
afs_Ne1k <- simulate_afs(Ne = 1000)
afs_Ne5k <- simulate_afs(Ne = 5000)
afs_Ne6k <- simulate_afs(Ne = 6000)
afs_Ne10k <- simulate_afs(Ne = 10000)
afs_Ne20k <- simulate_afs(Ne = 20000)
afs_Ne30k <- simulate_afs(Ne = 30000)
# plot all simulated AFS vectors, highlighting the observed AFS in black
plot(afs_observed, type = "b", col = "black", lwd = 3,
xlab = "allele count bin", ylab = "count", ylim = c(0, 13000))
lines(afs_Ne1k, lwd = 2, col = "blue")
lines(afs_Ne5k, lwd = 2, col = "green")
lines(afs_Ne6k, lwd = 2, col = "pink")
lines(afs_Ne10k, lwd = 2, col = "purple")
lines(afs_Ne20k, lwd = 2, col = "orange")
lines(afs_Ne30k, lwd = 2, col = "cyan")
legend("topright",
legend = c("observed AFS", "Ne = 1000", "Ne = 5000",
"Ne = 6000", "Ne = 10000", "Ne = 20000", "Ne = 30000"),
fill = c("black", "blue", "green", "pink", "purple", "orange", "cyan"))
# the true Ne was 6543!