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

我已经通过Gemma得到了关联分析的结果，如下。

image.png

prefix.log.txt 中包含了一个总的PVE，这不是我们想要的。

image.png

那么，如何计算这些位点的表型解释率？

据了解，有些关联分析软件是可以同时得到这个信息的，比如Tassel。
image.png
image.png

参考：Whole-genome resequencing of wild and domestic sheep identifies genes associated with
morphological and agronomic traits

有人说GAPIT的结果有这个信息。
image.png

我们知道PVE=R^2，在GAPIT结果中确实有一列是SNP的R方。但从值来看，应该不是PVE。

image.png
官方没有具体解释：
image.png

有人回答如下计算方法，但同时有人反对：
image.png

如果是GEMMA出来的结果，用上面这个公式是比较方便的。唯一不确定的是gemma中的af不是maf，不过从公式来看，不管是maf还是1-maf，结果不影响。
image.png

于是，我用了一下：

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

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

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

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

posted @ 2021-07-29 15:29  生物信息与育种  阅读(2948)  评论(0)  编辑  收藏  举报