Post date: Oct 23, 2015 9:54:18 PM
We are now redoing analyses to identify high Fst 20kb windows from the natural populations and experiment (between survivors), and test for an excess relative to null expectations. First, I am re-estimating Fst for windows as defined by the T. cristinae ecotypes. The analyses and inputs are in /labs/evolution/projects/timema_speciation/tcristinae/data, but the results will end up in /labs/evolution/projects/timema_speciation/tcristinae/fst (i.e. nei* files). The scripts calcNeiExperiment.pl generates the two R scripts that perform the calculations (sourceExp0.R for initial pops, and sourceExp1.R for survivors). The second of these scripts is pasted below:
nl<-8151262
nc<-3
llist<-matrix(scan("experimentLocList.txt",n=nl*4,sep=" "),nrow=nl,ncol=4,byrow=T)
scafs<-unique(llist[,3])
p1<-matrix(scan("p_timemaAsurv.txt",n=nl*nc,what="character",sep=" "),nrow=nl,ncol=nc,byrow=T)
p2<-matrix(scan("p_timemaCsurv.txt",n=nl*nc,what="character",sep=" "),nrow=nl,ncol=nc,byrow=T)
pavg<-(as.numeric(p1[,3]) + as.numeric(p2[,3]))/2
pi1<-2 * as.numeric(p1[,3]) * (1 - as.numeric(p1[,3]))
pi2<-2 * as.numeric(p2[,3]) * (1 - as.numeric(p2[,3]))
piavg<-2 * pavg * (1 - pavg)
dxy<-as.numeric(p1[,3]) * (1 - as.numeric(p2[,3])) + (1 - as.numeric(p1[,3])) * as.numeric(p2[,3])
regs<-read.table("tcrist20kbRegionsGenomeV03.txt",header=T)
for(i in 1:dim(regs)[1]){
x<-which(llist[,3]==regs[i,3]);pi1s<-pi1[x];pi2s<-pi2[x];dxys<-dxy[x];piavgs<-piavg[x]
lls<-llist[x,4]
window<-20000
lb<-regs[i,2]-10000
ub<-regs[i,2]+10000
a<-which(lls >= lb & lls < ub)
if(length(a) >= 1){
pi<-sum((pi1s[a]+pi2s[a])/2) * 1/window
D<-sum(dxys[a]) * 1/window
pit<-sum(piavgs[a]) * 1/window
fst<-(pit - pi)/pit
Da<-D-pi
nei<-cbind(regs[i,],length(a),pi,D,Da,fst)
write.table(round(nei,8),file="nei20k_AsurvCsurv.txt",row.names=F,col.names=F,quote=F,append=TRUE)
}
else{write.table(t(rep(NA,9)),file="nei20k_AsurvCsurv.txt",row.names=F,col.names=F,quote=F,append=TRUE)}
}
Testing for DSRs...