Post date: Mar 03, 2017 3:53:8 AM
I recalculated genomic change (and a null distribution) for FHA with the defined locus mel (see notes here). Right now the results and code are in gompert:/home/zgompert/Local/timemaGenomicChange as kingspeak is down for maintenance.
There are three plots summarizing genomic change in FHA for mel and one mel size region from each of 40 scaffolds that are as big or bigger than mel. These 40 points provide a null distribution. We could generate a more substantial null distribution by re-sampling within these 40 scaffolds, but in some ways this is nicer (as we discussed, there are fully independent).
In each case the points (or histogram) shows the null, with the brown dot or line showing mel. The plots are for Fst (histogram), diversity vs. change and residual change (shown against diversity). The last two are showing the same thing, just in slightly different ways. In each case, mel is the most extreme value. Thus, the p-value for standing out beyond the null is 1/41 = 0.024.
stripy doesn't stand out (and is not included), but we don't expect it to.
Here is the R code.
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)
pbar<-(p11+p13)/2
piavg<-2 * pbar * (1 -pbar)
pisub<-(div+div13)/2
## 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)
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
}
fst<-rep(NA,nsc)
scmns<-rep(NA,nsc)
scvar<-rep(NA,nsc)
scvar13<-rep(NA,nsc)
for(i in 1:nsc){
PIavg<-piavg[ids[,2]==slist[i]][bnds[i,1]:bnds[i,2]]
PIsub<-pisub[ids[,2]==slist[i]][bnds[i,1]:bnds[i,2]]
Dp<-dp[ids[,2]==slist[i]][bnds[i,1]:bnds[i,2]]
Div<-div[ids[,2]==slist[i]][bnds[i,1]:bnds[i,2]]
Div13<-div13[ids[,2]==slist[i]][bnds[i,1]:bnds[i,2]]
fst[i]<-mean(PIavg-PIsub)/mean(PIavg)
scmns[i]<-mean(Dp)
scvar[i]<-mean(Div)
scvar13[i]<-mean(Div13)
}
fstMel<-mean(piavg[mel]-pisub[mel])/mean(piavg[mel])
fstStr<-mean(piavg[stripy]-pisub[stripy])/mean(piavg[stripy])
pdf("fstFHA.pdf",width=6,height=6)
hist(fst,xlim=c(0,0.006),main="",xlab="Fst",col="gray")
abline(v=fstMel,col="brown",lwd=2)
dev.off()
o<-lm(scmns ~ scvar)
dpMel<-mean(dp[mel])
divMel<-mean(div[mel])
div13Mel<-mean(div13[mel])
pdf("fhaDivChange.pdf",width=7,height=7)
plot(scvar,scmns,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(divMel,dpMel,col="brown",pch=19)
dev.off()
## residual plots
o<-lm(scmns~ scvar)
pdf("fhaDivChangeResid.pdf",width=7,height=7)
plot(scvar, o$residuals,xlab="mean diversity",ylab="residual change",cex.lab=1.4,ylim=c(-0.002,0.0062),xlim=c(0.12,0.21),pch=19)
ex<-o$coefficients[1] + o$coefficients[2] * divMel
points(divMel,dpMel-ex,col="brown",pch=19)
abline(h=0,lty=2)
dev.off()