Untitled

RUN1

Current runs from farm:

s	run number
0.0001	1026262
0.00001	1026362
0.000001	1026462
0.01	1026562
0.001	1026157

Sims:

~/src/fwdpp/examples/bneck_selection_dist $N $neutral_theta $selected_theta $rho $s $pp $h $g1 $N2 $N3 $g2 $n $nreps $split $see

seed=$(( $SLURM_ARRAY_TASK_ID*$RANDOM ))
N=15000
neutral_theta=12
selected_theta=24
rho=$(($neutral_theta + $selected_theta ))
s=$1 #selection <- CHECK WHICH VERSION OF CODE YOU ARE USING. s is either fixed s or mean of exponential!
pp=0.001 #proportion of selected mutations that are positively selected
h=0.01 #dominance
g1=150000 #generations between N and N2?
N2=12000 #bneck size
N3=1500000 #new NE at end
n=23 #number samples
g2=90
nreps=1
split=1 #0 combines neutral and selected SNPs in one ms block; 1 splits them into two blocks. changing this changes the downstream

seldata=data.frame();
header<-read.table(colClasses = "character","~/Desktop/header")
sel<-10^(-(2:6))
nums<-c(1026562,1026157,1026262,1026362,1026462)
for(i in 1:5){
temp<-read.table(colClasses = "numeric",paste("~/Desktop/neutral.",nums[i],sep=""))
temp=cbind(rep(sel[i],length(temp$V1)),temp)
seldata=rbind(seldata,temp)
}

colnames(seldata)=c("selection",header)
ns<-seldata$selection*15000
boxplot(seldata$pi~ns,ylab=expression(pi),xlab="Ns",main="Neutral mutations"); abline(h=12)

#boxplot(seldata$tajd~ns,xlab="Ns",ylab="Taj D",main="Neutral mutations")

seldata=data.frame();
header<-read.table(colClasses = "character","~/Desktop/header")
sel<-10^(-(2:6))
nums<-c(1026562,1026157,1026262,1026362,1026462)
for(i in 1:5){
temp<-read.table(colClasses = "numeric",paste("~/Desktop/selected.",nums[i],sep=""))
temp=cbind(rep(sel[i],length(temp$V1)),temp)
seldata=rbind(seldata,temp)
}

colnames(seldata)=c("selection",header)
ns<-seldata$selection*15000
boxplot(seldata$pi~ns,ylab=expression(pi),xlab="Ns",main="Selected mutations"); abline(h=24)

RUN 2

Smaller bneck size

s	run number
0.0001	1042154
0.00001	1042254
0.000001	1042354
0.01	1041954
0.001	1042054

Sims:

~/src/fwdpp/examples/bneck_selection_dist $N $neutral_theta $selected_theta $rho $s $pp $h $g1 $N2 $N3 $g2 $n $nreps $split $see

seed=$(( $SLURM_ARRAY_TASK_ID*$RANDOM ))
N=15000
neutral_theta=12
selected_theta=24
rho=$(($neutral_theta + $selected_theta ))
s=$1 #selection <- CHECK WHICH VERSION OF CODE YOU ARE USING. s is either fixed s or mean of exponential!
pp=0.001 #proportion of selected mutations that are positively selected
h=0.01 #dominance
g1=150000 #generations between N and N2?
N2=1500 #bneck size
N3=1500000 #new NE at end
n=23 #number samples
g2=90
nreps=1

seldata=data.frame();
header<-read.table(colClasses = "character","~/Desktop/header")
sel<-10^(-(2:6))
nums<-c(1041954,1042054, 1042154,1042254,1042354)
for(i in 1:5){
temp<-read.table(colClasses = "numeric",paste("~/Desktop/neutral.",nums[i],sep=""))
temp=cbind(rep(sel[i],length(temp$V1)),temp)
seldata=rbind(seldata,temp)
}

colnames(seldata)=c("selection",header)
ns<-seldata$selection*15000
boxplot(seldata$pi~ns,ylab=expression(pi),xlab="Ns",main="Neutral mutations"); abline(h=12)

#boxplot(seldata$pi~ns,ylab=expression(pi),xlab="Ns",main="Neutral mutations",add=T); abline(h=12)

boxplot(seldata$tajd~ns,xlab="Ns",ylab="Taj D",main="Neutral mutations")

rho=rep(36,length(seldata$pi))
lowrho=cbind(seldata,rho)

seldata=data.frame();
header<-read.table(colClasses = "character","~/Desktop/header")
sel<-10^(-(2:6))
nums<-c(1041954,1042054, 1042154,1042254,1042354)
for(i in 1:5){
temp<-read.table(colClasses = "numeric",paste("~/Desktop/selected.",nums[i],sep=""))
temp=cbind(rep(sel[i],length(temp$V1)),temp)
seldata=rbind(seldata,temp)
}

colnames(seldata)=c("selection",header)
ns<-seldata$selection*15000
boxplot(seldata$pi~ns,ylab=expression(pi),xlab="Ns",main="Selected mutations"); abline(h=24)

RUN 3

Smaller bneck size, \(4\times \rho\)

s	run number
0.0001	1042468
0.00001	1042768
0.000001	1042668
0.001	1042568

Sims:

~/src/fwdpp/examples/bneck_selection_dist $N $neutral_theta $selected_theta $rho $s $pp $h $g1 $N2 $N3 $g2 $n $nreps $split $see seed=$(( $SLURM_ARRAY_TASK_ID*$RANDOM ))
N=15000
neutral_theta=12
selected_theta=24
rho=$(( 4*($neutral_theta + $selected_theta) ))
s=$1 #selection <- CHECK WHICH VERSION OF CODE YOU ARE USING. s is either fixed s or mean of exponential!
pp=0.001 #proportion of selected mutations that are positively selected
h=0.01 #dominance
g1=150000 #generations between N and N2?
N2=1500 #bneck size
N3=1500000 #new NE at end
n=23 #number samples
g2=90
nreps=1

seldata=data.frame();
header<-read.table(colClasses = "character","~/Desktop/header")
sel<-10^(-(3:6))
nums<-c(1042568,1042468, 1042768,1042668)
for(i in 1:4){
temp<-read.table(colClasses = "numeric",paste("~/Desktop/neutral.",nums[i],sep=""))
temp=cbind(rep(sel[i],length(temp$V1)),temp)
seldata=rbind(seldata,temp)
}

colnames(seldata)=c("selection",header)
ns<-seldata$selection*15000
boxplot(seldata$pi~ns,ylab=expression(pi),xlab="Ns",main="Neutral mutations"); abline(h=12)

#boxplot(seldata$pi~ns,ylab=expression(pi),xlab="Ns",main="Neutral mutations",add=T,col="red",lwd=2); abline(h=12)

boxplot(seldata$tajd~ns,xlab="Ns",ylab="Taj D",main="Neutral mutations")

rho=rep(144,length(seldata$pi))
seldata=cbind(seldata,rho)
rhodata=rbind(seldata,lowrho)

seldata=data.frame();
header<-read.table(colClasses = "character","~/Desktop/header")
sel<-10^(-(2:6))
nums<-c(1041954,1042054, 1042154,1042254,1042354)
for(i in 1:5){
temp<-read.table(colClasses = "numeric",paste("~/Desktop/selected.",nums[i],sep=""))
temp=cbind(rep(sel[i],length(temp$V1)),temp)
seldata=rbind(seldata,temp)
}

colnames(seldata)=c("selection",header)
ns<-seldata$selection*15000
boxplot(seldata$pi~ns,ylab=expression(pi),xlab="Ns",main="Selected mutations"); abline(h=24)

Comparison

library(ggplot2)
ggplot(data=rhodata,aes(factor(selection),pi))+geom_hline(aes(yintercept=12))+geom_boxplot(data=rhodata,aes(fill=factor(rhodata$rho)))

ggplot(data=rhodata,aes(factor(selection),tajd))+geom_hline(aes(yintercept=0))+geom_boxplot(data=rhodata,aes(fill=factor(rhodata$rho)))