Post date: Mar 03, 2017 8:48:7 PM
I recalculated genomic change (and a null distribution) for the OGA (science experiment) data with the defined loci mel and stripy (see notes here). Right now the results and code are in gompert:/home/zgompert/Local/timemaGenomicChange as kingspeak is down for maintenance.
For mel, I focused on total change across all bushes, and for stripy on divergence between A and C. Across all bushes (like for the Ecology Letters bugs). The signal for mel is whopping (in contrast to the Ecology Letters). Melanistic bugs were common and this may in part reflect dispersal, but the signal is very strong (p < 1/40 = 0.025). stripy has the greatest change (divergence) and change for diversity between A and C. The results are quantitatively less strong, but still clear (p < 1/21 = 0.048).
Here is the R code:
## read in allele freq. and id data
p0<-scan("pEM_OGA2010.txt")
p1a<-scan("pEM_MR1_A.txt")
p2a<-scan("pEM_MR2_A.txt")
p3a<-scan("pEM_MR3_A.txt")
p4a<-scan("pEM_MR4_A.txt")
p5a<-scan("pEM_MR5_A.txt")
p1c<-scan("pEM_MR1_C.txt")
p2c<-scan("pEM_MR2_C.txt")
p3c<-scan("pEM_MR3_C.txt")
p4c<-scan("pEM_MR4_C.txt")
p5c<-scan("pEM_MR5_C.txt")
p<-list(p1a,p1c,p2a,p2c,p3a,p3c,p4a,p4c,p5a,p5c)
pmat<-as.matrix(cbind(p1a,p1c,p2a,p2c,p3a,p3c,p4a,p4c,p5a,p5c))
pbar<-apply(pmat,1,mean)
pC<-apply(pmat[,seq(2,10,1)],1,mean)
pA<-apply(pmat[,seq(1,9,1)],1,mean)
ids<-read.table("snpids.txt",header=F)
## define mel and stripy
mlb<-4139489
mub<-6414835
slb<-6414836
mel<-c(which(ids[,2]==702.1 & ids[,3] < mlb),
which(ids[,2]==128 & ids[,3] < mub))
stripy<-which(ids[,2]==128 & ids[,3] >= slb)
nsnpMel<-length(mel)
nsnpStripy<-length(stripy)
## for mel, use total change
slist<-as.numeric(names(which(table(ids[,2]) > nsnpMel)))
drop<-which(slist == 702.1 | slist == 128)
slist<-slist[-drop]
nsc<-length(slist)
bnds<-matrix(NA,nrow=nsc,ncol=2)
for(i in 1:nsc){
n<-sum(ids[,2]==slist[i])
bnds[i,1]<-sample(1:(n-nsnpMel),1)
bnds[i,2]<-bnds[i,1]+nsnpMel
}
scmns<-rep(NA,nsc)
scvar<-rep(NA,nsc)
dp<-abs(pbar-p0)
div<-2 * p0 * (1-p0)
for(i in 1:nsc){
Dp<-dp[ids[,2]==slist[i]][bnds[i,1]:bnds[i,2]]
Div<-div[ids[,2]==slist[i]][bnds[i,1]:bnds[i,2]]
scmns[i]<-mean(Dp)
scvar[i]<-mean(Div)
}
o<-lm(scmns ~ scvar)
dpMel<-mean(dp[mel])
divMel<-mean(div[mel])
pdf("ogaDivChange.pdf",width=7,height=7)
plot(scvar,scmns,xlab="mean diversity",ylab="mean change",xlim=c(0.12,0.18),ylim=c(0.027,0.1),pch=19,cex.lab=1.4,cex.axis=1.2)
abline(o$coefficients)
points(divMel,dpMel,col="brown",pch=19)
dev.off()
pdf("ogaDivChangeResid.pdf",width=7,height=7)
plot(scvar,o$residuals,xlab="mean diversity",ylab="residual change",xlim=c(0.12,0.18),ylim=c(-0.003,0.06),pch=19,cex.lab=1.4,cex.axis=1.2)
ex<-o$coefficients[1] + o$coefficients[2] * divMel
points(divMel,dpMel-ex,col="brown",pch=19)
abline(h=0,lty=2)
dev.off()
###########################################
## for stripy, use divergence
slist<-as.numeric(names(which(table(ids[,2]) > nsnpStripy)))
drop<-which(slist == 702.1 | slist == 128)
slist<-slist[-drop]
nsc<-length(slist)
bnds<-matrix(NA,nrow=nsc,ncol=2)
for(i in 1:nsc){
n<-sum(ids[,2]==slist[i])
bnds[i,1]<-sample(1:(n-nsnpStripy),1)
bnds[i,2]<-bnds[i,1]+nsnpStripy
}
scmns<-rep(NA,nsc)
scvar<-rep(NA,nsc)
dp<-abs(pC-pA)
div<-2 * p0 * (1-p0)
for(i in 1:nsc){
Dp<-dp[ids[,2]==slist[i]][bnds[i,1]:bnds[i,2]]
Div<-div[ids[,2]==slist[i]][bnds[i,1]:bnds[i,2]]
scmns[i]<-mean(Dp)
scvar[i]<-mean(Div)
}
o<-lm(scmns ~ scvar)
dpStr<-mean(dp[stripy])
divStr<-mean(div[stripy])
pdf("ogaDivChangeAvC.pdf",width=7,height=7)
plot(scvar,scmns,xlab="mean diversity",ylab="mean change",xlim=c(0.12,0.143),ylim=c(0.0056,0.008),pch=19,cex.lab=1.4,cex.axis=1.2)
abline(o$coefficients)
points(divStr,dpStr,col="forestgreen",pch=19)
dev.off()
pdf("ogaDivChangeResidAvC.pdf",width=7,height=7)
plot(scvar,o$residuals,xlab="mean diversity",ylab="residual change",xlim=c(0.12,0.143),ylim=c(-0.0003,0.0015),pch=19,cex.lab=1.4,cex.axis=1.2)
ex<-o$coefficients[1] + o$coefficients[2] * divStr
points(divStr,dpStr-ex,col="forestgreen",pch=19)
abline(h=0,lty=2)
dev.off()
###########################################