Libraries

library(tidyr)
library(magrittr)
## 
## Attaching package: 'magrittr'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
library(dplyr)
## 
## Attaching package: 'dplyr'
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(data.table)
## 
## Attaching package: 'data.table'
## 
## The following objects are masked from 'package:dplyr':
## 
##     between, last

Colors

cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
  1. Run sims in fwdpp (bneck.sh)
  2. Split neutral and selected sites into files (split.sh) (??)
  3. Concatenate neutral sites
  4. Run ms_slide
  5. profit!

Functions

#get file, do cleaning/filtering/selecting etc.
clean<-function(filename){
  fread(paste("~/Projects/b_project/bneck/results/windows/",filename,sep="")) %>% 
  setnames(c("window","pi","singletons")) %>%
   filter(window>39,window<99) 
}

#compare two files, plot the smoothed splines
compare_windows<-function(file1,file2,main,type1,type2){
  simdata1<-clean(file1)
  simdata2<-clean(file2)

#pi
stat="pi"
spline1=smooth.spline(simdata1[,get(stat)]~simdata1$window)
spline2=smooth.spline(simdata2[,get(stat)]~simdata2$window)
lim=min(spline1$y/max(spline1$y),spline2$y/max(spline2$y))
plot(spline1$y/max(spline1$y)~spline1$x,ylim=c(lim-0.05,1),type="l",col=cbbPalette[2],lwd=3,main=main,ylab=stat,xaxt="n",xlab="bp")
axis(side=1,at=4:10*10,labels=(4:10*1000))
lines(spline2$y/max(spline2$y)~spline2$x,col=cbbPalette[3],lwd=3)
legend("bottomright",box.lwd=0,lty=1,col=c(cbbPalette[2],cbbPalette[3]),legend=c(type1,type2),lwd=3)

#singletons
stat="singletons"
spline1=smooth.spline(simdata1[,get(stat)]~simdata1$window)
spline2=smooth.spline(simdata2[,get(stat)]~simdata2$window)
lim=min(spline1$y/max(spline1$y),spline2$y/max(spline2$y))
plot(spline1$y/max(spline1$y)~spline1$x,ylim=c(lim-0.05,1),type="l",col=cbbPalette[2],lwd=3,main=main,ylab=stat,xaxt="n",xlab="bp")
axis(side=1,at=4:10*10,labels=(4:10*1000))
lines(spline2$y/max(spline2$y)~spline2$x,col=cbbPalette[3],lwd=3)
legend("bottomright",box.lwd=0,lty=1,col=c(cbbPalette[2],cbbPalette[3]),legend=c(type1,type2),lwd=3)
}

h=0.01 lots o’ rescessive

s=0.0001 in maize and teo

compare_windows("windows.neutral_sims.3921128","windows.neutral_sims.3880982","taxa for s=0.0001","maize","teosinte")

ain’t no difference

s=0.001 in maize and teo

compare_windows("windows.neutral_sims.3906126","windows.neutral_sims.3926161","taxa for s=0.001","maize","teosinte")

No effect or difference on singletons, but stronger effect on maize than teosinte!

s=0.01 in maize and teo

compare_windows("windows.neutral_sims.3911128","windows.neutral_sims.3896080","taxa for s=0.01","maize","teosinte")

No effect and no difference on singletons, big effect but virtually the same (very slightly stronger in maize) for pi.

h=0.1 kinda rescessive

s=0.0001 in maize and teo

compare_windows("windows.neutral_sims.3946720","windows.neutral_sims.3951721","taxa for h=0.1 s=0.0001","maize","teosinte")

s=0.001 in maize and teo

compare_windows("windows.neutral_sims.3956721","windows.neutral_sims.3941720","taxa for h=0.1 s=0.001","maize","teosinte")