笔记 GWAS 操作流程5-1：根红苗正的GWAS分析软件：GEMMA

笔记 GWAS 操作流程5-1：根红苗正的GWAS分析软件：GEMMA

1. GEMMA软件介绍

这个肯定厉害了，是大家闺秀，是名门望族，是根红苗正的GWAS分析软件。

GEMMA名称来源：

• G： Genome-wide
• E：Efficient
• MM：Mixed-model
• A：Association

GEMMAX主要特点：

快，话说同样的检测方法，GEMMA跑了3.3小时，而EMMA估计要跑27年？？？

2. GEMMA语法特点

相对于plink的语法，GEMMA语法更简练，一个杠，一个字母。比如：

• 表型数据：-p
• 协变量：-c
  而plink的语法的是两个杠，一个单词：
• 表型数据：--pheno
• 协变量：--covar

GEMMA支持plink的二进制文件：

• 读取plink文件：-bfile

表型数据格式：
一列，注意顺序和基因组的ID顺序一致，如果是多个性状，那就是多列，没有ID列。

协变量格式：
协变量是数字类型，如果有因子，需要转化为虚拟变量的形式。没有ID列，第一列是cov1，第二列是cov2……，注意，如果有一协变量，第一列为1（截距），第二列为协变量。

比如下面的数据是一个协变量，第一列为截距。

3. GEMMA分析一般线性模型没有协变量

首先将plink格式转化为二进制的plink格式：

plink --file b --make-bed --out c

然后将表型数据提取单独一列：

awk '{print $3}' phe.txt >p.txt

然后进行一般线性模型关联分析：

gemma-0.98.1-linux-static -bfile c -p p.txt -lm 1

结果和plink的linear结果对比：

plink的结果：

gemma的结果：

两者结果完全一致。

事实上，加上协变量的分析，gemma和plink的结果也是一样的，因为都是应用的是一般线性模型。

4. GEMMA分析混合线性模型

第一步：先生成G矩阵

gemma-0.98.1-linux-static -bfile c -gk 2 -p p.txt 

代码解释：

• -bfile 读取plink的二进制文件
• -gk 2 标准化的方法计算G矩阵
• -p 读取表型数据（这个不能省掉）

GEMMA 0.98.1 (2018-12-10) by Xiang Zhou and team (C) 2012-2018
Reading Files ... 
## number of total individuals = 1500
## number of analyzed individuals = 1500
## number of covariates = 1
## number of phenotypes = 1
## number of total SNPs/var        =    10000
## number of analyzed SNPs         =     3946
Calculating Relatedness Matrix ... 
================================================== 100%
**** INFO: Done.

G矩阵在output文件夹下：result.sXX.txt

第二步：使用混合线性模型进行GWAS分析

gemma-0.98.1-linux-static -bfile c -k output/result.sXX.txt -lmm 1 -p p.txt 
代码解释：
* -bfile plink的二进制文件
* -k 读取G矩阵的文件
* -lmm 1 使用Wald的方法进行SNP检验
* -p 表型数据

GEMMA 0.98.1 (2018-12-10) by Xiang Zhou and team (C) 2012-2018
Reading Files ... 
## number of total individuals = 1500
## number of analyzed individuals = 1500
## number of covariates = 1
## number of phenotypes = 1
## number of total SNPs/var        =    10000
## number of analyzed SNPs         =     3946
Start Eigen-Decomposition...
pve estimate =0.120599
se(pve) =0.0279188
================================================== 100%
**** INFO: Done.

第三步：查看结果文件
结果在output文件夹下：result.assoc.txt

5. GEMMA中LM模型和LMM模型的结果比较

setwd("/home/dengfei/gwas/qmsim/dat/plink_file/10_gemma_analysis_lmm/output")

mm_re = fread("result.assoc.txt")
head(mm_re)

lm_re = fread("../../06_gemma_analysis_lm/output/result.assoc.txt")
head(lm_re)
lm_re1 = lm_re[!is.na(p_wald)]
head(lm_re1)

dim(mm_re)
dim(lm_re1)

head(mm_re)
head(lm_re1)

re1 = merge(mm_re,lm_re1,by="rs")
head(re1)

# Pvalue 比较
cor(re1$p_wald.x,re1$p_wald.y)
plot(re1$p_wald.x,re1$p_wald.y)

# Beta回归系数比较
cor(re1$beta.x,re1$beta.y)
plot(re1$beta.x,re1$beta.y)

Pvalue比较

> cor(re1$p_wald.x,re1$p_wald.y)
[1] 0.5278895

Beta回归系数比较：

> cor(re1$beta.x,re1$beta.y)
[1] 0.8357564

注意：
这里面，混合线性模型分析时，一般要加上PCA的结果，防止群体结构造成的假阳性。下一章节，我们介绍LMM模型，如何加入协变量。