Post date: Aug 18, 2015 4:23:0 PM
We decided to report results on Fst for the trait-associated SNPs from the GBS-based GWA, so I re-ran (to make sure I know what was what) and re-summarized the results. The R code I used is here (/labs/evolution/projects/timema_speciation/timema_wgexperiment/timemaTraitsQTN/traitFstRandomTest.R):
n1<-22112
## trait pips
TrBL<-as.matrix(read.table("pip_2013_BL_.txt",header=F))
TrBW<-as.matrix(read.table("pip_2013_BW_.txt",header=F))
TrHW<-as.matrix(read.table("pip_2013_HW_.txt",header=F))
TrlatGB<-as.matrix(read.table("pip_2013_latGB_.txt",header=F))
TrlatL<-as.matrix(read.table("pip_2013_latL_.txt",header=F))
TrlatRG<-as.matrix(read.table("pip_2013_latRG_.txt",header=F))
Trfpent<-as.matrix(read.table("pip_2013_fpent_.txt",header=F))
Trfhepta<-as.matrix(read.table("pip_2013_fhept_.txt",header=F))
Trfnona<-as.matrix(read.table("pip_2013_fnona_.txt",header=F))
Trmpent<-as.matrix(read.table("pip_2013_mpent_.txt",header=F))
Trmhepta<-as.matrix(read.table("pip_2013_mhept_.txt",header=F))
Trmnona<-as.matrix(read.table("pip_2013_mnona_.txt",header=F))
Trpatln<-as.matrix(read.table("pip_2013_patln_.txt",header=F))
Tlist<-list(TrBL,TrBW,TrHW,TrlatGB,TrlatL,TrlatRG,Trfpent,Trfhepta,Trfnona,Trmpent,Trmhepta,Trmnona,Trpatln)
# get high pip SNPs
hiPip<-NULL
bnd<-0.05
#for(i in c(1:6,13)){
for(i in c(4:6,13)){
#for(i in c(7:12)){
#for(i in 13){
#for(i in 1:13){
x<-Tlist[[i]][which(Tlist[[i]][,4] > bnd),]
if(length(x) > 0){hiPip<-rbind(hiPip,x)}
}
uHiPip<-hiPip[!duplicated(hiPip[,3]),]
## fst
parname<-c("pi","Dxy","Da","Fst")
d1<-matrix(scan("nei20k_tcrist_MR1",n=n1*7,sep=" "),nrow=n1,ncol=7,byrow=T)
d2<-matrix(scan("nei20k_tcrist_R12",n=n1*7,sep=" "),nrow=n1,ncol=7,byrow=T)
d3<-matrix(scan("nei20k_tcrist_L",n=n1*7,sep=" "),nrow=n1,ncol=7,byrow=T)
d4<-matrix(scan("nei20k_tcrist_HV",n=n1*7,sep=" "),nrow=n1,ncol=7,byrow=T)
parm<-7
## 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,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)
Note that different sets of traits were analyzed by commenting/uncommenting the first for loop. Here are the results. I have results for different sets of traits and cut-offs. First are the results for all 13 traits with any SNPs with PIP > 5%, this is followed by the results for those with PIP > 10%. In each case I list the number of SNPs, the number of those that were on an LG in the GWAS analysis (the analysis appears to have included SNPs with and without LGs even for the version 3 of the genome) and then the number of those that have Fst estimates (i.e. were on a LG for version 2 and on big enough scaffold for a 20 kb window). The last number is the actual number analyzed. In each case I give the mean Fst (obs), mean null, 95th quantile of null, enrichment (i.e. obs/mean(null)), and the p-value for obs > null. After the 13 traits I report the results for the six CHC traits, all morphological traits (i.e. everything but CHCs), and color and stripe. I tried to run just stripe but we don't end up with any SNPs with PIP > 5% that we have Fst for (which means the color and stripe result is really just color). I hadn't recalled this detail from before, but it means we really aren't testing stripe (not without dropping below PIP > 5%, which would be hard to justify). In general, we don't see much, but there is a trend towards higher Fst when you focus on the morphology and color traits.
-all 13 traits, pip > 0.05; 59 total, 43 with LGs, 34 with Fst
obs = 0.0101
null = 0.0101
nullq95 = 0.0109
enrich = 0.997
p = 0.4965
-all 13 traits, pip > 0.1; 25 total, 16 with LGs, 11 with Fst
obs = 0.0095
null = 0.0101
nullq95 = 0.0115
enrich = 0.944
p = 0.7731
-CHCs, pip > 0.05; 33 total, 26 with LGs, 23 with Fst
obs = 0.0097
null = 0.0101
nullq95 = 0.0111
enrich = 0.961
p = 0.764
-morphology, pip > 0.05; 26 total, 17 with LGs, 11 with Fsts
obs = 0.0108
null = 0.0101
nullq95 = 0.0114
enrich = 1.073
p = 0.149
- lateral color and stripe, pip > 0.05; 17 total, 11 with LGs, 4 with Fsts
obs = 0.0113
null = 0.0101
nullq95 = 0.0122
enrich = 1.121
p = 0.119
- stripe, pip > 0.05; 4 total, 4 with LGs, 0 with Fst
obs = NA
null = NA
nullq95 = NA
enrich = NA
p = NA