Post date: May 25, 2015 7:21:32 PM
I classified all 20kb windows with pi > 95th empirical quantile within each population as high pi (note I dropped SM as we have done in the past, so we have 9 populations). I then enumerated the number of populations (0 to 9) that each 20 kb window was high pi in. I used a randomization test (1000 randomizations) to generate the null expectations for parallel high pi. The results are in /labs/evolution/projects/timema_radiation/popgen/popgen/piParTab.csv. The file includes the number of 20 kb windows that were hi pi in 1, ... , 9 populations (obs), as well as the enrichment (observed / mean of the null), the median 95th quantile and 99th quantile of the null distribution, and the p-value (ps) for observed > expected (one-tailed).
The R code for this (piQuantParallelism.R) is here:
nl<-25496
pi<-matrix(NA,nrow=nl,ncol=10)
load("hmm2s_BCBOGpar4.rwsp")
pi[,1]<-nei[,4]
load("hmm2s_BCTURpar4.rwsp")
pi[,2]<-nei[,4]
load("hmm2s_BMCG3par4.rwsp")
pi[,3]<-nei[,4]
load("hmm2s_BMTpar4.rwsp")
pi[,4]<-nei[,4]
load("hmm2s_BSpar4.rwsp")
pi[,5]<-nei[,4]
load("hmm2s_CRpar6.rwsp")
pi[,6]<-nei[,4]
load("hmm2s_HFRSpar6.rwsp")
pi[,7]<-nei[,4]
load("hmm2s_LPpar6.rwsp")
pi[,8]<-nei[,4]
load("hmm2s_SMpar6.rwsp")
pi[,9]<-nei[,4]
load("hmm2s_VPpar6.rwsp")
pi[,10]<-nei[,4]
## summarize hi pi
## minus SM, only 'independent pairs'
pi2<-pi[,-9]
qs<-apply(pi,2,quantile,probs=0.95)
pihi<-pi
for(i in 1:9){
pihi[,i]<-pi2[,i] >= qs[i]
}
null<-matrix(NA,nrow=nl,ncol=1000)
x<-apply(pihi == 1,1,sum)
for(i in 1:1000){
perm<-pihi
for(j in 1:9){
perm[,j]<-sample(perm[,j],nl,replace=F)
}
null[,i]<-apply(perm==1,1,sum)
}
mn<-apply(null,1,mean)
q95<-apply(null,1,quantile,probs=0.95)
q99<-apply(null,1,quantile,probs=0.99)
q99.9<-apply(null,1,quantile,probs=0.999)
q1<-apply(null,1,quantile,probs=1)
obs<-rep(NA,9)
exp<-matrix(NA,nrow=9,ncol=1000)
for(i in 1:9){
obs[i]<-sum(x==i)
exp[i,]<-apply(null==i,2,sum)
}
ps<-rep(NA,9)
enrich<-rep(NA,9)
for(i in 1:9){
ps[i]<-mean(exp[i,] >= obs[i])
enrich[i]<-obs[i]/mean(exp[i,])
}
expq<-apply(exp,1,quantile,probs=c(0.5,0.95,0.99))
pdf("piPar.pdf",width=8,height=5)
par(mar=c(5.5,5.5,0.5,0.5))
barplot(obs,width=0.8,space=0.25,names=1:9,ylab="no. of 20 kb windows",xlab="no. of loci",cex.lab=1.6,cex.axis=1.3)
dev.off()
out<-cbind(obs,enrich,t(expq),ps)
rownames(out)<-1:9
write.csv(out,"piParTab.csv",quote=FALSE)