Post date: Jun 25, 2015 8:8:12 PM
Still nothing significant for pattern or the other traits (taken together) in term of Fst in the natural populations (I even tried excluding the LAxPRC based on the ABC results). However, the test for pattern was basically garbage. Nearly all of the top pip SNPs were on LG8, but on scaffolds that moved off of LG 8 (and off of an LG period) with the revised genome. Stripe probably would have mapped to other current LG 8 had these orphan scaffolds not been there. So, I am not sure we can really say much. I am not sure what to do about this, but I will think a bit more.
Here is the R code:
n1<-22112
## trait pips
TrBW<-as.matrix(read.table("pip_2013FHA_all_residBW_MAF01_sparse.txt",header=F))
TrHW<-as.matrix(read.table("pip_2013FHA_all_residHW_MAF01_sparse.txt",header=F))
TrlatGB<-as.matrix(read.table("pip_2013FHA_all_resid_latGB_MAF01_sparse.txt",header=F))
TrlatL<-as.matrix(read.table("pip_2013FHA_all_resid_latL_MAF01_sparse.txt",header=F))
TrlatRG<-as.matrix(read.table("pip_2013FHA_all_resid_latRG_MAF01_sparse.txt",header=F))
Trfpent<-as.matrix(read.table("pip_2013FHA_females_penta_MAF01_sparse.txt",header=F))
Trfhepta<-as.matrix(read.table("pip_2013FHA_females_hepta_MAF01_sparse.txt",header=F))
Trfnona<-as.matrix(read.table("pip_2013FHA_females_nona_MAF01_sparse.txt",header=F))
Trmpent<-as.matrix(read.table("pip_2013FHA_males_mpenta_MAF01_sparse.txt",header=F))
Trmhepta<-as.matrix(read.table("pip_2013FHA_males_mhepta_MAF01_sparse.txt",header=F))
Trmnona<-as.matrix(read.table("pip_2013FHA_males_mnona_MAF01_sparse.txt",header=F))
Trpatln<-as.matrix(read.table("pip_2013FHA_all_resid_percent_stripe_MAF01_sparse.txt",header=F))
Tlist<-list(TrBW,TrHW,TrlatGB,TrlatL,TrlatRG,Trfpent,Trfhepta,Trfnona,Trmpent,Trmhepta,Trmnona,Trpatln)
## fst
parname<-c("pi","Dxy","Da","Fst")
d1<-matrix(scan("../popgen/nei20k_tcrist_MR1",n=n1*7,sep=" "),nrow=n1,ncol=7,byrow=T)
d2<-matrix(scan("../popgen/nei20k_tcrist_R12",n=n1*7,sep=" "),nrow=n1,ncol=7,byrow=T)
d3<-matrix(scan("../popgen/nei20k_tcrist_L",n=n1*7,sep=" "),nrow=n1,ncol=7,byrow=T)
d4<-matrix(scan("../popgen/nei20k_tcrist_HV",n=n1*7,sep=" "),nrow=n1,ncol=7,byrow=T)
parm<-7
# get high pip pattern SNPs
hiPip<-NULL
bnd<-0.05
#bnd<-0.2
hiPip<-Trpatln[which(Trpatln[,4] > bnd),]
uHiPip<-hiPip[which(!duplicated(hiPip[,3])),]
# get high pip SNPs
hiPip<-NULL
bnd<-0.05
for(i in 1:12){
x<-Tlist[[i]][which(Tlist[[i]][,4] > bnd),]
if(length(x) > 0){hiPip<-rbind(hiPip,x)}
}
uHiPip<-hiPip[-duplicated(hiPip[,3]),]
## 20kb window for each hi pip SNPs
n<-dim(uHiPip)[1]
pipWin<-rep(NA,n)
for(i in 1:n){
x<-which(d1[,1]==uHiPip[i,2] & abs(d1[,2]-uHiPip[i,3]) < 20000)
if (length(x) == 1){
pipWin[i]<-x
}
else if (length(x) > 1){
x<-which(d1[,1]==uHiPip[i,2] & abs(d1[,2]-uHiPip[i,3]) < 10000)
if (length(x) == 1){
pipWin[i]<-x
}
else{
cat ("mult snps for ",i,"\n")
}
}
}
pipWin<-pipWin[is.na(pipWin)==FALSE]
A<-as.matrix(cbind(d1[,parm],d2[,parm],d3[,parm],d4[,parm]))
mnFst<-apply(A[,-3],1,mean,na.rm=TRUE)
obs<-mean(mnFst[pipWin])
null<-rep(NA,10000)
for(i in 1:10000){
a<-sample(1:n1,length(pipWin),replace=TRUE)
null[i]<-mean(mnFst[a])
}
p<-mean(null > obs)