Post date: Oct 31, 2019 3:39:9 AM
I want to see whether or not allowing for epistasis for color affects the form of selection, in other words, does epistasis push individuals in parts of phenotype space that change selection. To this end, I am estimating selection gradients on GEBVs estimated from the additive and additive plus mapit epistasis models from gemma.
1. I first calcualted the GEBVs. This was done in /uufs/chpc.utah.edu/common/home/u6000989/projects/timema_genarch/color_exp_2019_gbs/Gemma_mel/output with the calcBV.R script shown below:
## calculate GEBVs for RG and GB with and w/o mapit epistasis SNPs
G_add<-as.matrix(read.table("../g_tchum.txt",header=FALSE))
G_epi<-as.matrix(read.table("../g_epi_tchum.txt",header=FALSE))
for(k in 1:dim(G_add)[1]){
G_add[k,]<-G_add[k,]-mean(G_add[k,])
}
for(k in 1:dim(G_epi)[1]){
G_epi[k,]<-G_epi[k,]-mean(G_epi[k,])
}
## drop fixed SNPs
fixSnp_add<-which(apply(G_add,1,sd)==0)
fixSnp_epi<-which(apply(G_epi,1,sd)==0)
## get list of files from GEMMA
f_bv_rg<-list.files(pattern=glob2rx("o_mod_g_tchum_ph1_ch*bv.txt"))
f_bv_gb<-list.files(pattern=glob2rx("o_mod_g_tchum_ph2_ch*bv.txt"))
f_eff_rg<-list.files(pattern=glob2rx("o_mod_g_tchum_ph1_ch*param.txt"))
f_eff_gb<-list.files(pattern=glob2rx("o_mod_g_tchum_ph2_ch*param.txt"))
f_bv_epi_rg<-list.files(pattern=glob2rx("oe_epi_tchum_ph1*bv.txt"))
f_bv_epi_gb<-list.files(pattern=glob2rx("oe_epi_tchum_ph2*bv.txt"))
f_eff_epi_rg<-list.files(pattern=glob2rx("oe_epi_tchum_ph1*param.txt"))
f_eff_epi_gb<-list.files(pattern=glob2rx("oe_epi_tchum_ph2*param.txt"))
## function to calcualte GEBVs for subset of individuals not in each training set,
calcGEBV<-function(bvf=NA,efff=NA,G=NA){
N<-length(bvf)
gebv<-vector("list",N)
for(j in 1:N){
bv<-scan(bvf[j])
eff<-read.table(efff[j],header=TRUE)
L<-dim(eff)[1]
I<-length(bv)
gebv[[j]]<-rep(NA,I)
mav<-eff$alpha + eff$beta * eff$gamma
for(i in 1:437){
gebv[[j]][i]<-sum(mav*G[,i])
}
}
gebv<-matrix(unlist(gebv),nrow=I,ncol=N)
return(gebv)
}
gebv_rg_add<-apply(calcGEBV(bvf=f_bv_rg,efff=f_eff_rg,G=G_add[-fixSnp_add,]),1,mean)
gebv_gb_add<-apply(calcGEBV(bvf=f_bv_gb,efff=f_eff_gb,G=G_add[-fixSnp_add,]),1,mean)
gebv_rg_epi<-apply(calcGEBV(bvf=f_bv_epi_rg,efff=f_eff_epi_rg,G=G_epi[-fixSnp_epi,]),1,mean)
gebv_gb_epi<-apply(calcGEBV(bvf=f_bv_epi_gb,efff=f_eff_epi_gb,G=G_epi[-fixSnp_epi,]),1,mean)
write.table(file="gebv_add.txt",round(cbind(gebv_rg_add,gebv_gb_add),5),row.names=FALSE,col.names=FALSE)
write.table(file="gebv_epi.txt",round(cbind(gebv_rg_epi,gebv_gb_epi),5),row.names=FALSE,col.names=FALSE)
This generates the two gebv* files that I will use for the selection gradients.