育种数据分析之放飞自我分享 http://blog.sciencenet.cn/u/yijiaobai 关注:生物统计,数量遗传,混合线性模型,生物信息,R,Perl,Python,GWAS,GS相关方法,文章及代码

博文

GWAS分析中SNP解释百分比PVE | 第四篇,MLM模型中如何手动计算PVE?

已有 1434 次阅读 2021-12-25 12:25 |个人分类:GWAS|系统分类:科研笔记

「系列部分:」

GWAS分析中SNP解释百分比PVE | 第一篇,SNP解释百分比之和为何大于1?

GWAS分析中SNP解释百分比PVE | 第二篇,GLM模型中如何计算PVE?

GWAS分析中SNP解释百分比PVE | 第三篇,MLM模型中如何计算PVE?

今天介绍一下如何手动计算MLM模型GWAS的PVE结果。因为GAPIT中的MLM模型又PVE结果,但是常用的GEMMA、GCTA的GWAS结果并没有PVE,本篇介绍一下如何根据GWAS结果手动计算,用R语言进行演示。

1. 参考文献

首先是这个论坛的内容:How to determine the percent phenotypic variation explained (PVE) by a selected SNP?

Feofanova, Elena. (2020). Re: How to determine the percent phenotypic variation explained (PVE) by a selected SNP?. Retrieved from: https://www.researchgate.net/post/How_to_determine_the_percent_phenotypic_variation_explained_PVE_by_a_selected_SNP/5e1600dd661123743209bc17/citation/download.

里面介绍了计算方法:

其中参考的文献是:

Shim, H., Chasman, D.I., Smith, J.D., Mora, S., Ridker, P.M., Nickerson, D.A., Krauss, R.M., and Stephens, M. (2015). A multivariate genome-wide association analysis of 10 LDL subfractions, and their response to statin treatment, in 1868 Caucasians. PLoS One 10, e0120758.

里面的附件1:Supplementary Information: S1: Computing proportion of variance in phenotype explained by a given SNP (PVE).

整个公式如下:

最后,完整的公式如下:其中:

  • 为GWAS中的effect值
  • MAF 为SNP的MAF次等位基因频率
  • 为GWAS中effect值的标准误(se)
  • N 为GWAS中该SNP参与分析的个体数

2. GAPIT中MLM模型分析PVE值

gaipit中的MLM模型代码如下:

# GWAS 分析
library(data.table)
source("http://zzlab.net/GAPIT/GAPIT.library.R")
source("http://zzlab.net/GAPIT/gapit_functions.txt")

myGd = fread("mdp_numeric.txt",header=T,data.table = F)
myGM = fread("mdp_SNP_information.txt",header = T,data.table=F)

myY = fread("dat_plink.txt",data.table = F)
head(myY)

covar = fread("cov_plink.txt",data.table = F)[,-1]
names(covar)[1] = "Taxa"
head(covar)

myGAPIT = GAPIT(Y = myY[,c(1,3)],GD = myGd, GM = myGM,
                # PCA.total=3,
                CV = covar,Random.model=TRUE,
                model="MLM")

对结果进行处理,计算PVE值,结果如下:

a2 = fread("05_lmm_gapit/GAPIT.MLM.V3.GWAS.Results.csv") %>% arrange(P.value) # GAPIT MLM
a2$PVE = a2$Rsquare.of.Model.with.SNP - a2$Rsquare.of.Model.without.SNP
head(a2)
a2 = a2 %>% select(1,lmm_effect = effect,maf = maf,lmm_p.value = P.value,
                   lmm_FDR_p = `FDR_Adjusted_P-values`,
                   lmm_Rsquare_with_snp = Rsquare.of.Model.with.SNP,
                   lmm_Rsquare_without_snp = Rsquare.of.Model.without.SNP,
                   lmm_Rs = Rs)
head(a2)

这里,GWAS结果中,因为没有effect的se的值,所以无法手动运算,下面我们看一下GEMMA和GCTA的fast-GWA,用同样的数据,进行GWAS分析,并手动计算PVE值,和GAPIT中的MLM模型的PVE值进行对比。

3. GEMMA进行MLM模型的GWAS分析

GEMMA进行GWAS分析,分为两步:

  • 第一步:构建G矩阵
  • 第二部:进行MLM的GWAS分析
# 构建G矩阵
gemma-0.98.1-linux-static -bfile ../geno/b -gk 2 -p p1.txt

# 进行LMM分析
gemma-0.98.1-linux-static -bfile ../geno/b -k output/result.sXX.txt -lmm 1 -p p1.txt  -c cov.txt

结果如下:结果中:

  • beta为effect
  • se为se
  • p_wald 为P值
  • n_miss 为总个体数的缺失,n为总个体数减去缺失
  • af为maf次等位基因频率

所以上面结果,读到R语言中,用下面公式进行计算PVE:这里的N为1000,计算结果如下:

a4$pve = (2*(a4$beta^2*a4$af*(1-a4$af)))/(2*a4$beta*a4$af*(1-a4$af) + a4$se^2*2*1000*a4$af*(1-a4$af))
head(a4)

比较GAPIT的MLM模型的PVE和手动根据GEMMA的MLM计算的PVE结果,可以看到SNPM98663,它的PVE在GAPIT中是0.01815,在GEMMA中是0.01988,结果有些差异,下面我们看一下相关系数。

「比较相关系数:」

re = merge(a2,a4,by.x="SNP",by.y = "rs")
head(re)
re %>% select(P.value,p_wald) %>% cor
re %>% select(effect,beta) %>% cor
re %>% select(PVE,pve) %>% cor
re %>% select(PVE,pve) %>% plot

GAPIT和GEMMA的P值比较结果:0.9986,基本一致。

GAPIT和GEMMA的effect值比较结果:0.9996,基本一致。

GAPIT和GEMMA的PVE值比较结果:0.9991,基本一致。

两款软件的PVE的散点图:可以看到,上面的手动计算方法,和GAPIT的MLM模型的PVE结果完全一致。

4. GCTA进行MLM模型的GWAS分析

GCTA进行GWAS分析,分为两步:

  • 第一步:构建G矩阵
  • 第二步:生产稀疏矩阵
  • 第三步:进行MLM的GWAS分析
gcta --bfile ../geno/b --make-grm --out geno_grm --make-grm-alg 1
gcta --grm geno_grm --make-bK-sparse 0.05 --out sp_grm
gcta --bfile ../geno/b --grm-sparse sp_grm --fastGWA-mlm --pheno dat_plink.txt --qcovar cov_plink.txt --out a1

生成的结果如下:

下面,我们读入到R语言中,手动计算PVE值。

a5 = fread("10_fast_GWA_MLM/a1.fastGWA") %>% arrange(P)
head(a5)

a5$pve = (2*(a5$BETA^2*a5$AF1*(1-a5$AF1)))/
    (2*a5$BETA*a5$AF1*(1-a5$AF1) + a5$SE^2*2*1000*a5$AF1*(1-a5$AF1))
head(a5)

和GAPIT的MLM模型比较PVE结果:

re = merge(a2,a5,by.x="SNP",by.y = "SNP")
head(re)

re %>% select(P.value,P) %>% cor
re %>% select(effect,BETA) %>% cor
re %>% select(PVE,pve) %>% cor
re %>% select(PVE,pve) %>% plot

结果如下:

P值相关系数为0.96,effect相关系数为0.98,PVE相关系数为0.98,基本一致。

5. GEMMA和GCTA手动计算PVE结果可行

所以,经过上面的测试,我们可以得到结论:

  • 对于GEMMA和GCTA软件,计算的GWAS结果,可以根据公式计算PVE
  • 结果和GAPIT结果一致

所以,网站上面各种搜索GEMMA如何计算PVE,GCTA如何计算PVE,EMMA如何计算PVE的各种问题,可以休矣。

6. 讨论

读到此,你是否有一种豁然开朗的感觉,GWAS分析中显著SNP如何计算解释百分比(PVE)的相关问题,终于解决了。

当然,有些GWAS方法是没有给出se的,比如FAMcpu等,那就不能用这种方法进行手动计算了。

需要注意的是,PVE的方法,之和远远大于1,是因为显著SNP之间存在LD,因为SNP代表的是基因,如果存在LD较高,那就是基因被代表了很多次,所以PVE就会偏高,我们不能说8个SNP解释了表型60%的变异,因为这8个SNP可能是连锁的,他们之和被高估了。

另外,从理论上来说,PVE的上限是遗传力(h2),比如GEMMA的结果中:给出的PVE是所有SNP的PVE之和,从算法上来说,就是Va/(Va+Ve),就是遗传力。所以这里,给出所有PVE之和的上限就是遗传力。

所以,在描述结果是,如果你的性状遗传力为0.3,那就表示你所有的SNP的解释百分比之和理论上限是30%,如果你计算的10个显著性的SNP的PVE之和为40%,然后还说自己的SNP多么牛叉,多么重要,这明显是不合适的,里面有很大重复估计的PVE在里面。

当然,相对于GLM的PVE计算(也就是R语言的单标记回归计算R-squared),MLM的计算方法重复估计偏低一点。之前的博客中有比较,同样的数据,GLM的PVE之和为50,而MLM的PVE之和为25。

最后,如果想要更严谨的计算多个SNP的解释百分比,或者一个区段内显著SNP的解释百分比(PVE),可以将该区段作为随机因子,在LMM模型中估算其方差组分,然后计算Vsnp/Vtotal的比值,这应该会降低假阳性,是更严谨的方法。具体文献见:

Citation: Tang Z, Xu J, Yin L, Yin D, Zhu M, Yu M, Li X, Zhao S and Liu X (2019) Genome-Wide Association Study Reveals Candidate Genes for Growth Relevant Traits in Pigs. Front. Genet. 10:302. doi: 10.3389/fgene.2019.00302

里面将显著的SNP区段作为block,进行方差组分的估计,进而计算PVE:之前,在星球内,有朋友问我如何计算PVE,我当时给了三个方法:

第一种:是使用R语言的回归分析去做,这个也是GLM的GWAS计算PVE的方法 第二种:是根据effect、se,maf去计算,这个也是MLM的GWAS计算PVE的方法 第三种:是将显著的区段(block)放到LMM模型中,计算PVE,这个就是上面文献计算的方法。

「2021年圣诞节的周六,花了一上午进行了PVE不同方法的整理,相信会对有相关疑惑的人有所帮助,我浏览了各种论坛,都没有找到一个确切的方案。这里,我用实际数据进行了测试,总结了几种方法,所以,你看到的应该是互联网上第一篇使用的方案,赶紧点个赞吧!」

欢迎关注我的公众号:育种数据分析之放飞自我。主要分享R语言,Python,育种数据分析,生物统计,数量遗传学,混合线性模型,GWAS和GS相关的知识。

下一篇,将相关的资料,整理为一个pdf,将相关论文,资料,代码,数据进行一个整理,方便需要的朋友进行重演和测试,欢迎继续关注。




https://wap.sciencenet.cn/blog-2577109-1318085.html

上一篇:GWAS分析中SNP解释百分比PVE | 第三篇,MLM模型中如何计算PVE?
下一篇:R语言相关性分析和相关性分析可视化常用方法汇总

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2022-6-29 08:43

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部