Post date: Jan 06, 2017 8:29:25 PM
Given the diversity change relationship we documented here, I wanted to see whether 128-left (or 702.1 or 128 in general) change more than you would expect for its diversity levels. I did this in two ways:
1. Inferring Fst, which standardizes divergence by diversity.
2. Examining the residuals from regressing change on diversity
In either case 128-left stands out: Fst, residual mean change, residual 95th q change
p11<-scan("pEM_fha2011.txt")
p13<-scan("pEM_fha2013.txt")
ids<-read.table("snpids.txt",header=F)
dp<-abs(p13-p11)
div<-2 * p11 * (1-p11)
div13<-2 * p13 * (1-p13)
scmns<-tapply(X=dp,INDEX=ids[,2],mean)
scvar<-tapply(X=div,INDEX=ids[,2],mean)
scvar13<-tapply(X=div13,INDEX=ids[,2],mean)
scn<-tapply(X=is.na(dp)==FALSE,INDEX=ids[,2],sum)
##
pbar<-(p11+p13)/2
piavg<-2 * pbar * (1 -pbar)
pisub<-(div+div13)/2
fstn<-tapply(X=(piavg-pisub),INDEX=ids[,2],mean)
fstd<-tapply(X=piavg,INDEX=ids[,2],mean)
fst<-fstn/fstd
## 128 left and right
left128<-which(ids[,2]==128 & ids[,3] < lb)
right128<-which(ids[,2]==128 & ids[,3] > ub)
sc128<-which(ids[,2]==128)
sc702<-which(ids[,2]==702.1)
fst128l<-mean(piavg[left128]-pisub[left128])/mean(piavg[left128])
fst128r<-mean(piavg[right128]-pisub[right128])/mean(piavg[right128])
fst128<-mean(piavg[sc128]-pisub[sc128])/mean(piavg[sc128])
fst702<-mean(piavg[sc702]-pisub[sc702])/mean(piavg[sc702])
pdf("fstFHA.pdf",width=6,height=6)
hist(fst[a],xlim=c(0,0.006),main="",xlab="Fst")
abline(v=fst128,col="blue")
abline(v=fst128l,col="green")
abline(v=fst128r,col="orange")
abline(v=fst702,col="red")
dev.off()
## drop small scaffolds, too noisy
a<-which(scn > 500)
x<-cbind(scn[a],scmns[a])
## plot of mean change
pdf("fhaScafChange.pdf",width=7,height=7)
plot(x[,1],x[,2],xlab="no. SNPs",ylab="mean freq. change",cex.lab=1.4,cex.axis=1.2,ylim=c(0,0.028))
## scaf 128
points(x[3,1],x[3,2],col="blue",pch=19)
## scaf 702.1
points(x[13,1],x[13,2],col="red",pch=19)
lb<-5000000
ub<-7000000
left128<-which(ids[,2]==128 & ids[,3] < lb)
right128<-which(ids[,2]==128 & ids[,3] > ub)
mnl<-mean(dp[left128])
mnr<-mean(dp[right128])
points(length(left128),mnl,col="green",pch=19)
points(length(right128),mnr,col="orange",pch=19)
legend(15000,0.028,c("702.1","128","128-left","128-right"),pch=19,col=c("red","blue","green","orange"))
dev.off()
## plot of 95thq change
sc95q<-tapply(X=dp,INDEX=ids[,2],quantile,probs=0.95)
x<-cbind(scn[a],sc95q[a])
pdf("fhaScafChange95.pdf",width=7,height=7)
plot(x[,1],x[,2],xlab="no. SNPs",ylab="95th q. freq. change",cex.lab=1.4,cex.axis=1.2,ylim=c(0,0.1))
## scaf 128
points(x[3,1],x[3,2],col="blue",pch=19)
## scaf 702.1
points(x[13,1],x[13,2],col="red",pch=19)
lb<-5000000
ub<-7000000
left128<-which(ids[,2]==128 & ids[,3] < lb)
right128<-which(ids[,2]==128 & ids[,3] > ub)
mnl<-quantile(dp[left128],probs=0.95)
mnr<-quantile(dp[right128],probs=0.95)
points(length(left128),mnl,col="green",pch=19)
points(length(right128),mnr,col="orange",pch=19)
legend(15000,0.028,c("702.1","128","128-left","128-right"),pch=19,col=c("red","blue","green","orange"))
dev.off()
## variation vs change
varl<-mean(div[left128])
varr<-mean(div[right128])
varl13<-mean(div13[left128])
varr13<-mean(div13[right128])
o<-lm(scmns[a] ~ scvar[a])
pdf("fhaDivChange.pdf",width=7,height=7)
plot(scvar[a],scmns[a],xlab="mean diversity",ylab="mean change",xlim=c(0.12,0.21),ylim=c(0.013,0.028),pch=19,cex.lab=1.4,cex.axis=1.2)
abline(o$coefficients)
points(scvar[a][13],scmns[a][13],col="red",pch=19)
points(scvar[a][3],scmns[a][3],col="blue",pch=19)
points(varr,mnr,col="orange",pch=19)
points(varl,mnl,col="green",pch=19)
legend(0.124,0.028,c("702.1","128","128-left","128-right"),pch=19,col=c("red","blue","green","orange"))
dev.off()
o<-lm(sc95q[a] ~ scvar[a])
pdf("fhaDiv95Change.pdf",width=7,height=7)
plot(scvar[a],sc95q[a],xlab="mean diversity",ylab="95th q. change",xlim=c(0.12,0.21),ylim=c(0.013,0.1),pch=19,cex.lab=1.4,cex.axis=1.2)
abline(o$coefficients)
points(scvar[a][13],sc95q[a][13],col="red",pch=19)
points(scvar[a][3],sc95q[a][3],col="blue",pch=19)
points(varr,mnr,col="orange",pch=19)
points(varl,mnl,col="green",pch=19)
legend(0.124,0.1,c("702.1","128","128-left","128-right"),pch=19,col=c("red","blue","green","orange"))
dev.off()
## residual plots
o<-lm(scmns[a] ~ scvar[a])
mnl<-mean(dp[left128])
mnr<-mean(dp[right128])
pdf("fhaDivChangeResid.pdf",width=7,height=7)
plot(scvar[a], o$residuals,xlab="mean diversity",ylab="residual change",cex.lab=1.4,ylim=c(-0.002,0.0037),xlim=c(0.12,0.21))
points(scvar[a][13],o$residuals[13],col="red",pch=19)
points(scvar[a][3],o$residuals[3],col="blue",pch=19)
ex<-o$coefficients[1] + o$coefficients[2] * varr
points(varr,mnr-ex,col="orange",pch=19)
ex<-o$coefficients[1] + o$coefficients[2] * varl
points(varl,mnl-ex,col="green",pch=19)
abline(h=0,lty=2)
legend(0.195,-0.0002,c("702.1","128","128-left","128-right"),pch=19,col=c("red","blue","green","orange"))
dev.off()
## with 95th quantile
o<-lm(sc95q[a] ~ scvar[a])
mnl<-quantile(dp[left128],probs=0.95)
mnr<-quantile(dp[right128],probs=0.95)
pdf("fhaDivChange95Resid.pdf",width=7,height=7)
plot(scvar[a], o$residuals,xlab="mean diversity",ylab="residual change (95th q.)",cex.lab=1.4,ylim=c(-0.01,0.019),xlim=c(0.12,0.21))
points(scvar[a][13],o$residuals[13],col="red",pch=19)
points(scvar[a][3],o$residuals[3],col="blue",pch=19)
ex<-o$coefficients[1] + o$coefficients[2] * varr
points(varr,mnr-ex,col="orange",pch=19)
ex<-o$coefficients[1] + o$coefficients[2] * varl
points(varl,mnl-ex,col="green",pch=19)
abline(h=0,lty=2)
legend(0.195,-0.0002,c("702.1","128","128-left","128-right"),pch=19,col=c("red","blue","green","orange"))
dev.off()
Here is the full code (see genomicChange.R in /uufs/chpc.utah.edu/common/home/u6000989/projects/timema_fluct/genomic_change_dark_morph/popgen_fha/).