GWAS实战教程之利用PLINK进行GWAS分析

这一期内容是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。

‍

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分析就先介绍到这里，下期我会带大家解读一下结果，敬请期待！