--- title: "R Notebook" output: html_notebook --- ## Illustrating the coalescent ```{r} options(tidyverse.quiet = TRUE) library(tidyverse) library(phyclust) library(ape) n_sample <- 100 ms_out <- ms(nsam = 2*n_sample, opts = "-T") ms_tree <- read.tree(text = ms_out) plot(ms_tree, show.tip.label = FALSE) ``` ## Comparing the efficiency of coalescent simulation with direct simulation Here's a simple comparison of simulating drift directly and simulating it with the coalescent. I've set \$N_e\$ to 100 and \$n_sample\$ to 25. ```{r} ## These are the parameters you can play with ## n_e <- 1000 n_sample <- 25 ## There's no easy way to tell when the forward simulation has run long enough ## to reach stationarity. I just fix the number of iterations to the expected ## fixation time and hope that that's close enough ## n_gen <- 4*n_e ## The number of times to repeat each simulation ## n_reps <- 10000 ## This is the forward simulation. ## start_time <- Sys.time() for (rep in 1:n_reps) { p <- numeric(n_gen) p[1] <- 0.5 for (i in 2:n_gen) { k <- rbinom(1, 2*n_e, p[i-1]) p[i] <- k/(2*n_e) } sample <- rbinom(1, 2*n_sample, p[n_gen]) } end_time <- Sys.time() cat("Elapsed time for forward simulation: ", round(end_time - start_time, 2)) ## Here's the coalescent simulation. ## start_time <- Sys.time() for (rep in 1:n_reps) { ms_out <- ms(nsam = 2*n_sample, opts = "-T") } end_time <- Sys.time() cat("\nElapsed time for coalescent simulation: ", round(end_time - start_time, 2)) ``` ## The coalescent with migration As you've seen in lab this week, `ms()` also allows you to simulate the coalescent with migration. What you didn't see was the actual call to `ms()` that does it. Here's an example with Wright's island model: 5 islands with a sample size of 25 from each island and \$N_em = 0.1\$. ```{r} ms_out <- ms(nsam = 125, # The total number of alleles in the sample nreps = 1, # The number of coalescent samples to produce opts = "-T -I 5 25 25 25 25 25 0.1", # -T produce a tree with coalescent times # -I migration follows Wright's island model # 5 populations # 25 the number of samples from each population # 0.1 Ne*m ) ms_tree <- read.tree(text = ms_out) plot(ms_tree, show.tip.label = FALSE) ```