Nothing
library(knitr) opts_chunk$set(comment = NA, fig.width = 7, fig.height = 5, fig.align = "center")
For 20 samples from a panmictic population:
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")
For a neutral model with two populations and migration:
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:
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)
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!
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.