Exploring the coalescent

We use ms() from phyclust to simulate samples from the coalescent and various functions from ape to display the coalescent tree and calculate statistics from it. Let’s start with a really simple example of constructing a coalescent sample of size 10. The nsam argument is the number of samples, and the opts = "-T" argument tell ms() to return a coalescent tree. We then convert it, using read.tree() to a format that we can plot()

The coalescent in R

library(ape)
library(phyclust)

rm(list = ls())

## set random number seed so that same tree is produced every time this code
## is run
##
set.seed(1234)

ms_out <- ms(nsam = 10, opts = "-T")
ms_tree_1 <- read.tree(text = ms_out)
plot(ms_tree_1)

You’ll notice that the tip labels run from s1 through s10. Each corresponds to a different allele. That was pretty easy, but we can do a bit more. We can ask how “deep” the tree is, i.e., how far back in time we have to go to find the most recent common ancestor of all the alleles in our sample.

Branching times

branching.times(ms_tree_1)[1]

      11 
1.112965 

You’re probably wondering about the [1] after branching.times(). Well, branching.times() returns a vector of all of the branching times in the tree.

branching.times(ms_tree_1)

        11         12         13         14         15         16         17         18         19 
1.11296475 0.79145014 0.25495386 0.03760511 0.15911412 0.03034905 0.09305522 0.05205787 0.02415735 

The numbers 11-19 are the node numbers in the tree above.

plot(ms_tree_1)
nodelabels()

As you can see, node 11 (which is the first one in the list), is the deepest node, i.e., the most recent common ancestor of all of the alleles in our sample.

Scaling of branching times

The other thing you may be wondering is how to interpret the number 1.1129647 that corresponds to the time of the most recent common ancestor.¹ To explain that I have to point out something that you may have missed when we called ms(): We didn’t specify a population size. How can we get away with that, since I told you that drift always depends on the population size?

Well, remember what you know² about the coalesent process. The expected time to coalescence of all alleles is \(4N_e\). ms() reports branching time in “coalescent units” where 1 unit equals \(4N_e\). So if branching.times(ms_tree_1)[1] is 1.1129647 and the effective size of our population is 25, then the time (number of generations before the present) to the most recent common ancestor is 111.2964749. Notice that the time isn’t an integer the way you’d expect it to be. That’s because of an approximation that is generally used in implementing the coalescent in software. If you want to turn it into an integer, you can do this:

round(4*25*branching.times(ms_tree_1)[1])

 11 
111 

Whether you round it or not, you’ll see that the most recent common ancestor of the alleles in our sample didn’t occur exactly 100 generations ago as we’d expect. If, however, we repeat that simulation a bunch of times, you can see that we come pretty close - on average.

library(ggplot2)

N_e = 25

c_time <- numeric(0)
for (i in 1:1000) {
  ms_out <- ms(nsam = N_e, opts = "-T")
  ms_tree <- read.tree(text = ms_out)
  c_time[i] <- 4*N_e*branching.times(ms_tree)[1]
}
p <- ggplot(data.frame(t = c_time), aes(x = t)) +
  geom_histogram(bins = 50, alpha = 0.6) +
  geom_vline(xintercept = 4*N_e, linetype = "dashed", color = "red") + 
  theme_bw()
p

We expect a coalescence time of 100 and our observed average coalescence time is 95.2170163. Not too bad. Notice that the standard deviation of the coalescence time is 50.9947031. That’s a lot of variability in coalescence time, which re-emphasizes a point I’ve made many times: There is a lot of variability associated with drift. While the expectations we focus on provide a useful point of comparison, they can delude us into thinking that we know more than we actually do.

The coalescent with migration

What happens when we have a subdivided population so that we have several populations exchanging genes? Here’s how we simulate a finite island model with three populations, local population sizes of 10 in each population, and \(N_em = 2\).³

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

get_pop <- function(tree, sample_sizes, k) {
  idx <- as.numeric(gsub("^s(.*)$", "\\1", tree$tip.label))
  cum <- 0
  for (i in 1:length(sample_sizes)) {
    if (idx[k] < cum + sample_sizes[i] + 1) {
      retval <- i
      break
    } else {
      cum <- cum + sample_sizes[i]
    }
  }
  return(i)
}

tipcolors <- function(tree, sample_sizes) {
  my_colors <- brewer.pal(length(sample_sizes), "Set1")
  colors <- numeric(length(sample_sizes))
  for (i in 1:sum(sample_sizes)) {
    colors[i] <- my_colors[get_pop(tree, sample_sizes, i)]
  }
  return(colors)
}

## Note: sample sizes should be multiplied by 2 before calling
## island for diploid samples
##
island <- function(sample_sizes, ne_m) {
  opt_string <- paste("-T -I ", length(sample_sizes), " ", sep = "")
  for (i in 1:length(sample_sizes)) {
    opt_string <- paste(opt_string, sample_sizes[i], " ", sep = "")
  }
  opt_string <- paste(opt_string, ne_m, "\n", sep = "")
  ms_out <- ms(nsam = sum(sample_sizes), nreps = 1, opts = opt_string)
  return(ms_out)  
}

plot_island <- function(sample_sizes, ne_m) {
  ms_out <- island(2*sample_sizes, ne_m)
  ms_tree <- read.tree(text = ms_out)
  plot(ms_tree, 
       show.tip.label = FALSE)
  tiplabels(pch = 19, cex = 0.75, 
            col = tipcolors(ms_tree, 
                            2*sample_sizes))
  return(branching.times(ms_tree)[1])
}

time <- plot_island(sample_sizes = c(5, 5, 5), ne_m = 10)

The colors correspond to the population from which each sample was collected. plot_island() returns the scaled time to the most recent common ancestor, 10.5703715. We are implicitly assuming that the effective population size is the same in each local population.⁴ The scaling is in terms of the product of the local effective population size and the number of populations, i.e., 1 coalescent unit = \(4N_ek\), where \(k\) is the number of populations.

Notice how different the result is if I set \(N_em = 0.5\) instead of \(N_em = 2\).

time <- plot_island(sample_sizes = c(5, 5, 5), ne_m = 0.1)

The alleles are now nicely arranged into populations. Does it surprise you that the populations in which alleles occur form monophyletic groups when \(N_em \ll 1\) and that they don’t when \(N_em \gg 1\).⁵

Lab 7

You will be performing another set of simulations this week to compare coalescence times in two scenarios:

1. A finite island model of migration.
2. A one-dimensional stepping stone model of migration.

You saw both the finite island model and the one-dimensional stepping stone model last week. As a reminder, in a one-dimensional stepping stone model the populations are arranged in a lineary array and only those populations that are adjacent to one another can exchange migrants. If \(m\) is the migration rate, i.e., the fraction of a population composed of migrants then the fraction of population \(i\) derived from population \(i-1\) is \(\frac{m}{2}\) and the fraction of population \(i\) derived from population \(i+1\) is \(\frac{m}{2}\).⁶

Use the functions below to explore the difference in coalescence times for finite island models and one-dimensional stepping stone models of migration, and answer the following questions:

1. How does the number of populations included affect the avearge time to coalescence of all alleles? Try at least one case in which the number of populations is relatively small, i.e., 10 or fewer, and one in which it is moderately large, i.e., 25 or more.

2. How does the migration rate affect the average time to coalescence of all alleles? Does it depend on whether migration is better approximated by a finite island model or a one-dimensional stepping stone model?

3. Are there consistent differences between the average time to coalescence in a finite island model and a one-dimensional stepping stone model? If there are, what might explain the differences that you see?

Hints

run_simulations() will automatically keep the number of populations in both models the same to ensure that you can always make appropriate comparisons of the models. It will also set the local effective population size to 25, and it sets the number of replications to 1000 by default. You might want to run the function once with n_reps = 100 to get a feel for how long it will take the simulation to run on your computer before you begin the full-scale simulations

run_simulations() returns a data frame⁷ with two columns: Model and Time. If you’re interested, you can use that to examine the median time to coalescence (rather than the mean) and the variance (or standard deviation) in coalescence times. You’ll notice that run_simulation() automatically prints out the mean coalesence times for you.

In answering question #3 keep in mind that your answer might depend on how many populations there are.⁸

## one-dimensional stepping stone
##
one_d_stepping_stone <- function(n_e, m, n_pops) {
  stopifnot((m > 0) && (m < 1))
  m_matrix <- diag(1 - m, nrow = n_pops)
  m_matrix[1, 2] <- m
  m_matrix[n_pops, n_pops - 1] <- m
  for (i in 2:(n_pops - 1)) {
    m_matrix[i, i - 1] <- m/2
    m_matrix[i, i + 1] <- m/2
  }
  opt_string <- paste("-T -I ", n_pops, " ", sep = "")
  for (i in 1:n_pops) {
    opt_string <- paste(opt_string, 2*n_e, " ", sep = "")
  }
  opt_string <- paste(opt_string, "-ma", sep = "")
  ## transpose m_matrix before converting to vector to get vector from matrix
  ## by row instead of by column
  ##
  m_vector <- 4*n_e*c(t(m_matrix))
  for ( i in 1:length(m_vector)) {
    opt_string <- paste(opt_string, " ", m_vector[i], sep = "")
  }
  ms_out <- ms(nsam = 2*n_e*n_pops, nreps = 1, opts = opt_string)
  return(ms_out)
}

run_simulation <- function(n_pops, m, n_reps = 1000) {
  N_e = 25
  pop_sizes <- rep(N_e, n_pops)
  df <- tibble(Model = NA, Time = NA)
  for (i in 1:n_reps) {
    ms_out <- island(pop_sizes, N_e*m)
    ms_tree <- read.tree(text = ms_out)
    df <- add_row(df,
                  Model = "Island",
                  Time = (4*N_e*n_pops)*branching.times(ms_tree)[1])
    ms_out <- one_d_stepping_stone(N_e, m, n_pops)
    ms_tree <- read.tree(text = ms_out)
    df <- add_row(df,
                  Model = "1-d",
                  Time = (4*N_e*n_pops)*branching.times(ms_tree)[1])
  }
  df <- filter(df, !is.na(Model))
  p <- ggplot(df, aes(x = Time, fill = Model)) +
    geom_histogram(position = "identity", bins = n_reps/10, alpha = 0.4) +
    ggtitle(paste("n_pops = ", n_pops, ", m = ", m)) +
    theme_bw()
  print(p)
  sum <- df %>% group_by(Model) %>%
    summarize(Mean_Time = mean(Time))
  print(sum)
  return(df)
}

---

1. Remember that when using the coalescent, time runs backwards. So the time of the most recent common ancestor is a time in the past.↩︎

2. or what I told you↩︎

3. Note: You’ll need to install RColorBrewer if you want to run this code.↩︎

4. There is a way to relax this assumption, but we won’t take advantage of it here.↩︎

5. Note: Read \(\ll\) as “much less than” and \(\gg\) as “much greater than.”↩︎

6. Note: For the first population in the array, i.e., population 1, the fraction of the population derived from population 2 is \(m\). Similarly, if there are \(k\) populations in the arry, the fraction of population \(k\) that is derived from population \(k-1\) is \(m\).↩︎

7. Well, actually a tibble↩︎

8. Why else would I have asked you to look both at one example in which tne number of populations is 10 or fewer and another in which the number of populations is 25 or more?↩︎