R语言制作曼哈顿图 Q-Q图等
qqman软件包简介
所述qqman包包括用于从GWAS结果创建曼哈顿地块和qq函数作图。gwasResults
软件包中包含的data.frame具有22个染色体上的16,470个SNP的模拟结果。看一下数据:
str(gwasResults)
'data.frame': 16470 obs. of 5 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 ...
$ zscore: num 0.107 0.0789 1.0666 0.2141 0.4653 ...
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
tail(gwasResults)
SNP CHR BP P zscore
16465 rs16465 22 530 0.5643702 0.5763624
16466 rs16466 22 531 0.1382863 1.4822028
16467 rs16467 22 532 0.3936999 0.8529268
16468 rs16468 22 533 0.1778749 1.3473271
16469 rs16469 22 534 0.2393020 1.1767332
16470 rs16470 22 535 0.2630441 1.1192251
每个染色体上有几个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(gwasResults)
我们还可以传入其他图形参数。我们添加一个标题(main=
),增加y轴限制(ylim=
),将点大小减小到60%(cex=
),并将轴标签的字体大小减小到90%(cex.axis=
)。在此过程中,让我们更改颜色(col=
),删除提示性的和在整个基因组范围内的重要性线,并为染色体提供我们自己的标签:
manhattan(gwasResults, main = "Manhattan Plot", ylim = c(0, 10), cex = 0.6,
cex.axis = 0.9, col = c("blue4", "orange3"), suggestiveline = F, genomewideline = F,
chrlabs = c(1:20, "P", "Q"))
现在,让我们看一个染色体:
manhattan(subset(gwasResults, CHR == 1))
让我们突出显示3号染色体上感兴趣的一些SNP。在此突出显示的100个SNP位于称为的字符向量中snpsOfInterest
。如果您尝试突出显示不存在的SNP,则会收到警告。
str(snpsOfInterest)
chr [1:100] "rs3001" "rs3002" "rs3003" "rs3004" "rs3005" ...
manhattan(gwasResults, highlight = snpsOfInterest)
我们可以结合突出显示和限制单个染色体,并使用xlim
图形参数放大感兴趣的区域(位置200-500之间):
manhattan(subset(gwasResults, CHR == 3), highlight = snpsOfInterest, xlim = c(200,
500), main = "Chr 3")
我们还可以根据SNP的p值对其进行注释。默认情况下,这仅注释每个超过annotatePval
阈值的染色体的最高SNP 。
manhattan(gwasResults, annotatePval = 0.01)
我们还可以注释所有满足阈值的SNP:
manhattan(gwasResults, annotatePval = 0.005, annotateTop = FALSE)
最后,该manhattan
函数可用于绘制任何值,而不仅仅是p值。在这里,我们仅将传递给p=
参数的函数称为要绘制的列的名称,而不是默认的“ P”列。在此示例中,让我们创建一个测试统计量(“ zscore”),绘制该值而不是p值,更改y轴标签,并删除默认的对数转换。我们还将删除全基因组和提示性行,因为这些行仅在您绘制-log10(p-values)时才有意义。
# Add test statistics
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
# Make the new plot
manhattan(gwasResults, p = "zscore", logp = FALSE, ylab = "Z-score", genomewideline = FALSE,
suggestiveline = FALSE, main = "Manhattan plot of Z-scores")
有关创建曼哈顿地块的一些注意事项:
- 运行
str(gwasResults)
。请注意,gwasResults
data.frame具有SNP,染色体位置,并命名为p值列SNP
,CHR
,BP
,和P
。如果你正在创建一个曼哈顿的情节和列名是不同的,你必须要列名传递给chr=
,bp=
,p=
,和snp=
参数。有关help(manhattan)
详细信息,请参见。 - 染色体列必须为数字。如果您具有“ X”,“ Y”或“ MT”染色体,则需要重命名这23、24、25等。您可以修改源代码(例如
fix(manhattan)
)以更改指定轴刻度的线标签(labs <- unique(d$CHR)
)将其设置为您想要的任何值。 - 如果要更改突出显示或建议性/基因组范围的线条的颜色,则需要修改源代码。搜索
col="blue"
,col="red"
或col="green3"
分别修改提示线,全基因组线,并突出显示的颜色。
创建QQ图
创建QQ图很简单-只需向qq()
函数提供p值向量即可。
qq(gwasResults$P)
我们通常可以提供许多其他图形参数。
qq(gwasResults$P, main = "Q-Q plot of GWAS p-values", xlim = c(0, 7), ylim = c(0,
12), pch = 18, col = "blue4", cex = 1.5, las = 1)