GWAS分析中SNP解释百分比PVE | 第三篇，MLM模型中如何计算PVE？
之前，想研究一下GWAS分析汇中PVE（表型方差解释百分比）的计算方法，写了两篇：
GWAS分析中SNP解释百分比PVE | 第一篇，SNP解释百分比之和为何大于1？
GWAS分析中SNP解释百分比PVE | 第二篇，GLM模型中如何计算PVE？
这里介绍第三篇：混合线性模型的GWAS，如何计算PVE。
1. R语言计算的PVE能否用于MLM模型？
昨天介绍了使用R语言计算显著SNP的表型方差解释百分比（PVE），它的步骤有三步：
- 第一步：将SNP和协变量（PCA和其它协变量）放到模型中，计算回归模型的R方（R-squared）这一步加上显著SNP
- 第二步：将协变量（PCA和其它协变量）放到模型中，计算回归模型的R方（R-squared）这一步去掉显著SNP
- 第三步：将第一步的R方减去第二步的R方，得到的值就是该SNP的表型变异解释百分比（PVE）
这种计算方法，也是可以用于MLM模型和其它GWAS模型的，它会存在过高估计的风险，但也是一种可行的方法。
代码演示：
第一步：
r$> mod_2 = lm(y ~ 1+pc1 + pc2 + pc3 ,data=dat1);summary(mod_2)Call:
lm(formula = y ~ 1 + pc1 + pc2 + pc3, data = dat1)Residuals:Min 1Q Median 3Q Max
-89.495 -19.318 -0.099 18.649 85.448Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 102.0467 0.8864 115.129 < 2e-16 ***
pc1 -111.8205 28.0295 -3.989 7.11e-05 ***
pc2 -21.6604 28.0295 -0.773 0.44
pc3 35.1350 28.0295 1.254 0.21
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 28.03 on 996 degrees of freedom
Multiple R-squared: 0.01783, Adjusted R-squared: 0.01487
F-statistic: 6.028 on 3 and 996 DF, p-value: 0.0004546第二步：
r$> mod_2 = lm(y ~ 1+pc1 + pc2 + pc3 ,data=dat1);summary(mod_2)Call:
lm(formula = y ~ 1 + pc1 + pc2 + pc3, data = dat1)Residuals:Min 1Q Median 3Q Max
-89.495 -19.318 -0.099 18.649 85.448Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) 102.0467 0.8864 115.129 < 2e-16 ***
pc1 -111.8205 28.0295 -3.989 7.11e-05 ***
pc2 -21.6604 28.0295 -0.773 0.44
pc3 35.1350 28.0295 1.254 0.21
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 28.03 on 996 degrees of freedom
Multiple R-squared: 0.01783, Adjusted R-squared: 0.01487
F-statistic: 6.028 on 3 and 996 DF, p-value: 0.0004546第三步：
r$> x1 - x2
[1] 0.03100801结果和GAPIT中的结果是一致的。
a3 = fread("04_glm_gapit/GAPIT.GLM.V3.GWAS.Results.csv") %>% arrange(P.value) # GAPIT GLM
a3$Rs = a3$Rsquare.of.Model.with.SNP - a3$Rsquare.of.Model.without.SNP
a3 = a3 %>% select(SNP,P.value,Rsquare.of.Model.without.SNP,Rsquare.of.Model.with.SNP,`FDR_Adjusted_P-values`,effect,Rs)
head(a3)
2. MLM的GWAS模型如何计算PVE？
GAPIT3.0增加了MLM模型计算PVE的方法：
下面是GAPIT论坛中，作者的回复，所以，我们可以在GAPIT软件中使用MLM模型，进行GWAS分析，并得到每个SNP的PVE值。
来源：https://groups.google.com/g/gapit-forum/c/acp8bCjZuJo Jiabo Wang
2021年9月13日 13:32:02收件人 GAPIT Forum Hi Manjeet,
In the GAPIT3, we can use “Random.model=TRUE” to calculate PVE
(phenotype variance explained) with the significant markers.
Sincerely, Jiabo
GWAS分析代码：
# 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")运行日志查看：
现在的GAPIT分析版本为：3.0
[1] “All packages are loaded already ! GAPIT.Version is 2020.10.24, GAPIT 3.0”
[1] “MLM”
结果文件查看：
这里关于rsquare_without_snp和rsquare_with_snp的解释如下：
“rsquare_with_snp” is the R-sqaure_LR of the model WITH the tested
SNP included in the model as an explanatory variable. In the context
of the models fitted with GAPIT, this is the mixed model with the
kinship matrix, PCs or any other covariates, and the SNP.
“rsquare_without_snp” is the R-sqaure_LR of the model WITHOUTthe
tested SNP included in the mdoel as an explanatory variable In the
context of the models fitted with GAPIT, this is the mixed model with
the kinship matrix, PCs or any other covariates.
SNP的PVE的计算方法是：
PVE = rsquare_with_snp - rsquare_without_snp
这里，我们对结果进行处理，得到：
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)
结果如下：
> head(a2)SNP Chromosome Position P.value maf nobs Rsquare.of.Model.without.SNP Rsquare.of.Model.with.SNP FDR_Adjusted_P-values effect PVE
1: M98663 20 73 0.00001208 0.2220 1000 0.06567 0.08382 0.3635 -7.475 0.01815
2: M73233 15 64 0.00012951 0.2255 1000 0.06567 0.07953 0.9039 -6.450 0.01385
3: M82957 17 57 0.00017835 0.1685 1000 0.06567 0.07895 0.9039 7.436 0.01328
4: M82958 17 57 0.00017835 0.1685 1000 0.06567 0.07895 0.9039 7.436 0.01328
5: M37588 8 51 0.00021549 0.2705 1000 0.06567 0.07861 0.9039 -6.367 0.01294
6: M36594 8 31 0.00026062 0.3520 1000 0.06567 0.07827 0.9039 -5.610 0.01260
可以看到，SNP名为M98663的PVE为0.01815
3. GLM和MLM模型PVE的比较
这里，就有问题了，PVE的计算方法，有GLM的和MLM的，应该用哪一个呢？
我们之前说过，使用R语言计算的PVE，会过高估计，如果在MLM模型中，我们将显著性的SNP提取，然后用R语言进行PVE分析，其实和GLM模型得到的PVE是一致的。
这里，我们进行比较一下两种PVE的计算结果：
GLM模型计算PVE结果
> # GAPIT 中的glm模型
> a1 = fread("04_glm_gapit/GAPIT.GLM.V3.GWAS.Results.csv") %>% arrange(P.value) # GAPIT GLM
> a1$PVE = a1$Rsquare.of.Model.with.SNP - a1$Rsquare.of.Model.without.SNP
> head(a1)SNP Chromosome Position P.value maf nobs Rsquare.of.Model.without.SNP Rsquare.of.Model.with.SNP FDR_Adjusted_P-values effect PVE
1: M57157 12 44 0.00000002658 0.3750 1000 0.01783 0.04884 0.0007997 7.518 0.03101
2: M57155 12 44 0.00000014999 0.3275 1000 0.01783 0.04543 0.0011567 7.242 0.02760
3: M57156 12 44 0.00000014999 0.3275 1000 0.01783 0.04543 0.0011567 7.242 0.02760
4: M98663 20 73 0.00000015377 0.2220 1000 0.01783 0.04538 0.0011567 -7.812 0.02755
5: M73233 15 64 0.00000051228 0.2255 1000 0.01783 0.04303 0.0029001 -7.526 0.02520
6: M8452 2 69 0.00000057828 0.4350 1000 0.01783 0.04279 0.0029001 -6.520 0.02496MLM模型计算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)SNP Chromosome Position P.value maf nobs Rsquare.of.Model.without.SNP Rsquare.of.Model.with.SNP FDR_Adjusted_P-values effect PVE
1: M98663 20 73 0.00001208 0.2220 1000 0.06567 0.08382 0.3635 -7.475 0.01815
2: M73233 15 64 0.00012951 0.2255 1000 0.06567 0.07953 0.9039 -6.450 0.01385
3: M82957 17 57 0.00017835 0.1685 1000 0.06567 0.07895 0.9039 7.436 0.01328
4: M82958 17 57 0.00017835 0.1685 1000 0.06567 0.07895 0.9039 7.436 0.01328
5: M37588 8 51 0.00021549 0.2705 1000 0.06567 0.07861 0.9039 -6.367 0.01294
6: M36594 8 31 0.00026062 0.3520 1000 0.06567 0.07827 0.9039 -5.610 0.01260
两者进行比较：
> re = merge(a1,a2,by="SNP")
> head(re)SNP Chromosome.x Position.x P.value.x maf.x nobs.x Rsquare.of.Model.without.SNP.x Rsquare.of.Model.with.SNP.x FDR_Adjusted_P-values.x effect.x PVE.x Chromosome.y
1: M100000 20 99 0.005662 0.1350 1000 0.01783 0.02541 0.1190 5.036 0.00758157 20
2: M10001 3 0 0.823235 0.0425 1000 0.01783 0.01788 0.9411 -0.716 0.00004923 3
3: M10003 3 0 0.490322 0.1175 1000 0.01783 0.01830 0.7910 -1.486 0.00046956 3
4: M10004 3 0 0.230065 0.3970 1000 0.01783 0.01925 0.5883 1.611 0.00142220 3
5: M10005 3 0 0.012427 0.2805 1000 0.01783 0.02402 0.1724 3.700 0.00618466 3
6: M10006 3 0 0.680577 0.3915 1000 0.01783 0.01800 0.8864 -0.547 0.00016722 3Position.y P.value.y maf.y nobs.y Rsquare.of.Model.without.SNP.y Rsquare.of.Model.with.SNP.y FDR_Adjusted_P-values.y effect.y PVE.y
1: 99 0.19093 0.1350 1000 0.06567 0.06728 0.9877 2.7332 0.001606674
2: 0 0.73907 0.0425 1000 0.06567 0.06578 0.9927 -1.1751 0.000104136
3: 0 0.34225 0.1175 1000 0.06567 0.06652 0.9927 -2.2101 0.000846925
4: 0 0.47611 0.3970 1000 0.06567 0.06615 0.9927 1.0501 0.000476678
5: 0 0.08458 0.2805 1000 0.06567 0.06847 0.9877 2.8215 0.002796028
6: 0 0.92109 0.3915 1000 0.06567 0.06568 0.9957 -0.1459 0.000009211
> re %>% select(PVE.x,PVE.y) %>% corPVE.x PVE.y
PVE.x 1.0000 0.7875
PVE.y 0.7875 1.0000
可以看到相关系数为0.78，相关性并不高。
两者PVE之和比较：
> sum(a1$PVE)
[1] 57.61
> sum(a2$PVE)
[1] 28.31
可以看到GLM计算的PVE之和为57.61，而MLM模型计算的PVE之和为28.31，很明显GLM的方法更虚高。
所以，在MLM模型的GWAS中，我们要选择MLM方法计算的PVE。
问题来了，如果不用GAPIT软件，该如何手动计算PVE值呢？
4. 其它GWAS分析软件如何计算PVE
我们知道，其它GWAS软件中是没有PVE的结果的，比如：
- GEMMA
- GCTA中的fast-GWA
下一节介绍一下如何用R语言进行演示MLM的PVE计算方法。
这应该是全网可以搜索的唯一的实现方法了。
欢迎继续关注。
GWAS分析中SNP解释百分比PVE | 第三篇，MLM模型中如何计算PVE？相关推荐
- GWAS分析中SNP解释百分比PVE | 第二篇，GLM模型中如何计算PVE？
上一篇,介绍了一下显著性的SNP,他们的解释表型变异百分比(PVE)之和,为何可能大于1. https://yijiaobani.blog.csdn.net/article/details/12209 ...
- GWAS分析中SNP解释百分比PVE | 第四篇，MLM模型中如何手动计算PVE？
系列部分: GWAS分析中SNP解释百分比PVE | 第一篇,SNP解释百分比之和为何大于1? GWAS分析中SNP解释百分比PVE | 第二篇,GLM模型中如何计算PVE? GWAS分析中SNP解释 ...
- GWAS分析中SNP解释百分比PVE | 第一篇，SNP解释百分比之和为何大于1？
关于GWAS分析中PVE的计算方法: 我查了一下,大体计算PVE的方法有三种:第一种回归分析或者方差分析的方法,计算R方(GLM模型),第二种是根据effect,se,maf计算PVE,第三种是根据L ...
- 中职计算机说课稿三篇,2020精选中职计算机说课稿3篇（15页）-原创力文档
2020 精选中职计算机说课稿 3 篇 2020 精选中职计算机说课稿篇 1 各位评委老师, 你们好! 我是来自 xxx 职业中专计算机专业的教 师 xxx ,今天我说课的题目是<获取屏幕图像& ...
- 职高学生计算机学情分析,精选中职计算机说课稿三篇
精选中职计算机说课稿三篇 中职计算机说课稿(一) 位评委老师,你们好!我是来自XXX职业中专计算机专业的教师XXX,今天我说课的题目是<电子表格基本操作>.本节课出自高等教育出版社出版的& ...
- 中职计算机说课稿三篇,精选中职计算机说课稿三篇-20210609060707.docx-原创力文档...
PAGE / NUMPAGES 精选中职计算机说课稿三篇 中职计算机说课稿(一) 位评委老师你们好!我是来自 XXX职业中专计算机专业的老师 XXX,今日我说课的题目是?<电子表格基本操作> ...
- 各大佬抨击ICML审稿太随意：LeCun三篇全没中，马毅说以后再也不投了
点击上方"视学算法",选择加"星标"或"置顶" 重磅干货,第一时间送达 梦晨 发自 凹非寺 量子位 | 公众号 QbitAI 顶会ICML结 ...
- [老老实实学WCF] 第三篇 在IIS中寄存服务
老老实实学WCF 第三篇 在IIS中寄宿服务 通过前两篇的学习,我们了解了如何搭建一个最简单的WCF通信模型,包括定义和实现服务协定.配置服务.寄宿服务.通过添加服务引用的方式配置客户端并访问服务.我 ...
- MLM模型中，是否应该按15%的比例mask？
翻译自 Should You Mask 15% in Masked Language Modeling? 摘要 MLM模型约定俗成按照15%的比例mask,主要基于两点:更多的mask比例对于学习更好 ...
最新文章
- Ubuntu 系统开机黑屏提示the root filesystem on /dev/sdb2 requires a manual fsck
- 柱状折线图2-双柱状重合堆积折线-重写图例点击事件
- VTK：图像理想高通用法实战
- Bug同样的shell脚本在win与linux系统下执行不一样
- 前端学习（3302）：类组件父组件和子组件createRef
- c语言用字符串统计一个整数中数字的个数_全国计算机等级考试二级C语言
- 按键精灵定义全局变量_按键精灵中如何定义和使用变量
- 数据绑定表达式语法(Eval,Bind区别)
- 那么多编程语言，为什么要选择C++？
- Silverlight实用窍门系列：51.Silverlight页面控件的放大缩小、Silverlight和Html控件的互相操作【附带源码实例】...
- 递归方式计算一个数的几次方
- Php处理输入法表情,php开发中手机输入法自带的表情、emoji表情、微信表情不显示问题，以及过虑emoji表情方法！...
- 人工智能和机器视觉技术学习培训设备
- Flash 芯片类型介绍
- 【计组】计算机乘法运算
- 电脑重装win10系统教程，简单易懂，不用U盘直接重装
- Microsoft visual studio关闭安全检查的几种方法（2015/2017）
- 打开我的计算机我的文档不见,win10系统我的文档不见了的设置教程
- go语言多package使用实战
- access工资明细表_利用ACCESS数据库报表功能制作工资条