Post date: Feb 13, 2015 7:46:47 PM
MORE TO ADD HERE
From parallelism.R in /home/zgompert/labs/evolution/projects/timema_radiation/popgen/popgen. Key results, significant parallelism
All pairs (minus SM)
obs expect obs/exp p-value
1 2159 9631.407 2.241625e-01 1.000
2 1076 3584.627 3.001707e-01 1.000
3 758 770.772 9.834296e-01 0.706
4 682 104.889 6.502112e+00 0.000
5 448 9.409 4.761399e+01 0.000
6 230 0.540 4.259259e+02 0.000
7 332 0.026 1.276923e+04 0.000
8 309 0.000 Inf 0.000
9 206 0.000 Inf 0.000
Quercus vs. Ceanothus
obs exp obs/exp p-value
1 2664 7419.151 0.3590707 1.000
2 1254 1286.005 0.9751128 0.882
3 638 108.616 5.8739044 0.000
4 486 4.434 109.6075778 0.000
5 261 0.051 5117.6470588 0.000
angiosperm vs. gymnospersm
obs exp obs/exp p-value
1 1897 6894.512 0.2751464 1.000
2 970 995.227 0.9746520 0.847
3 750 62.710 11.9598150 0.000
4 748 1.476 506.7750678 0.000
nl<-25496
hmmDa<-matrix(NA,nrow=nl,ncol=10)
load("hmm2s_BCBOGpar6.rwsp")
hmmDa[,1]<-fit[[1]]
load("hmm2s_BCTURpar6.rwsp")
hmmDa[,2]<-fit[[1]]
load("hmm2s_BMCG3par6.rwsp")
hmmDa[,3]<-fit[[1]]
load("hmm2s_BMTpar6.rwsp")
hmmDa[,4]<-fit[[1]]
load("hmm2s_BSpar6.rwsp")
hmmDa[,5]<-fit[[1]]
load("hmm2s_CRpar6.rwsp")
hmmDa[,6]<-fit[[1]]
load("hmm2s_HFRSpar6.rwsp")
hmmDa[,7]<-fit[[1]]
load("hmm2s_LPpar6.rwsp")
hmmDa[,8]<-fit[[1]]
load("hmm2s_SMpar6.rwsp")
hmmDa[,9]<-fit[[1]]
load("hmm2s_VPpar6.rwsp")
hmmDa[,10]<-fit[[1]]
pdf("noHighDa.pdf",width=7,height=7)
par(mar=c(5.5,5.5,0.5,0.5))
plot(x,type='l',xlab="locus",ylab="no. of pairs",cex.lab=1.5,cex.axis=1.2)
dev.off()
## minus SM, only 'independent pairs'
null<-matrix(NA,nrow=nl,ncol=1000)
hmmDa2<-hmmDa[,-9]
x<-apply(hmmDa2 == 1,1,sum)
for(i in 1:1000){
perm<-hmmDa2
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,])
}
pdf("noHighDaNull.pdf",width=7,height=7)
par(mar=c(5.5,5.5,0.5,0.5))
plot(x,type='l',xlab="locus",ylab="no. of pairs",cex.lab=1.5,cex.axis=1.2)
abline(h=c(mean(mn),mean(q95),mean(q99),mean(q99.9),mean(q1)),col="darkblue")
dev.off()
# by host
qc<-c(1,4,5,7,10)
ag<-c(2,3,6,9)
hmmDaQC<-hmmDa[,qc]
hmmDaAG<-hmmDa[,ag]
xQC<-apply(hmmDaQC == 1,1,sum)
xAG<-apply(hmmDaAG == 1,1,sum)
nullQC<-matrix(NA,nrow=nl,ncol=1000)
nullAG<-matrix(NA,nrow=nl,ncol=1000)
for(i in 1:1000){
perm<-hmmDaQC
for(j in 1:5){
perm[,j]<-sample(perm[,j],nl,replace=F)
}
nullQC[,i]<-apply(perm==1,1,sum)
perm<-hmmDaAG
for(j in 1:4){
perm[,j]<-sample(perm[,j],nl,replace=F)
}
nullAG[,i]<-apply(perm==1,1,sum)
}
obsQC<-rep(NA,5)
expQC<-matrix(NA,nrow=5,ncol=1000)
obsAG<-rep(NA,4)
expAG<-matrix(NA,nrow=4,ncol=1000)
for(i in 1:5){
obsQC[i]<-sum(xQC==i)
expQC[i,]<-apply(nullQC==i,2,sum)
}
psQC<-rep(NA,5)
enrichQC<-rep(NA,5)
for(i in 1:5){
psQC[i]<-mean(expQC[i,] >= obsQC[i])
enrichQC[i]<-obsQC[i]/mean(expQC[i,])
}
for(i in 1:4){
obsAG[i]<-sum(xAG==i)
expAG[i,]<-apply(nullAG==i,2,sum)
}
psAG<-rep(NA,4)
enrichAG<-rep(NA,4)
for(i in 1:4){
psAG[i]<-mean(expAG[i,] >= obsAG[i])
enrichAG[i]<-obsAG[i]/mean(expAG[i,])
}
out<-cbind(obsQC,apply(expQC,1,mean),enrichQC,psQC)
colnames(out)<-c("obs","exp","obs/exp","p-value")
rownames(out)<-1:5
write.csv(out,"parallelQC.csv")
out<-cbind(obsAG,apply(expAG,1,mean),enrichAG,psAG)
colnames(out)<-c("obs","exp","obs/exp","p-value")
rownames(out)<-1:4
write.csv(out,"parallelAG.csv")
mn<-apply(nullQC,1,mean)
q95<-apply(nullQC,1,quantile,probs=0.95)
q99<-apply(nullQC,1,quantile,probs=0.99)
q99.9<-apply(nullQC,1,quantile,probs=0.999)
q1<-apply(nullQC,1,quantile,probs=1)
pdf("noHighDaQC.pdf",width=7,height=7)
par(mar=c(5.5,5.5,0.5,0.5))
plot(xQC,type='l',xlab="locus",ylab="no. of pairs",cex.lab=1.5,cex.axis=1.2)
abline(h=c(mean(mn),mean(q95),mean(q99),mean(q99.9),mean(q1)),col="darkblue")
dev.off()
mn<-apply(nullAG,1,mean)
q95<-apply(nullAG,1,quantile,probs=0.95)
q99<-apply(nullAG,1,quantile,probs=0.99)
q99.9<-apply(nullAG,1,quantile,probs=0.999)
q1<-apply(nullAG,1,quantile,probs=1)
pdf("noHighDaAG.pdf",width=7,height=7)
par(mar=c(5.5,5.5,0.5,0.5))
plot(xAG,type='l',xlab="locus",ylab="no. of pairs",cex.lab=1.5,cex.axis=1.2)
abline(h=c(mean(mn),mean(q95),mean(q99),mean(q99.9),mean(q1)),col="darkblue")
dev.off()