plink对100个性状进行gwas分析

科技   科学   2024-08-01 22:05   河南  

大家好,我是邓飞。

GWAS分析时,3~5个性状是正常操作,要分析100个性状呢,手动修改参数,工作量是够了,但是程序员的修养体现在哪里了???

如果还是按照每个性状一个文件夹,每个文件夹中一个脚本,不断地修改脚本,一点也不高端,所以,遇到这种情况,批量处理就派上用场了。

之所以之前一直不用,因为10个性状一下,没有必要,费心思想还不如直接动手操作了,但是100个性状真的吓到我了,不满足才能有进步。就看了一下参数说明,然后五分钟搞定了。虽然五分钟搞定的事情,但是写博客20分钟记录一下还是有必要的,独乐乐不如众乐乐。

核心代码:

for i in {1..100};do echo "nohup plink --file b  --allow-no-sex --pheno mphe.txt --linear --out y_${i}_result --mpheno $i "|bash;done


下面开始详细介绍。

plink中其实没有多性状模型的参数,但是它有一个--mpheno,指定性状所在的列,我们可以借用。

数据来源,GWAS Cookbook的GWAS-dat2(领取方法:快来领取 | 飞哥的GWAS分析教程),用下面代码生成表型数据:

library(data.table)

dd = fread("phe.txt")
head(dd)

set.seed(123)
xx = rnorm(150000)
nn = matrix(xx,1500,100) %>% as.data.frame()
nn[1:10,1:10]


dd1 = cbind(dd,nn)
dd1[1:10,1:10]

fwrite(dd1,"mphe.txt",col.names = F,quote = F,sep = " ")

1. 表型数据

表型数据:模拟100个性状,整理为txt,第一列FID,第二列ID,第三列以后为性状

2. 基因型数据

3. 单个性状建模

用linear模型(GLM):

plink --file b --pheno mphe.txt --linear --allow-no-sex --out re1

结果文件:

$ ls re1*
re1.assoc.linear  re1.log  re1.nosex

GWAS分析结果:

注意,上面基因型没有质控,所以有P值为NA的情况,正常质控的数据不会存在这种情况。

4. plink批量分析多性状gwas

for i in {1..100};do echo "nohup plink --file b  --allow-no-sex --pheno mphe.txt --linear --out y_${i}_result --mpheno $i "|bash;done

上面代码就是多性状gwas分析,代码解析:

  • for 循环,1~100,表示100个性状,分别运行
  • 正常进行gwas分析
  • --mpheno 后面参数$i,是分别运行100次gwas分析
  • --out 结果文件中,分别保存100个性状的gwas分析
  • |bash;done,是用管道符的形式运行nohup

运行过程:

运行的结果:

随便找一个性状结果:

完全没问题。搞定!!!

上面的批量运行程序,不但可以是plink,也可以是gemma,gcta,GAPIT等软件,都可以按照这种写法,非常666!


拓展:

GCTA、GEMMA也是可以批量处理100个性状的GWAS分析的,然后批量绘制GWAS结果,批量对显著性位点基因注释,批量绘制LDblock图,批量导出结果……


推荐阅读:


1,快来领取 | 飞哥的GWAS分析教程


2,飞哥汇总 | 入门数据分析资源推荐


3,数量遗传学,分享几本书的电子版


4GWAS进阶分析为何要推荐Linux系统




育种数据分析之放飞自我
本公众号主要介绍动植物育种数据分析中的相关问题, 算法及程序代码.
 最新文章