我已经通过Gemma得到了关联分析的结果，如下。
prefix.log.txt 中包含了一个总的PVE，这不是我们想要的。
那么，如何计算这些位点的表型解释率？
据了解，有些关联分析软件是可以同时得到这个信息的，比如Tassel。
参考：Whole-genome resequencing of wild and domestic sheep identifies genes associated with
morphological and agronomic traits
有人说GAPIT的结果有这个信息。
我们知道PVE=R^2，在GAPIT结果中确实有一列是SNP的R方。但从值来看，应该不是PVE。
官方没有具体解释：
有人回答如下计算方法，但同时有人反对：
如果是GEMMA出来的结果，用上面这个公式是比较方便的。唯一不确定的是gemma中的af不是maf，不过从公式来看，不管是maf还是1-maf，结果不影响。
于是，我用了一下：
get_pve <- function(af,beta,se,N=217){
MAF=af
# MAF=1-af
PVE = (2*(beta^2)*MAF*(1-MAF))/(2*(beta^2)*MAF*(1-MAF)+((se^2)*2*N*MAF*(1-MAF)))
return(PVE)
}
结果有点偏大，值得商榷。
另外，我在一篇博文中，看到了类似GAPIT代码来计算PVE的。
https://aozhangchina.github.io/R/PVE/PVE.html
试了下，不好用。首先必须是在windows下（调用时弹框选择文件），其次要求hmp.txt文件，但是这个文件必须是单等位基因的。说实话，我没有耐心去改脚本。不过仍然感谢作者分享。
和几位网友交流，鉴于他们都是做人类疾病的，提供了几个计算方法。
-
一是孟德尔随机化书中的公式，这个比较准确。
image.png -
R包Twosamplemr
https://mrcieu.github.io/TwoSampleMR/articles/introduction.html
image.png R包gtx的grs.summary
https://www.rdocumentation.org/packages/gtx/versions/0.0.8/topics/grs.summary
人类做的很细致，这些方法在动植物研究中少见。不知可行否？
为更加了解PVE，可参考：全基因组关联分析项目设计——标记对表型的解释率
