## R code to calcualte SNP-specific Fst, infiles available upon request ## original hudson Fst = 1 - Hw/Hb hudsonFst2<-function(p1=NA,p2=NA,n1=NA,n2=NA){ numerator<-p1 * (1 - p1) + p2 * (1 - p2) denominator<-p1 * (1 - p2) + p2 * (1 - p1) fst<-1 - numerator/denominator out<-cbind(numerator,denominator,fst) return(out) } G<-matrix(scan("../variant/pntest_wildwgVarsA.txt",n=160 * 4391556,sep=" "),nrow=4391556,ncol=160,byrow=TRUE) pops<-read.table("pops.txt",header=F) plist<-unique(as.character(pops[,1])) pops<-as.character(pops[,1]) nloc<-dim(G)[1] npop<-length(plist) locinfo<-read.table("locids.txt",header=F) ## population allele frequencies af<-matrix(NA,nrow=nloc,ncol=npop) for (k in 1:npop){ af[,k]<-apply(G[,which(pops==plist[k])],1,mean)/2 } N<-table(pops) ## fst for each snp fstHVAxHVChu<-hudsonFst2(p1=af[,1],p2=af[,2],n1=N[1],n2=N[2]) fstMR1AxMR1Chu<-hudsonFst2(p1=af[,4],p2=af[,5],n1=N[4],n2=N[5]) fstR12AxR12Chu<-hudsonFst2(p1=af[,7],p2=af[,8],n1=N[7],n2=N[8]) fstLAxPRChu<-hudsonFst2(p1=af[,3],p2=af[,6],n1=N[3],n2=N[6])