GAPIT使用plink数据进行GWAS分析
GAPIT使用plink数据进行GWAS分析
邓飞
发布于 2022-12-13 17:51:46
发布于 2022-12-13 17:51:46
大家好，我是邓飞。hmp格式是一种基因型格式，但是现在更多的是vcf或者plink格式的数据，今天介绍一下plink格式的数据如何导入到GAPIT软件中进行分析。
GAPIT软件支持的基因型格式为：hmp格式，plink数据转化为hmp格式，中间经过了很多路。现在提供另外一种解决方案，不用将plink数据转化为hmp格式，进行GWAS分析。
另外，如果还没有安装GAPIT软件，可以参考这篇博文：如何安装GWAS软件包：GAPIT
1. GAPIT软件的基因型数据格式：hmp
「hmp格式：」
2. GAPIT软件的基因型数据格式：Numeric格式
查看GAPIT说明文档时，发现了GAPIT还支持Numeric format，即转化为0-1-2的格式，这样就好处理了，可以使用plink软件的recodeA参数，可以很容易的将plink数据转化为0-1-2的形式。
「基因型文件：」
- 第一列是ID
- 第二列以后是基因型分型0-1-2
- 有行头
「染色体位置文件：」
- 第一列是染色体名称，需要和基因型顺序一致
- 第二列是染色体编号
- 第三列是染色体物理位置
- 有行头
3. 将plink格式转化为0-1-2的格式
c为二进制的plink文件，运行下面命令，生成plink.raw文件。
plink --bfile c --recodeA --out --re # 生成0-1-2的基因型数据
plink --bfile c --recode --out file # 生成map数据，用于raw文件命名
然后准备两个文件：re.raw和file.map文件，用下面R代码，生成GAPIT运行的文件格式。
运行示例代码，比较GAPIT和GEMMA的MLM模型的结果：
source("http://zzlab.net/GAPIT/GAPIT.library.R")
source("http://zzlab.net/GAPIT/gapit_functions.txt")
myGd = read.table("mdp_numeric.txt",header=T)
myGM = read.table("mdp_SNP_information.txt",header = T)
myY = read.table("phe.txt")
head(myY)
myGAPIT = GAPIT(Y = myY[,c(1,3)],GD = myGd, GM = myGM,PCA.total=0)
library(data.table)
re1 = fread("GAPIT.MLM.V3.GWAS.Results.csv")
re2 = fread("../10_gemma_analysis_lmm/output/result.assoc.txt")
re = merge(re1,re2,by.x="SNP",by.y="rs")
cor(re$beta,re$effect)
cor(re$P.value,re$p_wald)
「结果：」
> cor(re$beta,re$effect)
[1] 0.9979531
> cor(re$P.value,re$p_wald)
[1] 0.9913544
可以看到，effect值相关系数为0.997，P值的相关系数为0.991，两者结果一致。
❝关注我的公众号：育种数据分析之放飞自我。主要分享R语言，Python，育种数据分析，生物统计，数量遗传学，混合线性模型，GWAS和GS相关的知识。 ❞
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2022-11-15，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录