Post date: Sep 18, 2015 8:56:58 PM
All of the results so far were based on the reversion lines samples at F45 or F46. In order to further test for AP and to ask whether loss of adatpation to lentil was ongoing in the reversion lines, I grabbed allele frequencies for all of the AP SNPs in each line from the LX LXR 35 and 45(6) samples.
perl makeRevFiles.pl l1APsnps.txt ../locusids.txt ../SeedBeetleSim/p_L1R-F46.txt ../SeedBeetleSim/p_L1R-F35.txt ../SeedBeetleSim/p_L1-F100.txt ../SeedBeetleSim/p_M-14.txt
perl makeRevFiles.pl l2APsnps.txt ../locusids.txt ../SeedBeetleSim/p_L2R-F45.txt ../SeedBeetleSim/p_L2R-F35.txt ../SeedBeetleSim/p_L2-F87.txt ../SeedBeetleSim/p_M-14.txt
perl makeRevFiles.pl l3APsnps.txt ../locusids.txt ../SeedBeetleSim/p_L3R-F45.txt ../SeedBeetleSim/p_L3R-F35.txt ../SeedBeetleSim/p_L3-F85.txt ../SeedBeetleSim/p_M-14.txt
I then examined whether the change between generations 35 and 45(6) could be explained by a null model of drift (at p = 0.05) and examined correlation in change between these generation and between the lentil lines and generation 45(6) samples from the reversion lines. All of the relevant files are in /labs/evolution/projects/cmaclentil/popgen/reversionSims/.
R code in rev.R:
## evidence of recent selection in reversion lines
## note columns are SonL SonM pLR4X pLR3X pLX pM14
l1<-read.table("revApDatal1.txt",header=FALSE)
l2<-read.table("revApDatal2.txt",header=FALSE)
l3<-read.table("revApDatal3.txt",header=FALSE)
#ne L1R L2R L3R
Ne<-c(533.0,512.9,762.3)
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
}
## null L1R
p0<-l1[,4]
p1<-l1[,3]
pv<-analyNullSimP(p0,p1,t=11,ne=Ne[1],lower=TRUE)
sum(pv < 0.025 | pv > 0.975)
# 21 sig
sigl1<-which(pv < 0.025 | pv > 0.975)
## null L2R
p0<-l2[,4]
p1<-l2[,3]
pv<-analyNullSimP(p0,p1,t=10,ne=Ne[2],lower=TRUE)
sum(pv < 0.025 | pv > 0.975)
# 32 sig
sigl2<-which(pv < 0.025 | pv > 0.975)
## null L3R
p0<-l3[,4]
p1<-l3[,3]
pv<-analyNullSimP(p0,p1,t=10,ne=Ne[3],lower=TRUE)
sum(pv < 0.025 | pv > 0.975)
# 25 sig
sigl3<-which(pv < 0.025 | pv > 0.975)
pdf("reversion.pdf",width=4,height=12)
par(mfrow=c(3,1))
par(mar=c(4,5.5,2,0.5))
cl<-1.6
ca<-1.3
cm<-1.3
cs<-c("gray","#1B9E77")
sym<-c(20,21)
## L1
## L1R-F46 - L1-F100
x<-l1[,3]-l1[,5]
## L1R-F46 - L1R-F35
y<-l1[,3]-l1[,4]
a<-which((x*y) > 0)
indx<-rep(1,length(x))
indx[a]<-2
indx2<-rep(1,length(x))
indx2[sigl1]<-2
plot(x,y,cex.lab=cl,cex.axis=ca,xlab=expression(paste(Delta[p]," L1R-F46 - L1-F100")),ylab=expression(paste(Delta[p]," L1R-F46 - L1R-F35")),pch=sym[indx2],col=cs[indx])
abline(h=0,lty=2)
abline(v=0,lty=2)
out<-cor.test(x,y)
r<-round(out$estimate,3)
title(main=paste("(a) L1R-F46 vs. L1R-F35, r = ",r," p < 0.001"),adj=0,cex.main=cm)
legend(-0.8,0.27,c("same, p < 0.05","same, p > 0.05","diff. p < 0.05","diff. p > 0.05"),pch=rep(rev(sym),2),col=rep(rev(cs),each=2))
## L2
## L2R-F45 - L2-F87
x<-l2[,3]-l2[,5]
## L2R-F45 - L2R-F35
y<-l2[,3]-l2[,4]
a<-which((x*y) > 0)
indx<-rep(1,length(x))
indx[a]<-2
indx2<-rep(1,length(x))
indx2[sigl2]<-2
plot(x,y,cex.lab=cl,cex.axis=ca,xlab=expression(paste(Delta[p]," L2R-F45 - L2-F87")),ylab=expression(paste(Delta[p]," L2R-F45 - L2R-F35")),pch=sym[indx2],col=cs[indx])
abline(h=0,lty=2)
abline(v=0,lty=2)
out<-cor.test(x,y)
r<-round(out$estimate,3)
title(main=paste("(b) L2R-F45 vs. L2R-F35, r = ",r," p < 0.001"),adj=0,cex.main=cm)
## L3
## L3R-F45 - L3-F8X
x<-l3[,3]-l3[,5]
## L3R-F45 - L3R-F35
y<-l3[,3]-l3[,4]
a<-which((x*y) > 0)
indx<-rep(1,length(x))
indx[a]<-2
indx2<-rep(1,length(x))
indx2[sigl3]<-2
plot(x,y,cex.lab=cl,cex.axis=ca,xlab=expression(paste(Delta[p]," L3R-F45 - L3-F85")),ylab=expression(paste(Delta[p]," L3R-F45 - L3R-F35")),pch=sym[indx2],col=cs[indx])
abline(h=0,lty=2)
abline(v=0,lty=2)
out<-cor.test(x,y)
r<-round(out$estimate,3)
title(main=paste("(c) L3R-F45 vs. L3R-F35, r = ",r," p < 0.001"),adj=0,cex.main=cm)
dev.off()