笔记 | GWAS 操作流程3：plink关联分析--完结篇

邓飞

发布于 2020-05-13 20:54:07

7K1

发布于 2020-05-13 20:54:07

飞哥注：具体数据文件，查看之前的系列博文。这里使用的是清洗后的数据。

「文件」

HapMap_3_r3_11.bed  HapMap_3_r3_11.bim  HapMap_3_r3_11.fam  HapMap_3_r3_11.log

1. 查看数据

这里，文件是bed文件，二进制不方便查看，我们将其转化为ped文件和map文件

注意，这里我使用的是ped和map格式，如果ped文件中有表型数据（第六列），如果想指定表型数据，用--pheno，包括三列：家系，个体，表型值。

plink --bfile HapMap_3_r3_11 --recode --out test

这里的数据：

基因型个体：110个
SNP个数：1073743

2. plink关联分析的类型

❝**参考：**https://www.jianshu.com/p/286050959dbd?utm_campaign=haruki ❞

2.1 阈值性状（1,2）

plink的语境叫“case and control”，即表型值数据是两类数据：1,2，其中0和-9都表示缺失。可以选择的方法有卡方检验和逻辑斯蒂回归（X2关联分析和logistic分析）。

「--assoc」，不允许有协变量
「--logistic」，允许有协变量，如果考虑协变量，速度变慢。比assoc速度慢。

2.2 连续性状（定量性状）

这里的性状是连续性状，也就是除了1,2，0，-9外还有其它数值，「--assoc」会进行T检验（Student's test），还可以用**--linear**进行分析。

「--assoc」，不允许有协变量，速度快
「--linear」，允许有协变量，速度慢

3. 控制假阳性

因为plink进行关联分析时常常面对的是大量的SNP数据，容易产生假阳性，因此需要矫正。

Bonferroni，使用0.05/n计算出矫正后的p值作为阈值，其中n为检测SNP的个数。缺点是SNP之间可能存在LD，而且这种方法过于保守，容易产生假阴性。
FDR，是一种最小化假阳性预测比例的方法

plink的解决方法是--adjust，生成多种类型的p值。

「文件说明：」

.*.adjusted (basic multiple-testing corrections)
Produced by --adjust.

A text file with a header line, and then one line per set or polymorphic variant with the following 8-11 fields:

CHR	Chromosome code. Not present with set tests.
'SNP'/'SET'	Variant/set identifier
UNADJ	Unadjusted p-value
GC	Devlin & Roeder (1999) genomic control corrected p-value. Requires an additive model.
QQ	P-value quantile. Only present with 'qq-plot' modifier.
BONF	Bonferroni correction
HOLM	Holm-Bonferroni (1979) adjusted p-value
SIDAK_SS	Šidák single-step adjusted p-value
SIDAK_SD	Šidák step-down adjusted p-value
FDR_BH	Benjamini & Hochberg (1995) step-up false discovery control
FDR_BY	Benjamini & Yekutieli (2001) step-up false discovery control
Variants/sets are sorted in p-value order. (As a result, if the QQ field is present, its values just increase linearly.)

UNADJ：原始p值
GC：基因组矫正P值（依赖加性模型）
QQ：P-value的QQ图
BONF：Bonferroni 矫正结果
HOLM：Holm-Bonferroni (1979) adjusted p-value，矫正结果
FDR_BY Benjamini & Yekutieli (2001) step-up false discovery control，FDR方法

4. 阈值性状关联分析

数据：观测值一列是1和2，可以用的方法有：--assoc和--logistic

4.1 accoc关联分析未矫正

plink --file test --assoc --out result

「结果文件说明：」

.assoc, .assoc.fisher (case/control association allelic test report)
Produced by --assoc acting on a case/control phenotype.

A text file with a header line, and then one line per variant typically with the following 9-10 fields:

CHR	Chromosome code
SNP	Variant identifier
BP	Base-pair coordinate
A1	Allele 1 (usually minor)
F_A	Allele 1 frequency among cases
F_U	Allele 1 frequency among controls
A2	Allele 2
CHISQ	Allelic test chi-square statistic. Not present with 'fisher'/'fisher-midp' modifier.
P	Allelic test p-value
OR	odds(allele 1 | case) / odds(allele 1 | control)
If the 'counts' modifier is present, the 5th and 6th fields are replaced with:

C_A	Allele 1 count among cases
C_U	Allele 1 count among controls
If --ci 0.xy has also been specified, there are three additional fields at the end:

SE	Standard error of odds ratio estimate
Lxy	Bottom of xy% symmetric approx. confidence interval for odds ratio
Hxy	Top of xy% approx. confidence interval for odds ratio

CHR 染色体
SNP SNPID
BP 物理位置
。。。
CHISQ 卡方值
P P值结果文件：

4.2 assoc关联分析矫正p值

 plink --file test --assoc --adjust --out result_adjust

结果生成三个文件：

result_adjust.assoc           result_adjust.log
result_adjust.assoc.adjusted

查看矫正后的值：

4.3 logistic不考虑协变量

plink --file test --logistic --out result_logistic

结果生成：

result_logistic.assoc.logistic  result_logistic.log

4.4 logistic 不考虑协变量矫正p值

plink --file test --logistic --adjust --out result_logistic_adjust

结果：

result_logistic_adjust.assoc.logistic           result_logistic_adjust.log
result_logistic_adjust.assoc.logistic.adjusted

查看矫正值

如果考虑协变量，加上--covar即可。

5. 可视化

「把文件改名字，方便后面代码里作图，这样不用修改代码了。」

cp result.assoc assoc_results
cp result_logistic.assoc.logistic logistic_results.assoc.logistic

注意，logistic里面会有NA，所以我们这里讲NA的行去掉。

awk '!/'NA'/' logistic_results.assoc.logistic >logistic_results.assoc2.logistic

5.1 曼哈顿图

R代码：注意，如果没有安装qqman，就运行代码：install.packages("qqman")，安装即可。

#install.packages("qqman",repos="http://cran.cnr.berkeley.edu/",lib="~" ) # location of installation can be changed but has to correspond with the library location
#library("qqman",lib.loc="~")
library(qqman)
results_log <- read.table("logistic_results.assoc_2.logistic", head=TRUE)
jpeg("Logistic_manhattan.jpeg")
manhattan(results_log,chr="CHR",bp="BP",p="P",snp="SNP", main = "Manhattan plot: logistic")
dev.off()

results_as <- read.table("assoc_results.assoc", head=TRUE)
jpeg("assoc_manhattan.jpeg")
manhattan(results_as,chr="CHR",bp="BP",p="P",snp="SNP", main = "Manhattan plot: assoc")
dev.off()

「assoc的卡方检验结果：」

「--logistic逻辑斯蒂回归结果：」

可以看到两者结果类似，assoc检测到了显著性位点，logistic没有检测到显著性位点。

5.2 QQ plot

R代码：

#install.packages("qqman",repos="http://cran.cnr.berkeley.edu/",lib="~" ) # location of installation can be changed but has to correspond with the library location
#library("qqman",lib.loc="~")
library(qqman)
results_log <- read.table("logistic_results.assoc_2.logistic", head=TRUE)
jpeg("QQ-Plot_logistic.jpeg")
qq(results_log$P, main = "Q-Q plot of GWAS p-values : log")
dev.off()

results_as <- read.table("assoc_results.assoc", head=TRUE)
jpeg("QQ-Plot_assoc.jpeg")
qq(results_as$P, main = "Q-Q plot of GWAS p-values : log")
dev.off()

「--assoc结果：」