歡迎來到"bio生物信息"的世界
之前的推文GWAS: 曼哈頓圖，QQ plot 圖，膨脹系數（ manhattan、Genomic Inflation Factor）寫過如何用qqman包做曼哈頓圖。
但眾所周知（據我說知），這個包畫圖很丑，尤其想highlight顯著位點時，屎一樣的綠色讓人一言難盡。
這話不是我說是，是下面這位小伙子說的
此外，該包還有一個明顯缺點，只能highlight一個顯著位點，如下所示：
很多時候，我們想組合多個不同表型的顯著位點，並對不同的顯著位點進行highlight時，用qqman包顯然就實現不了。
因此，今天的重點就是怎么在qqman的基礎上highlight多個顯著位點，並且用不同的顏色表示。
我們想實現的效果如下圖所示：
接下來，講重點。
如何用對曼哈頓圖多個顯著位點標志不同顏色。
1 准備關聯分析文件
准備包含SNP, CHR, BP, P的文件gwasResults，如下所示：
2 准備想要highlight的顯著位點
想highlight多少個locus就准備多少個文件。
以這里為例，我想highlight三個顯著locus，就准備三個文件。
每個文件的格式如下圖所示：
set1文件：
set2文件：
set3文件：
3 畫圖
畫圖分兩步走。
第一步：
完全復制如下代碼到你的R中，不需要修改任何參數。
manhattan1<-function (x, chr = "CHR", bp = "BP", p = "P", snp = "SNP", col = c("gray10",
"gray60"), chrlabs = NULL, suggestiveline = -log10(1e-05),
genomewideline = -log10(5e-08), highlight1 = NULL, highlight2 = NULL,highlight3 = NULL, logp = TRUE,
...)
{
CHR = BP = P = index = NULL
if (!(chr %in% names(x)))
stop(paste("Column", chr, "not found!"))
if (!(bp %in% names(x)))
stop(paste("Column", bp, "not found!"))
if (!(p %in% names(x)))
stop(paste("Column", p, "not found!"))
if (!(snp %in% names(x)))
warning(paste("No SNP column found. OK unless you're trying to highlight."))
if (!is.numeric(x[[chr]]))
stop(paste(chr, "column should be numeric. Do you have 'X', 'Y', 'MT', etc? If so change to numbers and try again."))
if (!is.numeric(x[[bp]]))
stop(paste(bp, "column should be numeric."))
if (!is.numeric(x[[p]]))
stop(paste(p, "column should be numeric."))
d = data.frame(CHR = x[[chr]], BP = x[[bp]], P = x[[p]])
if (!is.null(x[[snp]]))
d = transform(d, SNP = x[[snp]])
d <- subset(d, (is.numeric(CHR) & is.numeric(BP) & is.numeric(P)))
d <- d[order(d$CHR, d$BP), ]
if (logp) {
d$logp <- -log10(d$P)
}
else {
d$logp <- d$P
}
d$pos = NA
d$index = NA
ind = 0
for (i in unique(d$CHR)) {
ind = ind + 1
d[d$CHR == i, ]$index = ind
}
nchr = length(unique(d$CHR))
if (nchr == 1) {
options(scipen = 999)
d$pos = d$BP/1e+06
ticks = floor(length(d$pos))/2 + 1
xlabel = paste("Chromosome", unique(d$CHR), "position(Mb)")
labs = ticks
}
else {
lastbase = 0
ticks = NULL
for (i in unique(d$index)) {
if (i == 1) {
d[d$index == i, ]$pos = d[d$index == i, ]$BP
}
else {
lastbase = lastbase + tail(subset(d, index ==
i - 1)$BP, 1)
d[d$index == i, ]$pos = d[d$index == i, ]$BP +
lastbase
}
ticks = c(ticks, (min(d[d$CHR == i, ]$pos) + max(d[d$CHR ==
i, ]$pos))/2 + 1)
}
xlabel = "Chromosome"
labs <- unique(d$CHR)
}
xmax = ceiling(max(d$pos) * 1.03)
xmin = floor(max(d$pos) * -0.03)
def_args <- list(xaxt = "n", bty = "n", xaxs = "i", yaxs = "i",
las = 1, pch = 20, xlim = c(xmin, xmax), ylim = c(0,
ceiling(max(d$logp))), xlab = xlabel, ylab = expression(-log10))
dotargs <- list(...)
do.call("plot", c(NA, dotargs, def_args[!names(def_args) %in%
names(dotargs)]))
if (!is.null(chrlabs)) {
if (is.character(chrlabs)) {
if (length(chrlabs) == length(labs)) {
labs <- chrlabs
}
else {
warning("You're trying to specify chromosome labels but the number of labels != number of chromosomes.")
}
}
else {
warning("If you're trying to specify chromosome labels, chrlabs must be a character vector")
}
}
if (nchr == 1) {
axis(1, ...)
}
else {
axis(1, at = ticks, labels = labs, ...)
}
col = rep(col, max(d$CHR))
if (nchr == 1) {
with(d, points(pos, logp, pch = 20, col = col[1], ...))
}
else {
icol = 1
for (i in unique(d$index)) {
with(d[d$index == unique(d$index)[i], ], points(pos,
logp, col = col[icol], pch = 20, ...))
icol = icol + 1
}
}
if (suggestiveline)
abline(h = suggestiveline, col = "blue")
if (genomewideline)
abline(h = genomewideline, col = "red")
if (!is.null(highlight1)) {
if (any(!(highlight1 %in% d$SNP)))
warning("You're trying to highlight1 SNPs that don't exist in your results.")
d.highlight1 = d[which(d$SNP %in% highlight1), ]
with(d.highlight1, points(pos, logp, col = "darkmagenta", pch = 20,
...))
}
if (!is.null(highlight2)) {
if (any(!(highlight2 %in% d$SNP)))
warning("You're trying to highlight2 SNPs that don't exist in your results.")
d.highlight2 = d[which(d$SNP %in% highlight2), ]
with(d.highlight2, points(pos, logp, col = "dodgerblue", pch = 20,
...))
}
if (!is.null(highlight3)) {
if (any(!(highlight3 %in% d$SNP)))
warning("You're trying to highlight3 SNPs that don't exist in your results.")
d.highlight3 = d[which(d$SNP %in% highlight3), ]
with(d.highlight3, points(pos, logp, col = "coral", pch = 20,
...))
}
}
第二步：
第一步復制完以后，輸入如下代碼：
manhattan1(sw, highlight1=set1, highlight2=set2,highlight3=set3)
最后，坐等收圖。
這里解釋一下，如果你想highlight四個locus的話，只需要在第一步的末尾處新增加：
if (!is.null(highlight4)) {
if (any(!(highlight4 %in% d$SNP)))
warning("You're trying to highlight4 SNPs that don't exist in your results.")
d.highlight4 = d[which(d$SNP %in% highlight4), ]
with(d.highlight4, points(pos, logp, col = "coral", pch = 20,
...))
}
畫圖時改成manhattan1(sw, highlight1=set1, highlight2=set2,highlight3=set3,highlight4=set4)
其他的以此類推。
此文章感謝evermane修改的代碼。感興趣的請看鏈接的討論：https://github.com/stephenturner/qqman/issues/38