Post date: Oct 11, 2013 5:12:4 PM
We wanted to know whether SNPs associated with the stripe phenotype among the recaptured Timema in the within-generation transplant experiment were affected by selection. I posed this question in two different ways. First, I asked whether any of the 114 exceptional change SNPs from the Ceanothus treatment were associated with stripedness. I did this by calculating the allele frequency difference between striped and unstriped individuals for each of these 114 loci, and testing by permutation whether these differences were greater than expected by chance. Six loci were nominally associated with stripedness at p < 0.05, but not after FDR correction. Moreover, the aggregate allele frequency difference was not greater than expected by chance. I then asked whether the top 10 (arbitrary but perhaps justifiable cut-off given the expected simple architecture of the stripe phenotype and the human obsession with multiples of 10) trait-associated SNPs (i.e., SNPs with the largest allele frequency difference between striped and non-striped individuals) had greater evidence of being affected by selection (in aggregate) than expected by chance (permutation test). The answer is that they do, with p = 0.0012 (note three have P(exceptional) > 0.9 which has an overall frequency of about 1%). Thus, I think these results indicate that stripe-associated SNPs were affected by selection, but were not among the SNPs with the greatest evidence of (most affected by) selection. Here is the relevant R code for the analysis (see projects/timema_mortatlity_experiment/).
## analyses to test whether exceptional change snps are associated with color pattern in the transplant experiment survivors
G<-matrix(scan("mngenTimemaAll.txt",n=186575*500,sep=" "),nrow=186575,ncol=500,byrow=TRUE)
nloc<-186575
nind<-500
P<-read.table("transIdStripe.txt",header=TRUE)
StrtC<-read.table("genomeSelectionTrtC.txt",header=TRUE)
Xsnps<-which(StrtC$pExceptional > 0.975) ## list of 114 exceptional change snps
## subset G and P by survivor and exceptional
Gsub<-G[Xsnps,which(P$surve == 1)]
Psub<-P$stripe[which(P$surve == 1)]
## af difference between stripe and non-stripe
nsurv<-140
nsel<-114
af<-matrix(NA,nrow=nsel,ncol=2)
for(i in 1:nsel){
af[i,]<-tapply(X=Gsub[i,],INDEX=Psub,mean)/2
}
difaf<-abs(af[,1]-af[,2])
## null distribution
nrep<-10000
null<-matrix(NA,nrow=nsel,ncol=nrep)
for(j in 1:nrep){
ranPsub<-sample(Psub,length(Psub),replace=FALSE)
for(i in 1:nsel){
out<-tapply(X=Gsub[i,],INDEX=ranPsub,mean)/2
null[i,j]<-abs(out[1]-out[2])
}
}
pvalue<-numeric(nsel)
for(i in 1:nsel){
pvalue[i]<-mean(difaf[i] < null[i,])
}
sort(p.adjust(pvalue,"fdr"))
## six are nominally significan at p < 0.05, but none are after fdr correction
## all loci,
## af difference between stripe and non-stripe
nsurv<-140
afall<-matrix(NA,nrow=nloc,ncol=2)
for(i in 1:nloc){
afall[i,]<-tapply(X=Gsurv[i,],INDEX=Psub,mean)/2
}
difafall<-abs(afall[,1]-afall[,2])
## top 10 loci
A<-which(difafall > 0.157)
mnPexc<-mean(StrtC$pExceptional[A])
nrep<-10000
nullPexc<-numeric(nrep)
for(i in 1:nrep){
nullPexc[i]<-mean(sample(StrtC$pExceptional,10,replace=FALSE))
}
mean(mnPexc < nullPexc)
##p = 0.0012