bad R/ms simulation of bottlenecks

bad_bottleneck.R

library(tidyverse)
library(cowplot)
Ne=25000
mu=1E-8
win_size=20000
theta=win_size*4*mu*Ne
before<-as.numeric(sapply(1:20, function(i) 
  system(paste("~/src/msdir/ms 10 1000 -t", 4*theta, "-r", 4*theta, "1000 | ~/src/msdir/sample_stats | cut -f 2 | ~/src/msdir/stats | cut -f 1")
         ,intern=T)))

div<-as.numeric(sapply(1:100, function(i) 
  system(paste("~/src/msdir/ms 10 1000 -t", theta, "-r", theta, "1000 -eN", i/1000, "4 | ~/src/msdir/sample_stats | cut -f 2 | ~/src/msdir/stats | cut -f 1")
         ,intern=T)))
div=c(before,div)
bob=data.frame(c(-1*20:1*100,4*Ne*1:100/1000),div,rep(2,120))
colnames(bob)=c("time","div","epochs")

ggplot(data=bob,aes(y=div,x=time))+
  geom_point(size=2)+
  ylab("nucleotide diversity")+
  xlab("gens since bottleneck")