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

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

「文件」

HapMap_3_r3_11.bed HapMap_3_r3_11.bim HapMap_3_r3_11.fam HapMap_3_r3_11.log1. 查看数据

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

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

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

这里的数据：

基因型个体：110个SNP个数：10737432. 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 ratioCHR 染色体SNP SNPIDBP 物理位置。。。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.log4.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.logistic5.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结果：」

「--logistic结果：」

6. 注意

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

7. 总结

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

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

GWAS学习笔记系列：

笔记 | GWAS 操作流程1：下载数据

笔记 | GWAS 操作流程2-1：缺失质控

笔记 | GWAS 操作流程2-2：性别质控

笔记 | GWAS 操作流程2-3：MAF过滤

笔记 | GWAS 操作流程2-4：哈温平衡检验

笔记 | GWAS 操作流程2-5：杂合率检验

笔记 | GWAS 操作流程2-6：去掉亲缘关系近的个体

---来自腾讯云社区的---邓飞