用plink做GWAS（PCA、关联分析）并用R绘图

GWAS

主要是做质量控制、PCA分析和关联分析

一、观察初始数据

初始数据一般有不同格式的，最终都要转化为后缀为*.bed，*.bim 和 *.fam的三个文件，首先要学会如何看初始文件判断数据类型，判断好数据类型才能选择合适的关联分析方法。那么关于这三种文件可以看：初探PLINK文件格式（bed，bim，fam）.
一般来说fam文件里会有表型信息，或者是单独一个pheno文件给出表型信息。如果是单独pheno文件的话，需要先查看文件，找到需要的表型。这里我以做HDL_High_density_lipoprotein__mmol_l_关联分析为例，表型数据由T_pheno文件单独提供。我们先来查看一下T_pheno文件：
可以看出我们需要的表型HDL处于第四列，前两列是固定不变的FID和IID，实际表型从第三列数起，也就是说HDL是第二个表型，之后指明表型时可以用–mpheno 2或–pheno-name HDL_High_density_lipoprotein__mmol_l_ 两种参数指明。
同时，我们观察到表型数据是连续型，而非0和1这种二进制类型，所以我们选择用linear regression做关联分析。（binary trait主要用logistic regression，quantitive trait主要用linear regression，要先观察数据类型以免用错方法）

二、质量控制

过程参考了用plink做一套GWAS分析，选择了自己需要的方法。这里讲解一下很多人刚开始接触的时候不明白各种参数是做什么的，比如样本缺失率和位点缺失率的过滤为什么有–missing和–mind两种，他们之间的区别是：–missing是用来筛查有哪些样本和位点缺失了，输出文件是将有缺失的数据整合出来，提供给我们阅读。而–mind是用于对原始数据进行修改，将高于你设置的缺失率的数据删除，并不会告诉你删除了哪些数据，只会输出删除后的结果。所以我们一般用–missing查看缺失情况，用–mind来做具体删除过滤操作。（同理–freq和–maf）
我做的过程如下：

样本缺失率和位点缺失率过滤（产生.imiss和lmiss文件）

初始数据是example.bed、example.bim、example.fam、T_pheno。

plink --bfile example --missing

自动生成plink.imiss和plink.lmiss文件，可以less+文件名查看一下。

计算MAF

plink --bfile example --freq

自动生成plink.freq文件，可以less+文件名查看一下。

数据清理

这里为了不产生那么多文件，就把参数合在同一条指令，一次性输出过滤结果。具体数值是按照默认标准设定的。

plink –bfile example --mind 0.1 --maf 0.05 --geno 0.1 --hwe 0.01 --make-bed --out clean

各参数含义：
Missingness per individual --mind 0.1 >10% 排除
Missingness per marker --geno 0.1 > 10% 排除
Allele frequency --maf 0.05 MAF <= 0.05 排除
Hardy-Weinberg equilibrium --hwe 0.001 case组合control组会分开做，按照各自数据集的情况排除

最后输出的是clean.bed、clean.bim、clean.fam。
接下来我们就用过滤后的数据进行操作。

三、检查亲缘关系

需要无关个体数据的实验要做亲缘关系检查IBS

plink --bfile clean --genome --min 0.2 --pheno T_pheno --mpheno 2 --allow-no-sex

–genome是用于亲缘关系检查的参数，–pheno用于指明表型文件，–mpheno指明所需表型所在列，–allow-no-sex是忽略性别（原因是我的初始数据中性别均为未知，根据初始数据具体情况选择加不加这个参数）
输出文件为plink.genome，具体如何看genome文件，可以看GWAS学习笔记（一）：质量控制（QC）
我主要是看RT这一列，是UN则是无关个体，那么就无需进行过滤。如果有血缘关系的个体可以通过–remove参数进行删除。

四、PCA分析

plink --threads 16 --bfile clean --pca 10 --out pca --allow-no-sex --pheno T_pheno --pheno-name HDL_High_density_lipoprotein__mmol_l 

–pca参数用于对数据做pca分析，后面的10是取前十个pca结果的意思，一般不需要取太多，因为最后是用最显著的第一列和第二列作图。–threads是为了增加线程，–out参数后面跟输出文件名，可以自己设定。最终得到两个文件：.eigenval 和.eigenvec，我们用.eigenvec文件进行绘图。
绘图软件我选择的是Rstudio，安装了ggplot2这个包：

>install.packages('ggplot2')

先用文本编辑器打开.eigenvec文件，在第一行插入行，输入1 2 pc1 pc2 pc3 …
每个用空格间隔开，具体个数和你的列数相匹配，如果是之前设了10，那就要写到pc10，目的是给出列名，方便Rstudio查找具体列，修改完后保存文件。
导入包和数据：

>library(ggplot2)
>data<-read.table("D:\pca.eigenvec",header=T)

具体文件路径根据你的实际情况来。
绘制PCA图：

> ggplot(data,aes(x=pc1,y=pc2))+ geom_point()

在右边plots会显示图片，大致为这样：

五、关联分析

前面说过，观察了pheno文件是连续型表型数据，所以我们关联分析方法就选择线性回归（linear regression）。

plink --bfile clean --linear --pheno T_pheno --mpheno 2 --allow-no-sex

得到plink.assoc.linear文件
用Rstudio绘制曼哈顿图和QQ图：
安装qqman包：

>install.packages('qqman')

导入包和数据

>library(qqman)
>data<-read.table("D:\plink.assoc.linear",header=T)

绘制曼哈顿图

> color_set <- rainbow(9)
> jpeg("Linear_manhattan.jpeg")
> manhattan(results1,chr="CHR",bp="BP",p="P",snp="SNP",col=color_set, main = "Manhattan plot: linear")
> dev.off()

我选了彩虹颜色作为色系，也可以根据自己喜欢设置
画出的图大概是这样
最高点就是最显著
绘制qq图

> jpeg("QQ-Plot_linear.jpeg")
> qq(results1$P, main = "Q-Q plot of GWAS p-values : log")
> dev.off()

上面这个是最简单的绘制，如需要标注λ值可以参考下面这个：

> jpeg("QQ-Plot_linear.jpeg")
> p_value=results1$P
> z = qnorm(p_value/ 2)
> lambda = round(median(z^2, na.rm = TRUE) / 0.454, 3)
> qq(results1$P, main = "Q-Q plot of GWAS p-values : log",sub=paste("lamda=",lambda))
> dev.off()

画出的图大概是这样
观察偏离趋势

写在最后

至此，一套用plink做GWAS的流程就大概是这样。笔者是这两天才第一次接触学习生物信息，这个博客是根据自己的理解写的，如有错误欢迎指正。plink的操作我都是在linux下完成，而R绘图是在windows下，觉得文件传输有点麻烦，如有更好的办法可以和我交流，感激不尽！