前面自己写过一篇笔记，用的是tassel（https://mp.weixin.qq.com/s?__biz=MzU0ODk5MjM3OA==&mid=2247483855&idx=1&sn=eeb6b18437eec15bd5fe7668c5809f2c&chksm=fbb7e5daccc06ccca3629fb9a405311404a11bf01e8cc554f3c03ac4249d85e197ef649972ab&token=722272578&lang=zh_CN#rd），因为最近刚好帮别人分析一批GWAS的数据，所以也顺手尝试了GEMMA，发现比tassel快很多。所以果断用了GEMMA。
从GEMMA的介绍可以看出，还是相当强大的。
GEMMA（Genome-wide Efficient Mixed Model Association algorithm）是一款基于混合线性模型的GWAS分析软件。GEMMA相比较于其他基于混合线性模型的软件，它有如下优势：
-
快速：远远快于其他精确算法（EMMA和FaST-LMM）。
-
准确：EMMAX和GAPIT都采用固定零模型中的方差组分不变的策略来提高运算速度，这实际上就是一种近似算法，不如GEMMA准确。
-
方便：可直接使用plink二进制格式数据，无需进行复杂的数据格式转换。（尤其是不像tassel运行起来参数繁琐的很。）
-
功能全面：可进行单标记GWAS、多标记GWAS和多性状GWAS分析。
====GEMMA的下载和安装========
下载地址：https://github.com/genetics-statistics/GEMMA/releases
我下载的是目前的最新版本：0.98.4（gemma-0.98.4-linux-static-AMD64.gz）这个是编译过可以直接运行的。测试一下，也是可以运行没有问题的。
===准备输入文件======
GEMMA可以接受plink格式的输入文件，所以我们先准备plink格式的文件。
因为我处理的大多数是重测序的数据，所以得到的原始文件一般是vcf格式的，所以就从vcf格式开始转化了。
vcftools --vcf all.final.vcf --min-alleles 2 --max-alleles 2 --maf 0.05 --max-missing 0.3 --minQ 20 --recode --out test //做了一个初步过滤
vcftools --vcf test.recode.vcf --plink --out out
plink --file out --make-bed --out test--noweb
这样我们就获得了bed，bim，fam三个输入文件。
===表型数据====
因为GEMMA是通过识别plink格式的fam文件中的表型来进行关联分析的，所以我们需要预处理一下，把表型数据添加进去。添加到fam文件的6，7，8...等等列。
注意：样品名称的位置要一样。另外对于缺失值用NA或者-9。下图就是填上表型之后的图。如果有多个表型，可以一次填到后面的列，但是后面关联分析的时候利用-n指定用的是那个表型进行分析的。
===开始运行GEMMA===
首先获得kinship矩阵，使用混合线性模型进行分析。
/gpfs03/home/jingjing/software/gemma/./gemma -bfile test -gk 2 -o kin
其中：-gk 2 标准化的方法计算G矩阵
然后再output文件夹下会生成kin.sXX.txt文件。
下面进行关联分析：
/gpfs03/home/jingjing/software/gemma/./gemma -bfile test -k kin.sXX.txt -lmm 1 -o yield
其中：-k 读取G矩阵的文件 -lmm 1 使用Wald的方法进行SNP检验
在output文件夹中的yield.assoc.txt中是最后运行的结果。
其中：
chr：SNP所在染色体号
rs： SNP名称
ps： SNP物理位置
n_miss： SNP缺失个体数
allele1： 次等位基因
allele0： 主等位基因
af：SNP频率
beta： SNP效应值
se： beta估计标准误
l_remle： 计算该SNP效应时对应的lamda的remle估计值。
p_wald ：wald检验P值
其中，我们最关心的三个结果是chr, ps, p_wald，我们可以借助这三个结果画曼哈顿图和QQ图。l_remle比较难理解，需要懂模型才知道它的含义，但对分析来说，不是很重要。
所以我们提取其中的chr，ps，p_wald值进行画图。
格式如下：
然后在R中画图：
library(qqman)
png("manhattan.png")
manhattan(data, col = c(1:26))
png("QQ.png")
qq(data$P)
本文使用 文章同步助手 同步
网友评论