这一期内容是GWAS实战的重点部分,小陈会教大家如何简单使用PLINK这个软件完成一个常规的GWAS分析。
首先把咱们之前做的ped和map文件放到plink软件的目录下,这里我们可以使用dir这个指令查看,如下图所示:
然后执行如下指令:
plink.exe --file myWES_chr2 --make-bed --out myWES_chr2
最后会生成myWES_chr2.bed, myWES_chr2.bim和myWES_chr2.fam这三个文件,如下图所示:
接着,我们执行如下指令计算这些样本的遗传主成分:
plink.exe --bfile myWES_chr2 --pca 5 --out myWES_chr2
这里的参数--pca 5表示的是输出前五个主成分,然后我们可以得到myWES_chr2.eigenval和myWES_chr2.eigenvec这两个文件,其中.eigenvec文件储存着具体到每个样本的主成分数据,是后续矫正的部分:
通过主成分结果我们可以确定人群分层的情况,从而确定需要矫正的主成分个数(一般矫正的主成分多为5~10个)。
接下来,我们可以制作一下协变量文件,这里我们以sex,age和前五个主成分为协变量,注意PLINK要求协变量文件的前两列必须是FID和IID。
代码语言:javascript复制library(data.table)
sample_info <- fread('./sample_info.csv')
sample_info$sex <- -9
sample_info$sex[which(sample_info$`gender:ch1`=="Male")] <- 1
sample_info$sex[which(sample_info$`gender:ch1`=="Female")] <- 2
sample_info$smoking <- -9
sample_info$smoking[which(sample_info$`smoking_status:ch1`=="Non-smoker")] <- 1
sample_info$smoking[which(sample_info$`smoking_status:ch1`=="Smoker")] <- 2
sample_info$FID <- sample_info$id
sample_info$IID <- sample_info$id
sample_info$PID <- 0
sample_info$MID <- 0
pca <- fread("C:/Users/86151/Downloads/plink/myWES_chr2.eigenvec", header=F)
colnames(pca) <- c('FID','IID','PC1','PC2','PC3','PC4','PC5')
mycovar <- sample_info[,c("FID","sex","age:ch1")]
colnames(mycovar) <- c('FID','sex','age')
mycovar <- merge(mycovar, pca, by='FID')
mycovar <- mycovar[,c('FID','IID','sex','age','PC1','PC2','PC3','PC4','PC5')]
fwrite(mycovar, 'mycovar.tsv', sep='t')
接下来把做好的mycovar.tsv文件放到plink软件所在的目录底下,
最后,我们执行如下命令即可得到一个简单的GWAS summary结果:
plink.exe --bfile myWES_chr2 --maf 0.01 --hwe 2e-6 --covar mycovar.tsv --covar-number 1-7 --logistic --out myWES_chr2
注意第一个参数—bfile是输入文件的前缀,--maf是指最小等位基因的频率阈值(我们自动过滤掉maf小于0.01的SNP),--hwe是指哈迪温伯格平衡的p值阈值(我们自动过滤掉哈迪温伯格平衡的p值小于2e-6的SNP),--covar就是协变量文件名,--covar-number是指我们选择第1到7列为协变量(计算列数时自动跳过FID和IID这两列),--logistic表示使用逻辑回归模型(因为咱们的表型是吸烟和不吸烟这样的二分类变量),--out表示输出结果的前缀。最后我们会得到以.assoc.logistic结尾的文件,这就是我们的结果文件。
接下来,我们可以在CMD里输入如下指令查看结果文件的前几行数据:
more myWES_chr2.assoc.logistic
最后,使用“Q”健退出查看。
关于如何使用PLINK进行GWAS分析就先介绍到这里,下期我会带大家解读一下结果,敬请期待!