Post date: Sep 04, 2015 9:33:19 PM
I generated new plots for the null expectations for allele frequency differences (change) between M and L lines and L and LR lines. In doing so I caught what I am pretty sure was an error in the analysisNullSimP function, fst should be 1-(1-1/(2Ne))^t not 1-(1-2/(2Ne))^t. I corrected this. The code also generated files for the ABC analysis, that is the exceptional SNP in the lentil lines pLentilSnps.txt and for the lentil or reversion lines pLReversionSnps.txt. I added the estimates of Ne as a header for this files, which are currently in/home/zgompert/labs/evolution/projects/cmaclentil/popgen.
Updated code in /home/zgompert/labs/evolution/projects/cmaclentil/popgen/SeedBeetleSim/sims2.R:
## for p-values and first test of trade-offs
analyNullSimP<-function(p0=0.5,p1=0.5,t=100,ne=100,lower=TRUE){
fst<-1-(1-1/(2*ne))^t
alpha<-p0 * (1-fst)/fst
beta<-(1-p0) * (1-fst)/fst
pv<-pbeta(p1,alpha+0.001,beta+0.001,lower.tail=lower)
pv
}
bnd<-log(0.99995/0.00005) ## equivalent to p < 0.0001
## Ne values for branches separating populations
## M v L1 v L1R, logit quantile
pL1<-read.table("p_L1-F100.txt",header=FALSE,sep=",")
pL1R<-read.table("p_L1R-F46.txt",header=FALSE,sep=",")
p0<-pL1[,1]
p1<-pL1R[,1]
pv<-analyNullSimP(p0,p1,t=84,ne=455.5,lower=TRUE)
l1rpv<-log(pv/(1-pv))
pM<-read.table("p_M-14.txt",header=FALSE,sep=",")
p0<-pM[,1]
p1<-pL1[,1]
mpv<-analyNullSimP(p0,p1,t=231,ne=454.5,lower=TRUE)
l1pv<-log(mpv/(1-mpv))
## M v L2 v L2R, logit quantile
pL2<-read.table("p_L2-F87.txt",header=FALSE,sep=",")
pL2R<-read.table("p_L2R-F45.txt",header=FALSE,sep=",")
p0<-pL2[,1]
p1<-pL2R[,1]
pv<-analyNullSimP(p0,p1,t=84,ne=437.5,lower=TRUE)
l2rpv<-log(pv/(1-pv))
pM<-read.table("p_M-14.txt",header=FALSE,sep=",")
p0<-pM[,1]
p1<-pL2[,1]
mpv<-analyNullSimP(p0,p1,t=194,ne=525,lower=TRUE)
l2pv<-log(mpv/(1-mpv))
## M v L3 v L3R, logit quantile
pL3<-read.table("p_L3-F85.txt",header=FALSE,sep=",")
pL3R<-read.table("p_L3R-F45.txt",header=FALSE,sep=",")
p0<-pL3[,1]
p1<-pL3R[,1]
pv<-analyNullSimP(p0,p1,t=85,ne=937,lower=TRUE)
l3pv<-log(pv/(1-pv))
pM<-read.table("p_M-14.txt",header=FALSE,sep=",")
p0<-pM[,1]
p1<-pL3[,1]
mpv<-analyNullSimP(p0,p1,t=189,ne=549.5,lower=TRUE)
l3rpv<-log(mpv/(1-mpv))
## manhattan plots
pdf("manhatL.pdf",width=9,height=9)
par(mfrow=c(3,1))
par(mar=c(4,5.5,3,1.5))
plot(abs(l1pv),pch=19,col="darkgray",cex=0.5,ylab="logit quantile",xlab="",cex.lab=1.7,cex.axis=1.3)
abline(h=bnd,lty=3,col="red",lwd=1.4)
title(main="(a) L1",adj=0,cex.main=1.6)
plot(abs(l2pv),pch=19,col="darkgray",cex=0.5,ylab="logit quantile",xlab="",cex.lab=1.7,cex.axis=1.3)
abline(h=bnd,lty=3,col="red",lwd=1.4)
title(main="(b) L2",adj=0,cex.main=1.6)
plot(abs(l3pv),pch=19,col="darkgray",cex=0.5,ylab="logit quantile",xlab="locus (SNP)",cex.lab=1.7,cex.axis=1.3)
abline(h=bnd,lty=3,col="red",lwd=1.4)
title(main="(c) L3",adj=0,cex.main=1.6)
dev.off()
pdf("manhatRev.pdf",width=9,height=9)
par(mfrow=c(3,1))
par(mar=c(4,5.5,3,1.5))
plot(abs(l1rpv),pch=19,col="darkgray",cex=0.5,ylab="logit quantile",xlab="",cex.lab=1.7,cex.axis=1.3)
abline(h=bnd,lty=3,col="red",lwd=1.4)
title(main="(a) L1R",adj=0,cex.main=1.6)
plot(abs(l2rpv),pch=19,col="darkgray",cex=0.5,ylab="logit quantile",xlab="",cex.lab=1.7,cex.axis=1.3)
abline(h=bnd,lty=3,col="red",lwd=1.4)
title(main="(b) L2R",adj=0,cex.main=1.6)
plot(abs(l3rpv),pch=19,col="darkgray",cex=0.5,ylab="logit quantile",xlab="locus (SNP)",cex.lab=1.7,cex.axis=1.3)
abline(h=bnd,lty=3,col="red",lwd=1.4)
title(main="(c) L3R",adj=0,cex.main=1.6)
dev.off()
## tradeoffs
pdf("tradeoffsNullMod.pdf",width=6,height=18)
par(mfrow=c(3,1))
par(mar=c(4,5.5,3,1.5))
plot(l1pv,l1rpv,pch=19,col="darkgray",cex=0.5,xlab="logit quantile (M->L1)",ylab="logit quantile (L1->L1R)",cex.lab=1.7,cex.axis=1.3)
abline(h=c(bnd,-bnd),lty=3,col="red",lwd=1.4)
abline(v=c(bnd,-bnd),lty=3,col="red",lwd=1.4)
abline(h=0,lty=2)
abline(v=0,lty=2)
title(main="(a) L1 vs. L1R",adj=0,cex.main=1.6)
plot(l2pv,l2rpv,pch=19,col="darkgray",cex=0.5,xlab="logit quantile (M->L2)",ylab="logit quantile (L2->L2R)",cex.lab=1.7,cex.axis=1.3)
abline(h=c(bnd,-bnd),lty=3,col="red",lwd=1.4)
abline(v=c(bnd,-bnd),lty=3,col="red",lwd=1.4)
abline(h=0,lty=2)
abline(v=0,lty=2)
title(main="(b) L2 vs. L2R",adj=0,cex.main=1.6)
plot(l3pv,l3rpv,pch=19,col="darkgray",cex=0.5,xlab="logit quantile (M->L3)",ylab="logit quantile (L3->L3R)",cex.lab=1.7,cex.axis=1.3)
abline(h=c(bnd,-bnd),lty=3,col="red",lwd=1.4)
abline(v=c(bnd,-bnd),lty=3,col="red",lwd=1.4)
abline(h=0,lty=2)
abline(v=0,lty=2)
title(main="(c) L3 vs. L3R",adj=0,cex.main=1.6)
dev.off()
## counts and parallelism
sum(abs(l1pv) > bnd) # 314
sum(abs(l2pv) > bnd) # 231
sum(abs(l3pv) > bnd) # 188
sum((l1pv > bnd & l2pv > bnd) | (l1pv < -bnd & l2pv < -bnd)) # 56; = 17.8% and 24.2% shared
sum((l1pv > bnd & l3pv > bnd) | (l1pv < -bnd & l3pv < -bnd)) # 0
sum((l3pv > bnd & l2pv > bnd) | (l3pv < -bnd & l2pv < -bnd)) # 0
sum((l1pv > bnd & l2pv > bnd & l3pv > bnd) | (l1pv < -bnd & l2pv < -bnd & l3pv < -bnd)) # 0
sum(abs(l1rpv) > bnd) # 182
sum(abs(l2rpv) > bnd) # 192
sum(abs(l3rpv) > bnd) # 214
sum((l1rpv > bnd & l2rpv > bnd) | (l1rpv < -bnd & l2rpv < -bnd)) # 0
sum((l1rpv > bnd & l3rpv > bnd) | (l1rpv < -bnd & l3rpv < -bnd)) # 4
sum((l3rpv > bnd & l2rpv > bnd) | (l3rpv < -bnd & l2rpv < -bnd)) # 3
sum((l1rpv > bnd & l2rpv > bnd & l3rpv > bnd) | (l1rpv < -bnd & l2rpv < -bnd & l3rpv < -bnd)) # 0
## write combined allele frequencies for ABC for (i) candidate lentil snps (ii) candidate lentil and reversion snps
lsnps<-unique(which(abs(l1pv) > bnd | abs(l2pv) > bnd | abs(l3pv) > bnd))#673 snps
rsnps<-unique(which(abs(l1rpv) > bnd | abs(l2rpv) > bnd | abs(l3rpv) > bnd))#569 snps
## 215 SNPs in both data sets!
lrsnps<-unique(c(lsnps,rsnps)) #1024 snps
lids<-read.table("../locusids.txt",header=FALSE)
plsnps<-data.frame(lids[lsnps,1],pM[lsnps,1],pL1[lsnps,1],pL2[lsnps,1],pL3[lsnps,1])
plrsnps<-data.frame(lids[lrsnps,1],pM[lrsnps,1],pL1[lrsnps,1],pL2[lrsnps,1],pL3[lrsnps,1],pL1R[lrsnps,1],pL2R[lrsnps,1],pL3R[lrsnps,1])
write.table(plsnps,"../pLentilSnps.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
write.table(plrsnps,"../pLReversionSnps.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)
save(list=ls(),file="null.Rdata")