Post date: Oct 05, 2015 9:44:52 PM
I tested for elevated Fst for the 25 DSRs identified in cristinae in the 10 (mostly) non-cristinae taxon pairs. The analysis is described in dsrAnalysis.R which is in /labs/evolution/projects/timema_speciation/timema_radiation/popgen/popgen/. The results, as shown below, are that we actually had a pronounced reduction in Fst in these regions relative to null expectations. I have no idea why this would be (neither a real nor an artefactual reason):
"","obsv","expect","pv"
"BCBOG",0.0169667908,0.0401417101001667,1
"BCTUR",0.092651376,0.125152890826833,1
"BMCG3",0.0135530125,0.0326576725120167,1
"BMT",0.018613236,0.0329723604757333,1
"BS",0.0215572288,0.0327150495832667,0.998
"CR",0.0174027520833333,0.0391629241404333,1
"HFRS",0.1057624964,0.140777714822483,1
"LP",0.2894956748,0.717962868171083,1
"SM",0.3675375436,0.700163761338733,1
"VP",0.0174448296,0.0317644541732167,1
Here is the R code.
fstBCBOG<-as.matrix(read.table("nei20k_timemaBCBOG.txt",header=FALSE))
fstBCTUR<-as.matrix(read.table("nei20k_timemaBCTUR.txt",header=FALSE))
fstBMCG3<-as.matrix(read.table("nei20k_timemaBMCG3.txt",header=FALSE))
fstBMT<-as.matrix(read.table("nei20k_timemaBMT.txt",header=FALSE))
fstBS<-as.matrix(read.table("nei20k_timemaBS.txt",header=FALSE))
fstCR<-as.matrix(read.table("nei20k_timemaCR.txt",header=FALSE))
fstHFRS<-as.matrix(read.table("nei20k_timemaHFRS.txt",header=FALSE))
fstLP<-as.matrix(read.table("nei20k_timemaLP.txt",header=FALSE))
fstSM<-as.matrix(read.table("nei20k_timemaSM.txt",header=FALSE))
fstVP<-as.matrix(read.table("nei20k_timemaVP.txt",header=FALSE))
fst<-cbind(fstBCBOG[,7],fstBCTUR[,7],fstBMCG3[,7],fstBMT[,7],fstBS[,7],fstCR[,7],fstHFRS[,7],fstLP[,7],fstSM[,7],fstVP[,7])
pops<-c("BCBOG","BCTUR","BMCG3","BMT","BS","CR","HFRS","LP","SM","VP")
dsr<-read.csv("tcristDsrs.csv")
dsrindex<-rep(NA,25)
for(i in 1:25){
dsrindex[i]<-which(fstVP[,1]==dsr[i,1] & fstVP[,2]==dsr[i,2])
}
L<-dim(fstVP)[1]
mnfst<-apply(fst,1,mean,na.rm=TRUE)
obs<-mean(mnfst[dsrindex])
N<-1000
null<-rep(NA,N)
for(i in 1:N){
a<-sample(1:L,25,replace=FALSE)
null[i]<-mean(mnfst[a])
}
## each pair
nullm<-matrix(NA,nrow=N,ncol=10)
obsv<-apply(fst[dsrindex,],2,mean,na.rm=TRUE)
for(j in 1:10){
for(i in 1:N){
a<-sample(1:L,25,replace=FALSE)
nullm[i,j]<-mean(fst[a,j],na.rm=TRUE)
}
}
pv<-rep(NA,10)
expect<-apply(nullm,2,mean)
for(j in 1:10){
pv[j]<-mean(nullm[,j] >= obsv[j])
}
cbind(obsv,expect,pv)
rownames(out)<-pops
write.csv(out,"dsrFst.csv")