Post date: Jan 08, 2017 3:10:28 AM
We had previously looked at change as function of (and controlling for) diversity in the OGA between generation experiments for each replicate bush. This included five transplants each to A and C. I decided to also look at scaffold 128 and 702.1 for change between A and C. We don't have starting frequencies per se, so I just took average freq. on A vs. C. The code is in the commands.R file, but the most relevant bits (new code are below). Here are the results: mean dif, 95th quantile dif. 128 clearly stands out, the signal for 702.1 is weaker.
## 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)
ids<-read.table("snpids.txt",header=F)
## calculate allele frequency change
dp<-vector("list",10)
scmns<-vector("list",10)
scq95<-scmns
scn<-scmns
scvar<-scmns
for(i in 1:10){
dp[[i]]<-abs(p[[i]]-p0)
scmns[[i]]<-tapply(X=dp[[i]],INDEX=ids[,2],mean)
scq95[[i]]<-tapply(X=dp[[i]],INDEX=ids[,2],quantile,probs=0.95)
scvar[[i]]<-tapply(X=2 * p0 * (1 - p0),INDEX=ids[,2],mean)
scn[[i]]<-tapply(X=is.na(dp[[i]])==FALSE,INDEX=ids[,2],sum)
}
## divergence between C and A
dif<-abs((p1a+p2a+p3a+p4a+p5a)/5-(p1c+p2c+p3c+p4c+p5c)/5)
scdif<-tapply(X=dif,INDEX=ids[,2],mean)
scdif95<-tapply(X=dif,INDEX=ids[,2],quantile,probs=0.95)
## diversity by change
pdf("ogaDivMeanChange.pdf",width=4,height=10)
par(mfrow=c(5,2))
par(mar=c(4,4,2,2))
for(i in 1:10){
a<-which(scn[[i]] > 500)
x<-cbind(scvar[[i]][a],scmns[[i]][a])
plot(x[,1],x[,2],ylim=c(0.0,0.15),xlim=c(0.1,0.18),xlab="diversity (2pq)",ylab="mean freq. change",cex.lab=1.1,cex.axis=1)
o<-lm(x[,2] ~ x[,1])
abline(o$coefficients)
title(main=paste(trt[i],", ",re[i]))
## scaf 128
q128[i]<-x[3,2]
points(x[3,1],x[3,2],col="blue",pch=19)
## scaf 702.1
q702[i]<-x[15,2]
points(x[15,1],x[15,2],col="red",pch=19)
## left and right
mn128l[i]<-mean(dp[[i]][left])
mn128r[i]<-mean(dp[[i]][right])
var128l<-mean(2 * p0[left] * (1-p0[left]))
var128r<-mean(2 * p0[right] * (1-p0[right]))
points(var128l,mn128l[i],col="green",pch=19)
points(var128r,mn128r[i],col="orange",pch=19)
}
for(i in 1:10){ ## residuals
a<-which(scn[[i]] > 500)
o<-lm(scmns[[i]][a] ~ scvar[[i]][a])
x<-cbind(scvar[[i]][a],o$residuals)
plot(x[,1],x[,2],ylim=c(-0.02,0.07),xlim=c(0.1,0.18),xlab="diversity (2pq)",ylab="residual change",cex.lab=1.1,cex.axis=1)
title(main=paste(trt[i],", ",re[i]))
## scaf 128
q128[i]<-x[3,2]
points(x[3,1],x[3,2],col="blue",pch=19)
## scaf 702.1
q702[i]<-x[15,2]
points(x[15,1],x[15,2],col="red",pch=19)
## left and right
mn128l[i]<-mean(dp[[i]][left])
mn128r[i]<-mean(dp[[i]][right])
var128l<-mean(2 * p0[left] * (1-p0[left]))
var128r<-mean(2 * p0[right] * (1-p0[right]))
ex<-o$coefficients[1] + o$coefficients[2] * var128l
points(var128l,mn128l[i]-ex,col="green",pch=19)
ex<-o$coefficients[1] + o$coefficients[2] * var128r
points(var128r,mn128r[i]-ex,col="orange",pch=19)
abline(h=0,lty=2)
}
dev.off()
pdf("ogaDivDif95.pdf",width=5,height=10)
par(mfrow=c(2,1))
par(mar=c(4,4,2,2))
a<-which(scn[[i]] > 500)
x<-cbind(scvar[[i]][a],scdif95[a])
plot(x[,1],x[,2],ylim=c(0.054,0.11),xlim=c(0.1,0.18),xlab="diversity (2pq)",ylab="95q divergence",cex.lab=1.1,cex.axis=1)
o<-lm(x[,2] ~ x[,1])
abline(o$coefficients)
## scaf 128
points(x[3,1],x[3,2],col="blue",pch=19)
## scaf 702.1
points(x[15,1],x[15,2],col="red",pch=19)
## left and right
mn128l<-quantile(dif[left],probs=0.95)
mn128r<-quantile(dif[right],probs=0.95)
var128l<-mean(2 * p0[left] * (1-p0[left]))
var128r<-mean(2 * p0[right] * (1-p0[right]))
points(var128l,mn128l,col="green",pch=19)
points(var128r,mn128r,col="orange",pch=19)
# residuals
o<-lm(scdif95[a] ~ scvar[[1]][a])
x<-cbind(scvar[[i]][a],o$residuals)
plot(x[,1],x[,2],ylim=c(-0.013,0.028),xlim=c(0.1,0.18),xlab="diversity (2pq)",ylab="residual change",cex.lab=1.1,cex.axis=1)
## scaf 128
points(x[3,1],x[3,2],col="blue",pch=19)
## scaf 702.1
points(x[15,1],x[15,2],col="red",pch=19)
## left and right
mn128l<-quantile(dif[left],probs=0.95)
mn128r<-quantile(dif[right],probs=0.95)
var128l<-mean(2 * p0[left] * (1-p0[left]))
var128r<-mean(2 * p0[right] * (1-p0[right]))
ex<-o$coefficients[1] + o$coefficients[2] * var128l
points(var128l,mn128l-ex,col="green",pch=19)
ex<-o$coefficients[1] + o$coefficients[2] * var128r
points(var128r,mn128r-ex,col="orange",pch=19)
abline(h=0,lty=2)
dev.off()