#16 Populations diverge more than expected

[Anonymous]

I think migration rate is lower than expected from the setting parameters. Consider the following bash code that run 500 iterations of a two island model, calculate the one weir and Cockerham Fst value for each simulation and plot the histogram of Fst (via R) as well as the expected Fst value from Slatkin 1991 $\left(\frac{1}{1+N (m+\mu) \left(\frac{d}{d-1}\right)^2}\right)$.

Note that I used the last version (sfscode_20150910) for these simulations.

# Please set up the path to the SFScode executable and set a working directory
sfs="/Path/To/SFSexecutable"
workingDirectory="~"

# Parameters
nbiterations=500
N=1000
n=500
m=0.05
M=$(bc <<< "scale=2; 2*${N}*${m}")
mu=0.000001
theta=$(bc <<< "scale=2; 2*${N}*${mu}")
L=999

# Get to working directory
cd ${workingDirectory}

# Create quick R code to make the histogram
printf "SlatkinExpectedFst=1/(1+4*${N}*(${m}+${L}*${mu})*4)\nd=read.table(\"${PWD}/1.windowed.weir.fst\",header=TRUE)\nfor (i in 2:${nbiterations}) \n{\nd=rbind(d,read.table(paste0(\"${PWD}/\",i,\".windowed.weir.fst\"),header=TRUE))\n}\n\npdf(\"${PWD}/hist.pdf\")\nhist(d\$WEIGHTED_FST,breaks=50,xlim=range(c(SlatkinExpectedFst,d\$WEIGHTED_FST)))\nabline(v=SlatkinExpectedFst, col='red')\ndev.off()\n\nprint(paste0(\"mean fst in simulations = \",mean(d\$WEIGHTED_FST),\"  |  expected mean fst = \",SlatkinExpectedFst))" > plot_histogram.r

# Build individuals-population relationship files for vcftools
eval echo N0.{0..$(( n - 1 ))} | tr " " "\n" > 0_pop.txt
eval echo N1.{0..$(( n - 1 ))} | tr " " "\n" > 1_pop.txt

# Run simulations and calculate Weir and Cockerham Fst
for ((i=1; i<=nbiterations; i++));do
      ${sfs} 2 1 -TS 0 0 1 -TE 5 --VCF -o ${i}.vcf -N ${N} -n ${n} -L 1 ${L} -m A ${M} -t ${theta}
      vcftools --vcf ${i}.vcf --weir-fst-pop 0_pop.txt --weir-fst-pop 1_pop.txt --fst-window-size 9999 --fst-window-step 999 --out ${i} 
done

# Remove individuals-population relationship files
rm -f 0_pop.txt
rm -f 0_pop.txt

# Build histogram
Rscript plot_histogram.r
open ${PWD}/hist.pdf

Here is the resulting histogram

In red is the expected Fst value (~0.0012). As we can see, the fst values from the simulations are very high (with a mean at ~0.63) compared to the one predicted by classic models. Is there something I misunderstand or is there a bug somewhere?

Thank you,
Remi

[Anonymous]

The histogram failed to display apparently! I am adding the histogram in attachement. You can also find it here.

[Anonymous]

I'm trying to parse your code here... Would be helpful to include the output of this printf statement rather than trying to read through it. Are you calculating the expected value right? You have:

1/(1+4${N}(${m}+${L}${mu})4)

That doesn't seem quite right to me, it seems to have a couple extra 4s...? Also, remind me what d is in the equation from Slatkin 91.

Cheers,

Ryan

[Anonymous]

[Anonymous]

[Anonymous]

[Anonymous]