笔记 GWAS 操作流程3：plink关联分析_plink进行关联分析

文件

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

1
2

1. 查看数据

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

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

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

1
2

这里的数据：

• 基因型个体：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.)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

• 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

1
2

结果文件说明：

.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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

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

4.2 assoc关联分析矫正p值

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

1
2

结果生成三个文件：

result_adjust.assoc           result_adjust.log
result_adjust.assoc.adjusted 
1
2

查看矫正后的值：

4.3 logistic不考虑协变量

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

1
2

结果生成：

result_logistic.assoc.logistic  result_logistic.log

1
2

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

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

1
2

结果：

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

1
2
3

查看矫正值

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

5. 可视化

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

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

1
2
3

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

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

1
2

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() 
1
2
3
4
5
6
7
8
9
10
11
12

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()
1
2
3
4
5
6
7
8
9
10
11
12

–assoc结果：

–logistic结果：

6. 注意

注意，这是阈值性状的结果，分类性状，可以使用assoc和logistic，连续性状的话，如果没有协变量就用assoc，如果有协变量，就用linear即可。

7. 总结

这是使用plink计算GWAS分析的流程，包括数据的清洗，以及建模，以及出结果，以及可视化。

但是，这只是一个开始，在看一下GEMMA，R中的GAPIT的操作方法，最后能够通过编码进行GWAS的分析，才算是掌握了这个分析方法。