Let’s look at one more example. Imagine the following five
demographic models describing the changes in the population size of some
abstract population:
These models and the population size change events could be,
depending on the specifics of our project, defined using different time
units. If you work with radiocarbon-dated samples and ancient DNA you
might be used to thinking about population history in units of “years
before present”—in this case, the “present” time (bottom of the figure)
would be at time 0 and earlier events would be specified as “5000 years
before present”, etc.
On the other hand, if you were interested in modelling some
theoretical population, you might not want to concern yourselves with
this and would be quite happy starting our simulation at “generation 1”
and continue until some specified “generation X”.
In the previous section of the vignette we have shown that SLiM and
msprime can be both used as simulation backend. In the remainder of this
vignette we will present additional analysis, demonstrating the fact
that you can express time of your models in a forward or backward
direction and still run those models in slendr’s SLiM (forward
simulator) or msprime (backward simulator).
Specifically, we will define all five demographic models from the
figure above twice, in different orientations of time – forward or
backward. We will then compute allele frequency spectrum from each
simulation (five models, each forward or backward) and verify that both
forward and backward pairs give the same result regardless of whether we
use the SLiM back end or the msprime backend.
Model definition
N <- 1000
N_factor <- 5 # by what factor should Ne change
seq_len <- 50e6
rec_rate <- 1e-8
mut_rate <- 1e-8
First, let’s define the five population histories in a forward time
direction:
# constant Ne model
forward_const <- population("const", time = 1, N = N)
# decreasing step Ne model
forward_decr <- population("decr", time = 1, N = N, map = FALSE) %>%
resize(time = 2000, N = N / N_factor, how = "step")
# increasing step Ne model
forward_incr <- population("inc", time = 1, N = N) %>%
resize(time = 2000, N = N * N_factor, how = "step")
# exponential increase in size
forward_exp_incr <- population("exp_inc", time = 1, N = N) %>%
resize(time = 2000, end = 3000, N = N * N_factor, how = "exponential")
# exponential decrease in size
forward_exp_decr <- population("exp_decr", time = 1, N = N) %>%
resize(time = 2000, end = 3000, N = N / N_factor, how = "exponential")
Now let’s flip the flow of time around and define the population
histories in units of “time ago”:
# constant Ne model
backward_const <- population("const", time = 5000, N = N)
# decreasing step Ne model
backward_decr <- population("decr", time = 5000, N = N) %>%
resize(time = 3000, N = N / N_factor, how = "step")
# increasing step Ne model
backward_incr <- population("inc", time = 5000, N = N) %>%
resize(time = 3000, N = N * N_factor, how = "step")
# exponential increase in size
backward_exp_incr <- population("exp_inc", time = 5000, N = N) %>%
resize(time = 3000, end = 2000, N = N * N_factor, how = "exponential")
# exponential decrease in size
backward_exp_decr <- population("exp_decr", time = 5000, N = N) %>%
resize(time = 3000, end = 2000, N = N / N_factor, how = "exponential")
For convenience, we will write a helper function which will, for a
given population history object (ten objects just above), compile a
slendr model, run it in both SLiM and msprime back
ends, load the tree sequence, and compute allele frequency spectrum
(AFS) from it. This will save us lots of code repetition (always a good
thing), minimize the chance of copy-pasting bugs (definitely a good
thing!) and compare the AFS patterns across different models and the two
slendr simulation back ends (they should be the same according
to population genetics theory).
compile_run_afs <- function(model_name, pop, seed = 42) {
# maximum length of the simulation (necessary for forward models which start
# in generation 1)
simulation_length <- 5000
# define sampling times given the direction of time
if (attr(pop, "history")[[1]]$time == 1) {
sampling_time <- simulation_length
direction <- "forward"
} else {
sampling_time <- 0
direction <- "backward"
}
# compile model
model <- compile_model(pop, generation_time = 15, direction = direction, simulation_length = simulation_length)
samples <- schedule_sampling(model, times = sampling_time, list(pop, 50))
# run the model in SLiM
ts_slim <- slim(model, sequence_length = seq_len, recombination_rate = rec_rate,
samples = samples, random_seed = seed, verbose = FALSE)
# run the same model in msprim
ts_msprime <- msprime(model, sequence_length = seq_len, recombination_rate = rec_rate,
samples = samples, random_seed = seed, verbose = FALSE)
# load the SLiM tree sequence
ts_slim <- ts_recapitate(ts_slim, Ne = N, recombination_rate = rec_rate, random_seed = seed) %>%
ts_mutate(mut_rate, random_seed = seed)
# load the msprime tree sequence
ts_msprime <- ts_mutate(ts_msprime, mut_rate, random_seed = seed)
# compute the AFS from the SLiM and msprime tree sequences and bind the
# results (derived allele counts per frequency bin) in a data frame
msprime_afs <- ts_afs(ts_msprime, polarised = TRUE)[-1]
slim_afs <- ts_afs(ts_slim, polarised = TRUE)[-1]
rbind(
data.frame(simulator = "msprime", model = model_name, f = msprime_afs),
data.frame(simulator = "SLiM", model = model_name, f = slim_afs)
) %>%
group_by(simulator, model) %>%
mutate(n = 1:n(), direction = direction) %>%
ungroup()
}
Now that this function doesn’t do anything special. It simply wraps a
standard slendr simulation pipeline (1. create population(s);
2. compile a model object; 3. simulate the model; 4. load the results;
5. analyse the data) into a single function.
Now that we have the whole pipeline, we can generate all of our
simulated data and bind all the result into a single data frame:
afs <- bind_rows(
compile_run_afs("constant", forward_const),
compile_run_afs("constant", backward_const),
compile_run_afs("step contraction", forward_decr),
compile_run_afs("step contraction", backward_decr),
compile_run_afs("step increase", forward_incr),
compile_run_afs("step increase", backward_incr),
compile_run_afs("exponential decrease", forward_exp_decr),
compile_run_afs("exponential decrease", backward_exp_decr),
compile_run_afs("exponential increase", forward_exp_incr),
compile_run_afs("exponential increase", backward_exp_incr)
) %>%
mutate(model = factor(
model, levels = c("step contraction", "constant", "step increase",
"exponential decrease", "exponential increase"))
)
Finally, we can plot the allele frequency spectra from all models,
for the two time directions, and for the two simulation back ends:
ggplot(afs, aes(n, f, color = direction, linetype = simulator)) +
geom_line(stat = "identity") +
facet_wrap(~ model) +
labs(x = "number of derived alleles", y = "frequency",
title = "Site frequency spectra obtained from five demographic models",
subtitle = "Each model was specified in forward or backward direction of time and executed by
two different backend scripts in slendr (SLiM and msprime)") +
guides(color = guide_legend("direction of\ntime in slendr"),
linetype = guide_legend("slendr backend\nengine used")) +
scale_linetype_manual(values = c(3, 2)) +
scale_x_continuous(breaks = c(1, seq(20, 100, 20)), limits = c(1, 100)) +
theme(legend.position = "bottom")
Again, a couple of things to note:
The two models that involve population contraction (be it in a
single step or an exponential collapse), we see a decrease in low
frequency variants. On the other hand, models with population expansion
show a high proportion of low frequency variants. The model with
constant population size shows intermediate pattern. This is a nice
validation that the simulation works as expected as it is exactly what
is predicted by population genetics theory.
Both forward and backward time models give exactly the same
result. This is not unexpected as anything else would be a bug in
slendr. However, although unsurprising, we hope that it makes
it clear that you are free to use whatever time unit specification you
want to encode slendr models.
For non-spatial models, results obtained from msprime
and SLiM back ends are consistent with one another. Again, this makes
sense because although they are completely different software,
population genetic theory governing the shape of allele frequency
spectra applies equally to both. In fact, the same example as we have
just implemented is an important part of slendr unit test suite
which verifies that the simulations created with slendr are
correct.
