R Notebook

Illustrating the coalescent

options(tidyverse.quiet = TRUE)
library(tidyverse)
library(phyclust)

Loading required package: ape

Attaching package: ‘ape’

The following object is masked from ‘package:dplyr’:

    where


Attaching package: ‘phyclust’

The following object is masked from ‘package:lubridate’:

    ms

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.

## 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))

Elapsed time for forward simulation:        1.72

## 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))

Elapsed time for coalescent simulation:  5.64

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\).

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)