【原】根红苗正的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