Coalescent

Concept

• Population genetics: characterization of the genetic variation within/between populatitons.

• We typically sample individuals from the population, and look at the variation in the samples. In molecular population genetics, we may sequence a couple genes from each individual, and look at the similarity/difference between them.

  

• We are interested in what demographic and evolutionary processes (e.g. population size change & genetic drift, mutation, recombination, migration, selection etc.) could have shaped the observed pattern of variation.

• The coalescent is a stochastic process useful for modeling the molecular evolution.
  Traditional way of modeling the stochastic process is chronological order (forward modeling, e.g., diffusion model).

  

• 2 alleles: Expected time of coalescence and average pairwise differences.
  1. Randomly pick 2 alleles from a population
  2. Go back in time until you find the most recent common ancestor (MRCA).
  3. Repeat this many times, and then take the average time to the MRCA.
  4. It is easy to calculate this expected time to the MRCA:

     expected time to the MRCA of 2 alleles $\approx$ 2 N

  5. If you compare two randomly sampled DNA sequences, how many differences ($\pi$) are there on average?

• m-coalescent
  • Expected time of coalescence among m samples $\approx$ 4 N

    

  • Expected total branch lengths: $4 N a_m$, where $a_m = \sum_{i=1}^{m-1} \frac{1}{i}$ (Watterson's coefficient).

  • Number of segregating sites:

    

    Expected number of segregating sites:

  • Tajima's D
    • Tajima's test statistics
      $$D \propto \pi_{obs} - (S_{obs}/a_m)$$ Note that the two elements are the estimaters of $\theta$: $D = 0$ expected under the neutral evolution.
    • It is a summary statistics of genealogical shape.

      

Coalescent Simulation

• We'll use Hudson's ms .

• Basic command:

  ms sampleSize numSimulations -t theta
  ms 4 2 -t 5.0
  

  The output looks like:

  ms 4 2 -t 5 
  7923 57589 4894
  
  //
  segsites: 8
  positions: 0.0323 0.0415 0.0706 0.1730 0.2127 0.3199 0.6417 0.8121 
  10001011
  10001110
  10011011
  01100000
  
  //
  segsites: 8
  positions: 0.1010 0.1347 0.2711 0.3833 0.4321 0.8748 0.9091 0.9095 
  00000011
  11111100
  00000011
  00000011
  

• -T option

  The simulation occurs in two steps: 1st, genealogy is created and 2nd, mutations are sprinkled on top of the genealogy. This option specifies to do only the first step, and print out the gene trees in Newick format.

  ms 4 2 -T
  

  The output is:

  ms 4 2 -T 
  27531 45558 8825
  
  //
  (3:0.129,(1:0.089,(2:0.012,4:0.012):0.077):0.040);
  
  //
  ((1:0.026,3:0.026):0.521,(2:0.532,4:0.532):0.015);
  

  The unit of branch lengths is 4 $N_0$ (current pop size).

• Seq-Gen can create the samples of DNA sequences from the gene tree.

  ms 4 2 -T | grep '(' > trees
  seq-gen -mHKY -l 40 -s .2 < trees > seqs
  

  Explanation of options:
  
  
  

  scale: Seq-gen considers that the number of expected substitution in a branch is given by the product of the branch length and the “scale” value.

  • In combination with the unit used in the output of ms, the option, -s .2, means that $\theta$ per nucleotide site is 0.2. See the following example to get this.

  • Example:
    • $N_0 = 10^6, \mu = 5 x 10^{-8} \Rightarrow \theta = 4 N_0 \mu = 0.2$
    • Let's say a branch is $4 x 10^6$ generations.
    • With the given mutation rate per site ($5 x 10^{-8}$), you expect $0.2 (=4 x 10^6 * 5 x 10^{-8})$ mutations per site.
    • The tree created by ms will say the branch lengh is 1.0 (because it uses the unit of $4 N_0$. So if you give the scaling of 0.2, it will give the correct mutation rate to realize $\theta = 0.2$.

• Other factors which can be incoporated into simulations.
  • Crossing-over, gene conversion
  • Demography (population size changes)
  • Population substructure and migration
  For these, you'll need to check the official documentation.

Exercises

1. Create a subroutine, which calculates the Watterson's coefficient for a given sample size, m.

2. Create a short perl script which extracts the number of segregating sites from the output of ms (from lines: segsites: 8 etc) and estimates the $\theta$ for each simulation. Then calculate the average of the $\theta$ from 1 million simulations. Confirm that the Watterson's estimater of $\theta = S_{obs} / a_m$ is pretty good.