Example: Drive other programs

• Simple divergence model.

  

• Goal: Test if symmetric migration model or unidirectional model is better.
• Strategy:
  1. Choose a random number (0.01 - 1) for T1
  2. Run 1 ms after flipping a coin:
     head -> Simulate Model A

     ms 8 1 -t 8 -I 2 3 5 -ma x 8.0 16.0 x  -n 1 0.5 -n 2 0.25  -en $T1 1 1.0  -ej $T1 2 1
     

     tail -> Simulate Model B
     Use the migration matrix:

     -ma x 8.0 0.0 x
     

  3. Calculate the summary statistics (pi, number of segregating sites, D, thetaH, H).
  4. print out following info:
     
     • model (0 for model-A, 1 for model-B)
       
     • T1 (drawn from the prior distribution)
       
     • a vector of summary statistics
  5. Goto Step 1 and repeat this many many times.

  The final output will look like:

  M	T1	pi	ss	D	thetaH	H
  0	0.605	2.07	6	-0.48	1.35	0.71
  0	0.034	18.75	40	1.15	17.53	1.21
  1	0.045	8.32	24	-0.53	5.67	2.64
  :	:
  :	:
  

• Simulation code (msDrive.pl)

• Migration matrix
  • -ma x 8.0 16.0 x is equivalent to:
    
    $$\left( \begin{array}{c c} 4 N_0 m_{11} & 4 N_0 m_{12} 4 N_0... ...\left( \begin{array}{c c} x & 8.0 16.0 & x \end{array}\right)$$
    
    
    $m_{ij}$ is the fraction of subpopulation i which was in the j-th subpopulation in the previous generation.
    • A little tricky because two different subpop sizes.
    • We specified $\theta = 4 N_0 \mu = 8$ (-s 8). This means that the ancestral population size = 1, $N_1$ = 1/2, and $N_2$ = 1/4.
    • So actual numbers of migrants into the subpop 1 & 2 are $2 N_1 m_{12} = 8.0 / 2 / 2 = 2.0$ and $2 N_2 m_{21} = 16.0 / 4 / 2 = 2.0$
    • So this is a symmetric migration model in terms of actual number of migrants.
    • more details in p. 12 of ms manual .

Background Info/Explanation of msDrive.pl

• For loop:

  #!/usr/bin/perl -w
  for ($i = 0; $i < 10; $i++) {
    $num = rand(5);
    print "repetion $i: $num\n";
  }
  

  In this loop, $i goes from 0 to 9.

Exercise: modify msDrive.pl

Instead of setting the migration rate to a pre-fixed value, modify the program to estimate the migration rate (in addition to divergence time). In other words, you need to draw the migration from a random distribution.