根红苗正的GWAS软件：GEMMA
现在GWAS更多使用LMM模型，这个模型plink没法做，下面介绍GEMMA软件。学习plink软件做GWAS，更多的是学习数据质控和GWAS原理，真正应用广泛的还要是混合线性模型LMM或MLM，GEMMA是一个明星软件，当然也有其它软件，比如GAPIT、FamCPU、rMVP、GCTA，最近又新出了一个fastGWS软件，后面的教程都要包含在内。
GEMMA名称来源：
- G：Genome-wide
- E：Efficient
- MM：Mixed-model
- A：Association
GEMMAX主要特点：快
就是它跑3.3h，其它软件跑27年的那种！！！
GEMMA语法特点
相对于plink的语法，GEMMA语法更简练，一个杠，一个字母。比如：
- 表型数据：-p
- 协变量：-c，而plink的语法的是两个杠接一个单词，比如表型数据：--pheno；协变量：--covar
GEMMA支持plink的二进制文件：
- 读取plink文件：-bfile
GEMMA生成G矩阵：
gemma-0.98.1-linux-static -bfile c -gk 2 -p p.txt GEMMA分析MLM模型：
gemma-0.98.1-linux-static -bfile c -k output/result.sXX.txt -lmm 1 -p p.txt结果查看：
干货来了：
GEMMA软件默认是没有PVE的（snp解释表型变异百分比），根据PVE的公式：
分子分母都有MAF*(1-MAF)，可以删除，剩下的公式为：
PVE = 2* beta^2 /( 2*beta^2 + 2*N*se^2)
这里写了一个代码脚本，输入结果文件和样本个数，就可以生成pve的结果和显著性位点。
代码如下：
$ cat ~/bin/add_pve_from_gemma_result_and_tiqu_sig.R
#! /home/dengfei/bin/Rscript
args = commandArgs(T)
if(length(args) == 0){
cat("\n\n\tRscript add_pve_from_gemma_result_and_tiqu_sig.R result.assoc.txt N，geaam的gwas结果， N为分析的样本数\n
结果生成文件：添加pve的文件：result.assoc_add_pve.txt，还有显著性文件：result_signals.csv\n")
quit("no")
}
library(data.table)
library(tidyverse)
N = as.numeric(args[2])
# d1 = fread("output/result.assoc.txt")
dd = fread(args[1])
head(dd)
dd$pve = 2*dd$beta^2/(2*dd$beta^2+dd$se^2 * 2 * (N-dd$n_miss))
head(dd)
fwrite(dd,"result.assoc_add_pve.txt",sep=" ",quote=F)
# 导出显著性的SNP位点
dat1 = dd %>% dplyr::select(SNP = rs, CHROM = chr, POS = ps,REF=allele1, ALT=allele0, maf = af,Effect = beta, SE=se, pvalue = p_wald,pve) %>% drop_na(pvalue)
head(dat1)
p_th = 1/nrow(dat1)
dat2 = dat1 %>% filter(pvalue <p_th)
fwrite(dat2,"result_signals.csv")
代码调用说明：
$ add_pve_from_gemma_result_and_tiqu_sig.R
Rscript add_pve_from_gemma_result_and_tiqu_sig.R result.assoc.txt N，geaam的gwas结果， N为分析的样本数
结果生成文件：添加pve的文件：result.assoc_add_pve.txt，还有显著性文件：result_signals.csv
代码调用演示：
结果文件：
result.assoc_add_pve.txt
result_signals.csv本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2024-03-19，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读