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

GWAS
- 一、观察初始数据
- 二、质量控制
- - 样本缺失率和位点缺失率过滤（产生.imiss和lmiss文件）
  - 计算MAF
  - 数据清理
- 三、检查亲缘关系
- 四、PCA分析
- 五、关联分析
- 写在最后

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下，觉得文件传输有点麻烦，如有更好的办法可以和我交流，感激不尽！

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

【生信】全基因组关联分析（GWAS）原理
[生信]全基因组关联分析(GWAS)原理文章的文字/图片/代码部分/全部来源网络或学术论文,文章会持续修缮更新,仅供大家学习使用. 目录 [生信]全基因组关联分析(GWAS) 1.前提知识介绍 1. ...
GWAS | 全基因组关联分析 | PLINK | 实战 | 统计遗传学
参考:PLINK | File format reference vcftools plink的主要功能:数据处理,质量控制的基本统计,群体分层分析,单位点的基本关联分析,家系数据的传递不平衡检验,多 ...
GWAS全基因组关联分析流程（BWA+samtools+gatk+Plink+Admixture+Tassel）
我梳理了GWAS全基因组关联分析的整个流程,并提供了基本的命令,用到的软件包括BWA.samtools.gatk.Plink.Admixture.Tassel等,在此分享出来给大家提供参考. 一.BW ...
GWAS | 原理和流程 | 全基因组关联分析 | Linkage disequilibrium (LD)连锁不平衡 | 曼哈顿图 Manhattan_plot | QQ_plot |...
问题: linkage disequilibrium (LD)和 pairwise correlation的区别?似乎它们都能达到相同的目的. 先从直觉上理解一下GWAS的原理: 核心就是SNP与表型 ...
DNA 12. SCI 文章绘图之全基因组关联分析可视化(GWAS)
点击关注,桓峰基因桓峰基因生物信息分析,SCI文章撰写及生物信息基础知识学习:R语言学习,perl基础编程,linux系统命令,Python遇见更好的你 134篇原创内容公众号桓峰基因公众号推 ...
2020.10.21【转载】丨GWAS全基因组关联分析流程
感谢CSDN用户追梦生信人梳理了GWAS全基因组关联分析的整个流程,并提供了基本的命令,用到的软件包括BWA.samtools.gatk.Plink.Admixture.Tassel等,在此分享出 ...
全基因组关联分析（GWAS）实现途径之一
转自:简书:链接:https://www.jianshu.com/p/bf7e67680414 ################################################# 写在 ...
GWAS处理流程（全基因组关联分析）——对从ADNI数据库下载的SNP数据及进行质控（QC）
对从ADNI数据库下载的SNP数据及进行质控(QC) 简介一.先查看数据中的个体和SNP缺失情况 1.查看 2.生成绘图以可视化缺失结果. 二.QC第一步:删除缺失度大于某个数值的SNP和个体 1. ...
GWAS - 基因型与表型的关联分析流程
这里介绍的SNP数据与单一表型的相关分析流程一.准备plink文件 ped和map文件转换为bed.fam.bim plink --file mydata --out mydata --make-b ...

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

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

GWAS

一、观察初始数据

二、质量控制

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

计算MAF

数据清理

三、检查亲缘关系

四、PCA分析

五、关联分析

写在最后

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

最新文章

热门文章