Post date: Mar 09, 2015 2:39:58 AM
We have done some analyses already looking at mean Fst in 20kb windows (one or more) centered on two SNPs that were associated with stripe in the survivors from the Ecology Letters experiments. Now I want to refine these analyses by calculating posterior probabilities for at least one associated SNP in each 20kb window. Note that the windows are a bit different for the four T. cristinae populations from the Science paper and the other 10 pairs as they use different versions of the genome. The R code for this is below and the results will be in /labs/evolution/projects/timema_radiation/popgen/popgen/stripeData/.
## read in gwas results for stripe from gemma
nl<-3743817
pip<-scan("stripePips.txt",n=nl)
scafPos<-matrix(scan("stripeLocusList.txt",n=nl*2,sep=" "),nrow=nl,ncol=2,byrow=TRUE)
## T crist region list
reg<-as.matrix(read.table("../regionTcr20klist.txt",header=F))
nr<-dim(reg)[1]
regPip<-rep(NA,nr)
for(i in 1:nr){
scaf<-reg[i,1]
lb<-reg[i,2]-9999
ub<-reg[i,2]+10000
X<-which(scafPos[,1]==scaf & scafPos[,2] >= lb & scafPos[,2] <= ub)
regPip[i]<-1-prod(1-pip[X])
}
write.table(regPip,"regionPipsTcrist.txt",row.names=F,col.names=F,quote=F)
## Timema region list
reg<-as.matrix(read.table("../region20klist.txt",header=F))
nr<-dim(reg)[1]
regPip<-rep(NA,nr)
for(i in 1:nr){
scaf<-reg[i,1]
lb<-reg[i,2]-9999
ub<-reg[i,2]+10000
X<-which(scafPos[,1]==scaf & scafPos[,2] >= lb & scafPos[,2] <= ub)
regPip[i]<-1-prod(1-pip[X])
}
write.table(regPip,"regionPipsTimema.txt",row.names=F,col.names=F,quote=F)