大家好,我是邓飞。
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图,批量导出结果……
推荐阅读: