library(tidyverse)
library(parallel)
library(MASS)
set.seed(42)Simulations in population genetics
Many problems in population genetics cannot be solved by a mathematician, no matter how gifted. [It] is already clear that computer methods are very powerful. This is good. It […] permits people with limited mathematical knowledge to work on important problems […]
Exploring expected behavior of statistics
Fitting model parameters (i.e. ABC)
Ground truth for testing and method development
\(N\) = 500, \(p_0 = 0.5\)
Code
plot(p_traj, type = "l", ylim = c(0, 1),
xlab = "generations", ylab = "allele frequency")
abline(h = p_start, lty = 2, col = "red")
\(N\) = 500, \(p_0 = 0.5\) (20 replicates)
Code
reps <- replicate(20, simulate(N = 500, p_start = 0.5, generations = 2000))
matplot(reps, ylim = c(0, 1), xlab = "generations", ylab = "allele frequency", type = "l", lty = 1)
\(N\) = 1000, \(p_0 = 0.5\) (20 replicates)
Code
reps <- replicate(20, simulate(N = 1000, p_start = 0.5, generations = 2000))
matplot(reps, ylim = c(0, 1), xlab = "generations", ylab = "allele frequency", type = "l", lty = 1)
\(N\) = 5000, \(p_0 = 0.5\) (20 replicates)
Code
reps <- replicate(20, simulate(N = 5000, p_start = 0.5, generations = 2000))
matplot(reps, ylim = c(0, 1), xlab = "generations", ylab = "allele frequency", type = "l", lty = 1)
\(N\) = 10000, \(p_0 = 0.5\) (20 replicates)
Code
reps <- replicate(20, simulate(N = 10000, p_start = 0.5, generations = 20000))
matplot(reps, ylim = c(0, 1), xlab = "generations", ylab = "allele frequency", type = "l", lty = 1)
\(N\) = 10000, \(p_0 = 0.5\) (20 replicates)
Code
factors <- fractions(c(3, 2, 1, 1/2, 1/5, 1/10))
matplot(reps, ylim = c(0, 1), xlab = "generations", ylab = "allele frequency", type = "l", lty = 1)
abline(v = 10000 * factors, lwd = 5)
Expected allele frequency distribution
Code
p_start <- 0.5
N <- 10000
factors <- fractions(c(3, 2, 1, 1/2, 1/5, 1/10))
final_frequencies <- parallel::mclapply(
seq_along(factors), function(i) {
f <- factors[i]
t <- as.integer(N * f)
# get complete trajectories as a matrix (each column a single trajectory)
reps <- replicate(10000, simulate(N = N, p_start = p_start, generations = t))
# only keep the last slice of the matrix with the final frequencies
data.frame(
t = sprintf("t = %s * N", f),
freq = reps[t, ]
)
},
mc.cores = parallel::detectCores()
) %>% do.call(rbind, .)
final_frequencies$t <- fct_rev(fct_relevel(final_frequencies$t, sprintf("t = %s * N", factors)))
final_frequencies %>% .[.$freq > 0 & .$freq < 1, ] %>%
ggplot() +
geom_histogram(aes(freq, y = ..density.., fill = t), position = "identity", bins = 100, alpha = 0.75) +
labs(x = "allele frequency") +
coord_cartesian(ylim = c(0, 3)) +
facet_grid(t ~ .) +
guides(fill = guide_legend(sprintf("time since\nthe start\n[assuming\nN = %s]", N))) +
theme_minimal() +
theme(strip.text.y = element_blank(),
axis.title.x = element_text(size = 15),
axis.title.y = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15),
legend.title = element_text(size = 15),
legend.text = element_text(size = 15))
What is SLiM?
- Forward-time simulator
- It’s fully programmable!
- Massive library of functions for:
- Demographic events
- Various mating systems
- Selection, quantitative traits, …
- > 700 pages long manual!
SLiMgui – IDE for SLiM
What is msprime?
What is msprime?
- A Python module for writing coalescent simulations
- Extremely fast (genome-scale, population-scale data)
- You must know Python fairly well to build complex models
www.slendr.net
Why a new package? – spatial simulations!
Workaround for an RStudio bug
RStudio sometimes interferes with Python setup that we need for simulation. To fix this, go to Tools -> Global Options in your RStudio and set the following options:
Model visualization
plot_model(model)
What is a tree sequence?
- A record of full genetic ancestry of a set of samples
- An encoding of DNA sequence carried by those samples
- An efficient analysis framework
What we usually have
What we usually want
(As full as possible) a representation of our samples’ history:
How does this work?!
Tree-sequence tables (tskit docs)
A tree (sequence) can be represented by tables of:
- nodes (“chromosomes”)
- edges (branches) between nodes
- mutations on edges
- individuals (who carry nodes)
- populations
A set of such tables is a tree sequence.
Tree-sequence tables in practice
Code
ts_nodes(ts) %>%
as_tibble() %>%
dplyr::select(name, pop, node_id, time) %>% head(3)# A tibble: 3 × 4
name pop node_id time
<chr> <fct> <int> <dbl>
1 pop_1 pop 0 0
2 pop_1 pop 1 0
3 pop_2 pop 2 0Code
ts_edges(ts) %>%
as_tibble() %>%
dplyr::select(child_node_id, parent_node_id) %>% head(3)# A tibble: 3 × 2
child_node_id parent_node_id
<int> <int>
1 0 43032
2 0 48715
3 0 54705Let’s get back to the model we defined earlier
R interface to tskit
This R interface links to Python methods implemented in tskit.
Here is the magic
Tree sequences make it possible to directly compute many quantities of interest without going via conversion to a genotype table/VCF!
Extracting sample information
Each “sampled” individual in slendr has a symbolic name, a sampling time, and a population assignment:
Example: allele frequency spectrum
# sample 5 individuals
# (i.e. 10 chromosomes)
samples <-
ts_samples(ts) %>%
sample_n(5) %>%
pull(name)
# compute allele frequency
# spectrum from the given set
# of individuals
afs1 <- ts_afs(
ts, list(samples),
mode = "branch",
polarised = TRUE,
span_normalise = TRUE
)
afs1 [1] 40188.717 20286.332 13275.749 10084.000 8204.045 6736.871 5629.087
[8] 5140.529 4426.037 3913.154
plot(afs1, type = "b",
xlab = "allele count bin",
ylab = "frequency")
How can we compute statistics?
There is a duality between mutations and branch lengths in trees (more here).
But what if we want mutations?
Coalescent and mutation processes can be decoupled!
This means we can add mutations
after the simulation.
Example: allele frequency spectrum
# sample 5 individuals
# (i.e. 10 chromosomes)
samples <-
ts_samples(ts) %>%
sample_n(5) %>%
pull(name)
# compute allele frequency
# spectrum from the given set
# of individuals
afs2 <- ts_afs(
ts, list(samples),
polarised = TRUE
)
afs2 [1] 40822 19806 13012 9995 7906 6712 5604 5142 4540 3830
plot(afs2, type = "b",
xlab = "allele count bin",
ylab = "frequency")
lines(afs1, type = "b", col = "red")
legend("topright", legend = c("mutation-based AFS", "theoretical branch-based AFS"),
fill = c("black", "red"))
Workaround for an RStudio bug
RStudio sometimes interferes with Python setup that we need for simulation. To fix this, go to Tools -> Global Options in your RStudio and set the following options:
Exercise #2
The squirrels have split into three species with different demographic histories (s1, s2, s3 – model on the right). Species s2 and s3 are daughter populations which split in generation 2 from the original species s1.
Help the Oracle predict the future shape of the AFS of the squirrels after 10000 generations, assuming starting \(N_e\) = 6543? Species s1 will remain constant, species s2 will expand 3X, species s3 will get 3X smaller.
Let’s build another toy model
Code
o <- population("o", time = 1, N = 100)
c <- population("c", time = 2500, N = 500, parent = o)
a <- population("a", time = 3000, N = 2000, parent = c)
b <- population("b", time = 3500, N = 4000, parent = a)
x1 <- population("x1", time = 3800, N = 8000, parent = c)
x2 <- population("x2", time = 4000, N = 10000, parent = x1)
gf <- gene_flow(from = b, to = x1, start = 5500, end = 6000, rate = 0.3)
model <- compile_model(
populations = list(a, b, x1, x2, c, o), gene_flow = gf,
generation_time = 1, simulation_length = 7000
)
ts <- msprime(model, sequence_length = 100e6, recombination_rate = 1e-8) %>%
ts_mutate(mutation_rate = 1e-8)
plot_model(model, proportions = TRUE)
Program the following model in slendr
Diversity and divergence
Use ts_diversity() to compute diversity in each simulated population (CH, AFR, …).
Does the result match the expectation given demographic history?
What about divergence between all pairs of populations (ts_divergence())? Do your results recapitulate phylogenetic relationships? How about ts_f3() or ts_fst()?
Exercise #4 – solution (divergence)
Code
divergence <-
ts_divergence(ts, sample_names) %>%
arrange(divergence)
divergence# A tibble: 10 × 3
x y divergence
<chr> <chr> <dbl>
1 ALT VI 0.000100
2 AFR EUR 0.000465
3 EUR VI 0.000847
4 ALT EUR 0.000849
5 AFR VI 0.000871
6 AFR ALT 0.000872
7 AFR CH 0.00427
8 CH VI 0.00427
9 ALT CH 0.00427
10 CH EUR 0.00427 Code
divergence %>%
mutate(pair = paste(x, "-", y),
pair = fct_reorder(pair, divergence)) %>%
ggplot(aes(pair, divergence)) +
geom_point() +
labs(x = "", y = "divergence") +
ggtitle("Pairwise genetic divergence between populations") +
theme_minimal()
Computing \(f_4\) statistic – ts_f4()
Having simulated data from this model, can we detect gene flow from b into x1?
Computing \(f_4\) statistic for all x individuals
Code
ggplot(x_inds) +
geom_rect(aes(xmin = 5500, xmax = 6000, ymin = -Inf, ymax = Inf), alpha = 0.5, fill = "gray") +
geom_jitter(aes(time, f4, color = pop)) +
geom_line(data = . %>% group_by(pop, time) %>% summarise(f4 = mean(f4)),
aes(time, f4, color = pop), size = 2) +
geom_hline(yintercept = 0, linetype = 2, size = 1) +
theme_minimal()
Sometimes a tree sequence is too big
Extracting a tree #1734
tree <- ts_phylo(ts_small, 1734, quiet = TRUE)Code
# this needs another R package
# simply using the ape package and running `plot(tree)` will also work
library(ggtree)
nodes <- ts_nodes(tree) %>%
as_tibble %>%
dplyr::select(node = phylo_id, pop)
p_tree <- ggtree(tree) %<+% nodes +
geom_tiplab(aes(color = pop, fill = pop)) +
guides(color = "none") +
scale_x_continuous(limits = c(-7e6, 900e3)) +
geom_vline(xintercept = -700e3, linetype = 2, color = "red")
revts(p_tree)
Extracting a tree #1
tree <- ts_phylo(ts_small, 1, quiet = TRUE)Code
# this needs another R package
# simply using the ape package and running `plot(tree)` will also work
library(ggtree)
nodes <- ts_nodes(tree) %>%
as_tibble %>%
dplyr::select(node = phylo_id, pop)
p_tree <- ggtree(tree) %<+% nodes +
geom_tiplab(aes(color = pop, fill = pop)) +
guides(color = "none") +
scale_x_continuous(limits = c(-7e6, 900e3)) +
geom_vline(xintercept = -700e3, linetype = 2, color = "red")
revts(p_tree)
Extracting a tree #903
tree <- ts_phylo(ts_small, 903, quiet = TRUE)Code
# this needs another R package
# simply using the ape package and running `plot(tree)` will also work
library(ggtree)
nodes <- ts_nodes(tree) %>%
as_tibble %>%
dplyr::select(node = phylo_id, pop)
p_tree <- ggtree(tree) %<+% nodes +
geom_tiplab(aes(color = pop, fill = pop)) +
guides(color = "none") +
scale_x_continuous(limits = c(-7e6, 900e3)) +
geom_vline(xintercept = -700e3, linetype = 2, color = "red")
revts(p_tree)