## 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)