【GWAS】如何计算显著关联位点的表型解释率PVE（phenotypic variation explained）？

我已经通过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文件，但是这个文件必须是单等位基因的。说实话，我没有耐心去改脚本。不过仍然感谢作者分享。

---

---

和几位网友交流，鉴于他们都是做人类疾病的，提供了几个计算方法。

• 一是孟德尔随机化书中的公式，这个比较准确。
  

• R包Twosamplemr
  https://mrcieu.github.io/TwoSampleMR/articles/introduction.html
  

• R包gtx的grs.summary
  https://www.rdocumentation.org/packages/gtx/versions/0.0.8/topics/grs.summary

人类做的很细致，这些方法在动植物研究中少见。不知可行否？

为更加了解PVE，可参考：全基因组关联分析项目设计——标记对表型的解释率