bnecks

Dumb stuff:

options(scipen=999) # no printy scientific notation in our sims damnit
library(ggplot2)

Let’s do a bneck:

sims=1000000
bneck_end=runif(sims,0,2) # end comes first backward in time
bneck_start=bneck_end+0.005 # start comes 2nd so is bigger number
tbsfile<-cbind(bneck_end,bneck_start)
write.table(file="tbsfile",tbsfile,col.names=F,row.names=F,quote=F)
# we do with recombination, with theta=10 or N~80,000
#bneck then is 0.005*4*80,000 = 1,600 generations
f=pipe(paste("~/src/msdir/ms 10 ", sims, " -t 10 -r 10 1000 -eN tbs 0.1 -eN tbs 1 < tbsfile | ~/src/msdir/sample_stats | cut -f 6",sep=""))
taj_d=scan(f)

Let’s plot that shit:

qplot(bneck_end,taj_d,geom="smooth")

## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.