Post date: Dec 17, 2013 12:9:56 AM
I would like to know whether the expected phenotype (survival on alfalfa) is higher for alfalfa feeding populations than non-alfalfa feeding populations. I began by looking just at SLA and GLA, which were the populations in the experiment. This is a tricky question. First, we have estimates of genetic architecture of survival on alfalfa in GLA, SLA, and GLA and SLA combined. The latter is not on the binary scale, but uses corrected phenotypes (this could be weird), whereas it is not clear which population specific architecture is most appropriate to use. It makes biological sense to consider the architecture from each population for itself, but the effects are with respect to a latent risk scale (normal quantiles) that is forced to be the same across populations. Thus, the expected value for each population given its own architecture should basically always be the same. On the other hand it is not the case that the same variants were associated with survival in each population, so assuming the same architecture is not appropriate.
The only reasonable, albeit not perfect, approach is to use the parameters from the combined architecture estimates. This is what I did (though I played around with many alternatives). The estimated breeding values are quite similar across populations and from different host plant treatments. Thus, we find no evidence that the genetic effects we estimated are associated with fitness in the wild. Here is the R code I used:
## read genotype data
g_gla<-read.table("../allelefreq/mod_genloc_gla.txt",header=FALSE)
g_sla<-read.table("../allelefreq/mod_genloc_sla.txt",header=FALSE)
## read wild binary vectors
wild_gla<-read.table("../allelefreq/glaMsWild.txt",header=FALSE)
wild_sla<-read.table("../allelefreq/slaMsWild.txt",header=FALSE)
## create cenetered genoypte matrixes, centered using mean experiment genotype on Ms
G_gla<-g_gla[,-c(1:2)]
G_sla<-g_sla[,-c(1:2)]
Gc_glaxgla<-G_gla[,wild_gla==1]-apply(G_gla[,wild_gla==0],1,mean)
Gc_glaxsla<-G_gla[,wild_gla==1]-apply(G_sla[,wild_sla==0],1,mean)
Gc_slaxgla<-G_sla[,wild_sla==1]-apply(G_gla[,wild_gla==0],1,mean)
Gc_slaxsla<-G_sla[,wild_sla==1]-apply(G_sla[,wild_sla==0],1,mean)
## read sparse effect estimates for individual populations
beta<-read.table("results/effectsSparseSurvProbitBeta.txt",header=TRUE)
## breeding values
bv_gla_archGlaMsCenGla<-apply(Gc_glaxgla * beta$glaMs,2,sum,na.rm=TRUE)
bv_sla_archGlaMsCenGla<-apply(Gc_slaxgla * beta$glaMs,2,sum,na.rm=TRUE)
bv_gla_archSlaAcCenSla<-apply(Gc_glaxsla * beta$slaAc,2,sum,na.rm=TRUE)
bv_sla_archSlaAcCenSla<-apply(Gc_slaxsla * beta$slaAc,2,sum,na.rm=TRUE)
bv_sla_archSlaMsCenSla<-apply(Gc_slaxsla * beta$slaMs,2,sum,na.rm=TRUE)
## read total effect from the combined analysis
betaC<-read.table("results/effectsCombSurvTotal.txt",header=TRUE)
## standardize based on combined analysis
G_gla<-g_gla[,-c(1:2)]
G_sla<-g_sla[,-c(1:2)]
Gc_glaxcomb<-G_gla[,wild_gla==1]-apply(cbind(G_gla[,wild_gla==0],G_sla[,wild_sla==0]),1,mean)
Gc_slaxcomb<-G_sla[,wild_sla==1]-apply(cbind(G_gla[,wild_gla==0],G_sla[,wild_sla==0]),1,mean)
## breeding values
bv_glaMx_comb<-apply(Gc_glaxcomb * betaC$Ms,2,sum,na.rm=TRUE)
bv_slaMs_comb<-apply(Gc_slaxcomb * betaC$Ms,2,sum,na.rm=TRUE)
bv_glaAc_comb<-apply(Gc_glaxcomb * betaC$Ac,2,sum,na.rm=TRUE)
bv_slaAc_comb<-apply(Gc_slaxcomb * betaC$Ac,2,sum,na.rm=TRUE)