这个肯定厉害了,是「大家闺秀」,是「名门望族」,是「根红苗正」的GWAS分析软件。
「GEMMA名称来源:」
「GEMMAX主要特点:」
❝快,话说同样的检测方法,GEMMA跑了3.3小时,而EMMA估计要跑27年???
❞
相对于plink的语法,GEMMA语法更简练,一个杠,一个字母。比如:
-p
-c
而plink的语法的是两个杠,一个单词:--pheno
--covar
GEMMA支持plink的二进制文件:
-bfile
「表型数据格式:」一列,注意顺序和基因组的ID顺序一致,如果是多个性状,那就是多列,没有ID列。
「协变量格式:」协变量是数字类型,如果有因子,需要转化为虚拟变量的形式。没有ID列,第一列是cov1,第二列是cov2……,注意,如果有一协变量,第一列为1(截距),第二列为协变量。
比如下面的数据是一个协变量,第一列为截距。
「首先将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的结果也是一样的,因为都是应用的是一般线性模型。
「第一步:先生成G矩阵」
gemma-0.98.1-linux-static -bfile c -gk 2 -p p.txt
代码解释:
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
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模型,如何加入协变量。