library(RColorBrewer) L <- read.table("~/L14/wfabc/sd3/incomp/sel_L14incomp.txt",header=FALSE,sep=" ") ##for sd =.3 ## # significant and magnitutde # P to F4 mean(abs(L[(L[,4] > 0 | L[,5] <0),1])) #.7087367 mean(abs(L[,1])) # 3543269 length(L[(L[,4] > 0 | L[,5] <0),1]) #49 # F4 to F16A mean(abs(L[(L[,9] > 0 | L[,10] <0),6])) #0.4352378 mean(abs(L[,6])) # 0.1967997 length(L[(L[,9] > 0 | L[,10] <0),6]) #49 # F4 to F16B mean(abs(L[(L[,14] > 0 | L[,15] <0),11])) #0.4248731 mean(abs(L[,11])) # 0.2044151 length(L[(L[,14] > 0 | L[,15] <0),11]) #56 ## how many stay the same/flip for A same <- which(L[,4]>0 & L[,9]>0 | L[,5]<0 & L[,10]<0) flip <- which(L[,4]>0 & L[,10]<0 | L[,5]<0 & L[,9]>0) ## how many stay the same/flip for B same <- which(L[,4]>0 & L[,9]>0 | L[,5]<0 & L[,10]<0) flip <- which(L[,4]>0 & L[,10]<0 | L[,5]<0 & L[,9]>0) ## how many stay the same/flip for A B same <- which((L[,9]>0 & L[,14]>0) | (L[,10]<0 & L[,15]<0)) flip <- which(L[,4]>0 & L[,15]<0 | L[,5]<0 & L[,14]>0) (L[,4]>0 & L[,9]>0) | (L[,5]<0 & L[,10]<0) ## Summary of H's H <- read.table("~/L14/wfabc/sd3/incomp/het_L14incomp.txt",header=FALSE,sep=" ") ##for sd =.3 L <- read.table("~/L14/wfabc/sd3/incomp/sel_L14incomp.txt",header=FALSE,sep=" ") ##for sd =.3 ## Individual selection plots pdf('scatter_sh.pdf', width=7, height=21) par(mar=c(2.5,2,2.5,2), mfrow=c(3,1),oma=c(5,5.5,2,0)) plot(L[,1],H[,1], xlab='', ylab='', cex=1.75, cex.axis=2, ylim=c(0,1), xlim=c(-1,1)) abline(v=0,h=.5,lty=2) legend('bottomright',legend=(paste('cor =', round(cor(L[,1], H[,1]),3))), bty='n', cex=1.5) title(main="(a) P to F4",adj=0,cex.main=2, line=.4) plot(L[,6],H[,6], xlab='', ylab='', cex=1.75, cex.axis=2, ylim=c(0,1), xlim=c(-1,1)) abline(v=0,h=.5, lty=2) legend('bottomright',legend=(paste('cor =', round(cor(L[,6], H[,6]),3))), bty='n', cex=1.5) title(main="(b) F4 to F16A",adj=0,cex.main=2, line=.4) plot(L[,11],H[,11], xlab='',ylab='',cex=2, cex.axis=2, ylim=c(0,1), xlim=c(-1,1)) abline(v=0,h=.5,lty=2) legend('bottomright',legend=(paste('cor =', round(cor(L[,11], H[,11]),3))), bty='n',cex=1.5) title(main="(c) F4 to F16B",adj=0,cex.main=2, line=.4) mtext(text='Selection Coefficients', side=1, line=2, outer=TRUE, cex=3) mtext(text='Heterozygous Effect', side=2, line=2, outer=TRUE, cex=3) dev.off() S <- read.table("~/L14/wfabc/sd3/incomp/pp_incomp.txt",header=FALSE,sep=" ") ## which selection model 'won' for each SNP, posterior probabilities of each model model <- apply(S,1,function(x) which(x==max(x))) length(which(model==1)) ## 53 length(which(model==2)) ## 93 length(which(model==3)) ## 105 pdf('pp_box.pdf') par(mar=c(5,5,2,2)) boxplot(S, xlab='scenario', ylab='posterior probability', col='grey', xaxt='n', cex.lab=2, cex.main=2) axis(1, at=1:3, labels=c(1,2,3)) dev.off()