【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，可參考：全基因組關聯分析項目設計——標記對表型的解釋率