PVE(phenotypic variation explained),是指表型变异解释率,即等位变异可以在多大程度上影响表型,这个SNP解释了表型方差的百分之多少。在数量性状的GWAS 研究里经常会看到。
gwas关联分析可以使用很多软件,也可以计算各种模型。接下来主要说明TASSEL和GEMMA怎么得到解释率
TASSEL是可以直接得到pve值的,直接使用数据即可
1.1 GLM模型
marker_Rsq 即R2值,即表型解释率PVE
1.2 MLM模型及CMLM模型
MarkerR2 即R2值,即表型解释率PVE
注意:标记的R2(即表型的解释率),因为做的是单标记检验,每个标记都是独立计算的,所以多个标记的R2求和可能大于1,即每个标记的R2之间没有什么关系。
SNP解释百分比之和为何大于1?
在TASSEL中GLM或者MLM模型中,是单标记扫描,之所以SNP的R2(R square)之和会大于1,因为标记间存在LD,比如一个标记关联的基因能解释20%的变异,这个位点附近有6个标记都存在LD状态,那么这6个标记的解释百分比之和就会是120%。
TASSEL是可以直接得到pve值的,而有些模型和软件(如FarmCPU,Fast-LMM,EMMAX,GEMMA等)不给出R2
这种情况下,有两种方法:
法一:根据effect,se,maf计算表型解释率PVE
法二:可以选出显著的位点,使用简单回归模型(Im)进行拟合,获得R2值作为解释率。
如何根据effect,se,maf计算表型解释率PVE,公式如下:
其中:
β为GWAS中的effect值
MAF为SNP的MAF次等位基因频率
se为GWAS中effect值的标准误(se)
N为GWAS中该SNP参与分析的个体数
以GEMMA软件得到的计算结果为例
1.lm模型
这 11 列结果分别是:染色体、SNP id、碱基对在染色体上的位置、相应 SNP 的缺失个体数量、相应 SNP 的无缺个体数量、最小等位基因、最大等位基因、等位基因频率、beta 检测、beta 标准误、Wald 检验的 P 值。
再次对照公式
公式中:
β 为 GEMMA结果中的 beta,beta就是GWAS中的effect值
se 为 GEMMA结果中的 se
N 为 GEMMA结果中的 n_obs;如果结果中没有 n_obs,只有n_miss,就用样本总数 - n_miss
MAF 为 GEMMA结果中的 af
2.lmm模型
lmm模型计算公式与lm一致,唯一差别就是没有n_obs,只有n_miss
在计算时,N 为 样本总数 — n_miss ,即做个减法
其余都是一样的。
在R语言中根据公式计算即可
#读取GEMMA软件GWAS分析结果
a<-read.table("lm.pvalue.txt",header = T, sep = '\t')
#计算pve——lm模型
a$pve = (2*(a$beta^2*a$af*(1-a$af)))/(2*a$beta*a$af*(1-a$af) + a$se^2*2*a$n_obs*a$af*(1-a$af))
###或者计算pve——lmm模型
a$pve = (2*(a$beta^2*a$af*(1-a$af)))/(2*a$beta*a$af*(1-a$af) + a$se^2*2*(421-a$n_miss)*a$af*(1-a$af))
#只提取这4列
b<-a[c('chr','ps','p_score','pve')]
#输出结果文件
write.table(b, file=args$output, sep = '\t', row.names = F, quote
参考:https://zhuanlan.zhihu.com/p/448255236
https://zhuanlan.zhihu.com/p/449712020
https://www.sohu.com/a/378139452_278730
https://blog.csdn.net/wt141643/article/details/105834917
如果觉得我的文章对您有用,请随意打赏。你的支持将鼓励我继续创作!