GEMMA(Genome-wide Efficient Mixed Model Association algorithm)是一款基于混合线性模型的GWAS分析软件。相比较于其他基于混合线性模型的软件,它有如下优势:1.快速:远远快于其他精确算法(EMMA和FaST-LMM)。2.准确:EMMAX和GAPIT都采用固定零模型中的方差组分不变的策略来提高运算速度,其实就是一种近似算法,不如GEMMA准确。3.方便:可直接使用plink二进制格式数据,无需进行复杂的数据格式转换。4.功能全面:可进行单标记GWAS、多标记GWAS和多性状GWAS分析。
Step0. 安装Linux子系统:
虽然GEMMA是在Linux系统下运行的程序,使用Windows系统的同学也不用担心。Win10可以很方便地安装Ubuntu子系统!
Step1. 安装GEMMA:
推荐使用conda安装软件,不用下载源代码编译,简单快捷。
Step2. Win下使用plink1.90进行基因型过滤
plink ped文件要求前六列由[家系id],[个体id],[父亲id],[母亲id],[性别] 和 [表型]组成,从第七列开始才是基因分型结果。一般我们拿到的plink数据只有个体id号和分型数据,就需要在过滤的同时告知plink我们只有个体id号,没有其他5列。具体过滤代码如下:
plink --file rawfile --no-fid --no-parents --no-pheno --no-sex --mind 0.1 --maf 0.01 --geno 0.1 --hwe 0.000001 --make-bed --out qc_data
参数注解:
--file 原始文件名 --no-fid 没有家庭id --no-parents没有父母id号 --no-pheno没有表型 --no-sex 没有性别
--mind 0.1 删除基因型检出率<90%的个体
--maf 0.01 删除最小等位基因频率<1%的位点
--geno 0.1 删除基因型检出率<90%的位点
--hwe 0.000001 删除哈代-温伯格平衡检验<0.000001的位点
--make-bed 转换文件为二进制格式
--out 输出文件名
Step3. 运行GEMMA:
运行GEMMA需要以下文件:
1. plink二进制文件
GEMMA要求在plink二进制文件中需包含表型数据,fam文件的第六列就是表型数据。因为我们的二进制文件中没有表型,因此fam文件第6列为-9。使用GEMMA前需要将第6列替换为对应要分析的表型(注意表型要和id号对应起来,这里可以使用R语言的merge函数或者tidyverse包)。
2. 协变量文件
在GWAS分析中,我的所感兴趣的是标记和表型之间的关联性,但除了标记外,场、年、季等固定效应同样会对我们的目标性状产生影响,虽然这些不是我们的研究的主要目标,但会对我们研究的结果产生影响。一般会将这些效应作为协变量在模型中进行校正。在GEMMA中,可以不放入协变量文件,如果有协变量文件,第一列必须添加一列全为1的截距项。
具体运行代码如下:
#计算kinship矩阵
gemma -bfile qc_data -gk 2 -o kin
#将kinship矩阵移动至当前目录
mv ./output/kin.sXX.txt.
#进行LMM分析
gemma -bfile qc_data -k kin.sXX.txt -lmm 1 -c gemma_cov.txt -o GE_GWAS
参数注解:
-bfile:plink二进制文件名称
-gk:生成kinship矩阵。选项1代表中心化的kinship矩阵,选项2为标准化的kinship矩阵。
-o:输出文件名称
-k:kinship文件名称
-lmm:进行混合线性模型GWAS分析。选项代表不同的检验方式:1为Wald检验;2为似然比检验;3为score检验;4为以上三个检验全部执行。一般使用wald检验即可。
-c:协变量文件名称
这里说明下:GEMMA进行GWAS分析时要求基因型不能有缺失,如果存在缺失,GEMMA不会计算缺失率>5%的位点。对于缺失率在0-5%的位点,GEMMA会用该位点在群体中的平均值取替代该位点的缺失值。因为我们在plink过滤时,只删除了缺失率>10%的位点,因此最终算出来的结果会比plink过滤后的点要少。在此,还是推荐大家使用Beagle等填充软件对缺失基因型进行填后再进行GWAS分析。GEMMA从生成kinship到完成GWAS分析,速度非常快。不过GEMMA也有缺点,就是只会输出三个文件:kinship文件、结果文件和日志文件。对于画Manhattan图等就需要自己写代码或者借助一些其他的分析工具了。