PLINK-GWAS学习9------对于二元数据的关联分析

前期的准备文件：上一个教程中的HapMap_3_r3_13 (with .bed, .bim, and .fam. extensions) and covar_mds.txt
对于二元形状

plink --bfile HapMap_3_r3_13 --assoc --out assoc_results

注意，——assoc选项不允许纠正协变量，如主成分(PC的)/ MDS成分，这使得它不适合关联分析。
plink的–logistic命令允许加入协变量。协变量文件就是上一个教程中产生的协变量文件，包含是个成分。

plink --bfile HapMap_3_r3_13 --covar covar_mds.txt --logistic --hide-covar --out logistic_results

使用–hide-covar命令可以只在输出文件中显示SNPs的附加结果。
剔除NA结果，这会影响之后的作图步骤。

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

后续即可作图了。
需要注意的是，对于数量性状，需要将–logistic命令替换为–linear。

多重检验。
除了常规的全基因组显著性阈值5.0e-8之外，还有多种方法进行多种检测。
adjust

plink --bfile HapMap_3_r3_13 -assoc --adjust --out adjusted_assoc_results

这个文件给出了一个Bonferroni校正后的p值，以及FDR和方法校正的P值。

permutation

# Generate subset of SNPs
awk '{ if ($4 >= 21595000 && $4 <= 21605000) print $2 }' HapMap_3_r3_13.bim > subset_snp_chr_22.txt
# Filter your bfile based on the subset of SNPs generated in the step above.
plink --bfile HapMap_3_r3_13 --extract subset_snp_chr_22.txt --make-bed --out HapMap_subset_for_perm
# Perform 1000000 perrmutations.
plink --bfile HapMap_subset_for_perm --assoc --mperm 1000000 --out subset_1M_perm_result
# Order your data, from lowest to highest p-value.
sort -gk 4 subset_1M_perm_result.assoc.mperm > sorted_subset.txt
# Check ordered permutation results
head sorted_subset.txt

这个的具体用法不清楚，但是文章中的描述是比上述两种P值的校正方法（Bonferroni和FDR）更加高效。

制图，qq.man

install.packages("qqman") # location of installation can be changed but has to correspond with the library location 
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()  
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的结果