一步步教你如何绘制GWAS中的曼哈顿图以及QQ图
简介
GWAS研究中的曼哈顿图以及QQ图可以说是标配了，至于具体如何理解这两个图，可以参考我以前的博客。上周我利用自己的数据跑了个GWAS，我是利用GAPIT跑的，出的图感觉真心丑，大家可以看看：
所以我利用GAPIT出来的数据重新画了图，用的就是包qqman。
安装
Github主页上提供了三种安装方式：
# Install the stable release from CRAn
install.packages("qqman")
# 直接从Github上安装
devtools::install_github("stephenturner/qqman")
# 从Github上安装最新的开发版
devtools::install_github("stephenturner/qqman",ref="dev")
数据格式
qqman里内置了一套数据集gwasResults
library(qqman)
dim(gwasResults)
## [1] 16470 4
head(gwasResults)
## SNP CHR BP P
## 1 rs1 1 1 0.9148060
## 2 rs2 1 2 0.9370754
## 3 rs3 1 3 0.2861395
## 4 rs4 1 4 0.8304476
## 5 rs5 1 5 0.6417455
## 6 rs6 1 6 0.5190959
str(gwasResults)
## 'data.frame': 16470 obs. of 4 variables:
## $ SNP: chr "rs1" "rs2" "rs3" "rs4" ...
## $ CHR: int 1 1 1 1 1 1 1 1 1 1 ...
## $ BP : int 1 2 3 4 5 6 7 8 9 10 ...
## $ P : num 0.915 0.937 0.286 0.83 0.642 ...
可以看到该数据集共有16470行，4列，分别为SNP名、染色体、SNP位置以及P值。 该数据集不大，SNP标记十分少，通过每条染色体上的SNP标记可以看出
as.data.frame(table(gwasResults$CHR))
## Var1 Freq
## 1 1 1500
## 2 2 1191
## 3 3 1040
## 4 4 945
## 5 5 877
## 6 6 825
## 7 7 784
## 8 8 750
## 9 9 721
## 10 10 696
## 11 11 674
## 12 12 655
## 13 13 638
## 14 14 622
## 15 15 608
## 16 16 595
## 17 17 583
## 18 18 572
## 19 19 562
## 20 20 553
## 21 21 544
## 22 22 535
Manhattan plot
qqman提供了一个绘制Manhattan plot的函数manhattan()。
manhattan(gwasResults)
manhattan()提供了大量的参数设置：标题(main=)、y轴范围(ylim=)、控制点的大小(cex=)、x轴坐标轴标签的字体大小(cex.axis=)、颜色(col=)、阈值参考线控制(suggestiveline=, genomewideline=)、以及添加注释信息等。
manhattan(gwasResults, main="Manhattan plot", ylim=c(0, 10), cex=0.6, cex.axis=0.9, col = c("blue","orange"), suggestiveline = F, genomewideline = F, chrlabs = c(1:20, "P","Q"))
还提供针对一条染色体数据进行绘图
manhattan(subset(gwasResults, CHR==3))
对有显著性影响的SNP进行高亮(没有的话会忽略掉)
manhattan(gwasResults, highlight = snpsOfInterest)
针对特定染色体以及特定区间进行高亮
manhattan(subset(gwasResults, CHR==3), highlight=snpsOfInterest, xlim=c(200, 500), main="Chr 3")
manhattan注释的时候默认注释每条染色体中P值最小且超过我们设定的p值阈值的SNP点，默认的P值阈值为0.01
manhattan(gwasResults, annotatePval = 0.01)
我们可以通过自行设定参数注释所有超过P值阈值的SNP点
manhattan(gwasResults, annotatePval = 0.001, annotateTop = FALSE)
manhattan可以绘制任意的value值，不限于p-value，只需要在设置中指定value值就行了。
gwasResults <- transform(gwasResults, zscore=qnorm(P/2, lower.tail = FALSE))
head(gwasResults)
## SNP CHR BP P zscore
## 1 rs1 1 1 0.9148060 0.1069785
## 2 rs2 1 2 0.9370754 0.0789462
## 3 rs3 1 3 0.2861395 1.0666287
## 4 rs4 1 4 0.8304476 0.2141275
## 5 rs5 1 5 0.6417455 0.4652597
## 6 rs6 1 6 0.5190959 0.6447396
接下来利用value值zscore来绘制图形
manhattan(gwasResults, p = "zscore", logp = FALSE, ylab="Z-score", genomewideline = FALSE, suggestiveline = FALSE, main="Manhattan plot of Z-scores")
需要注意的几点：
- 数据集列名默认是SNP、CHR、BP和P，因此如何自己的数据集列名与此不一致的话，要么转为一致，要么绘图时指定chr=、bp=、p=和snp=。具体情况可见help(manhattan)。
- 染色体编号那一列必须是数值型，如果数据集中有"X"、“Y"以及"MT"等染色体，想要转为数值型编号再在绘图时指定坐标名
manhattan只提供了修改snp点的颜色参数，如果要修改阈值线、高亮、注释等的颜色，需要修改源码。
QQ图
qqman提供了绘制QQ图的函数qq()
qq(gwasResults)
通过参数设置可以指定点的类型、大小、颜色等
qq(gwasResults$P, main="Q-Q plot of GWAS p-value", xlim=c(0,7), ylim=c(0,12), pch=18, col = "blue4", cex=1.5, las=1)
SessionInfo
sessionInfo()
## R version 3.4.3 (2017-11-30)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 16299)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Chinese (Simplified)_China.936
## [2] LC_CTYPE=Chinese (Simplified)_China.936
## [3] LC_MONETARY=Chinese (Simplified)_China.936
## [4] LC_NUMERIC=C
## [5] LC_TIME=Chinese (Simplified)_China.936
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] qqman_0.1.3
##
## loaded via a namespace (and not attached):
## [1] compiler_3.4.3 backports_1.1.2 magrittr_1.5 rprojroot_1.3-2
## [5] tools_3.4.3 htmltools_0.3.6 yaml_2.1.16 Rcpp_0.12.14
## [9] calibrate_1.7.2 stringi_1.1.6 rmarkdown_1.8 knitr_1.18
## [13] stringr_1.2.0 digest_0.6.13 evaluate_0.10.1
