笔记 | GWAS 操作流程5-1：根红苗正的GWAS分析软件：GEMMA
笔记 | GWAS 操作流程5-1：根红苗正的GWAS分析软件：GEMMA
邓飞
发布于 2020-06-04 17:51:00
发布于 2020-06-04 17:51:00
笔记 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模型，如何加入协变量。
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2020-06-01，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录
腾讯云开发者
