Post date: Jul 02, 2018 2:21:29 AM
I generated some initial summaries of the gemma results for the NSF CAREER grant (more to do later). The R script (commands.R) is in the Gemma output subdirectory. The short of it is that heritabilities are quite high (with high PGE) for relative proportions of short and long chain CHCs but less so for intermediate length and degree of methylation. All interesting. Estimates of genetic architecture parameters (medians and 95% ETPIs). There are even individual SNPs with high PIPs and big model-averaged effects (see bit of code at the end for 26 carbon CHCs).
## params: 1=h, 2=pve, 3=rho, 4=pge, 5=pi, 6=n_gamma
J<-6
I<-500000
N<-11
hyp<-vector("list",N)
for(n in 1:N){
x<-paste("comb_trait",n,"_hyp.txt",sep="")
hyp[[n]]<-matrix(scan(x,n=J*I),nrow=I,ncol=J,byrow=TRUE)
}
## pve
pve<-matrix(NA,nrow=N,ncol=3)
for(n in 1:N){
pve[n,]<-quantile(hyp[[n]][,2],probs=c(0.5,0.025,0.975))
}
## pge
pge<-matrix(NA,nrow=N,ncol=3)
for(n in 1:N){
pge[n,]<-quantile(hyp[[n]][,4],probs=c(0.5,0.025,0.975))
}
## pve x pge
pveXpge<-matrix(NA,nrow=N,ncol=3)
for(n in 1:N){
pveXpge[n,]<-quantile(hyp[[n]][,4]*hyp[[n]][,2],probs=c(0.5,0.025,0.975))
}
## n-gamma
ng<-matrix(NA,nrow=N,ncol=3)
for(n in 1:N){
ng[n,]<-quantile(hyp[[n]][,6],probs=c(0.5,0.025,0.975))
}
nms<-c("C23","C25","C26","C27","C28","C29","C30","C31","M0","M1","M2")
pdf("chcHeritability.pdf",width=5,height=4)
par(mar=c(5,5,.5,.5))
barplot(pve[,1],names.arg=nms,ylab="Heritability",cex.lab=1.5,cex.axis=1.1,cex.names=1.1,xlab="CHC trait",ylim=c(0,1),las=2)
barplot(pveXpge[,1],add=TRUE,col="forestgreen",axes=FALSE)
legend(8.2,0.95,c("PVE","PVE X PGE"),fill=c("gray","forestgreen"))
dev.off()
rownames(pve)<-nms;colnames(pve)<-c("q50","q2.5","q97.5")
rownames(pge)<-nms;colnames(pge)<-c("q50","q2.5","q97.5")
rownames(ng)<-nms;colnames(ng)<-c("q50","q2.5","q97.5")
pve
# q50 q2.5 q97.5
#C23 0.6234156 0.41340893 0.8843749
#C25 0.7350265 0.26567905 0.9586744
#C26 0.8148146 0.67804200 0.9370611
#C27 0.4382373 0.12512326 0.7773366
#C28 0.3297600 0.09951869 0.6670562
#C29 0.3619402 0.04173193 0.8227529
#C30 0.3860645 0.18372040 0.6084682
#C31 0.5589508 0.35025748 0.7898377
#M0 0.2742083 0.01673294 0.6997185
#M1 0.3843119 0.03304737 0.8780734
#M2 0.1864618 0.01037629 0.5179232
pge
# q50 q2.5 q97.5
#C23 0.8250224 0.5624344 0.9898174
#C25 0.8226940 0.1650731 0.9925112
#C26 0.9653783 0.8438643 0.9985605
#C27 0.3009565 0.0000000 0.9514454
#C28 0.7602067 0.1094226 0.9886698
#C29 0.2541225 0.0000000 0.9394889
#C30 0.7937598 0.3857757 0.9898131
#C31 0.8518498 0.5903450 0.9927070
#M0 0.4188029 0.0000000 0.9667477
#M1 0.3812706 0.0000000 0.9627453
#M2 0.4004707 0.0000000 0.9663385
ng
# q50 q2.5 q97.5
#C23 5 2 18
#C25 15 3 80
#C26 11 5 22
#C27 31 0 196
#C28 3 1 61
#C29 25 0 279
#C30 3 2 11
#C31 4 2 12
#M0 13 0 149
#M1 13 0 127
#M2 9 0 258
## is this for real??
c26<-vector("list",5)
for(j in 1:5){
jj<-j-1
x<-paste("out_Trait3_Ch",jj,".param.txt",sep="")
c26[[j]]<-read.table(x,header=TRUE)
}
c26est<-matrix(0,nrow=41791,ncol=3)
for(j in 1:5){
c26est<-c26est+c26[[j]][,5:7]
}
c26est<-c26est/5