文件
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
image
这里的数据：
- 基因型个体：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
image
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
查看矫正后的值：
image
4.3 logistic不考虑协变量
plink --file test --logistic --out result_logistic
结果生成：
result_logistic.assoc.logistic result_logistic.log
image
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
查看矫正值
image
如果考虑协变量，加上--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的卡方检验结果：
image
–logistic逻辑斯蒂回归结果：
image
可以看到两者结果类似，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结果：
image
–logistic结果：
image
6. 注意
注意，这是阈值性状的结果，分类性状，可以使用assoc和logistic，连续性状的话，如果没有协变量就用assoc，如果有协变量，就用linear即可。
7. 总结
这是使用plink计算GWAS分析的流程，包括数据的清洗，以及建模，以及出结果，以及可视化。
但是，这只是一个开始，在看一下GEMMA，R中的GAPIT的操作方法，最后能够通过编码进行GWAS的分析，才算是掌握了这个分析方法。
