本文首發(fā)于“生信補給站”，ggplot2|玩轉(zhuǎn)Manhattan圖-你有被要求這么畫嗎？更多關(guān)于R語言，ggplot2繪圖，生信分析的內(nèi)容，敬請關(guān)注小號。
Manhattan圖算是GWAS分析的標(biāo)配圖了，可參考Bio|manhattan圖 進行繪制。 由于Manhattan點太多，后期AI/PS修改的話難度有點大，如果可以“個性化”繪制的話那是極好的！ 一 載入R包，數(shù)據(jù)1）載入數(shù)據(jù)處理的tidyverse包，使用qqman中g(shù)wasResults示例數(shù)據(jù)集 #載入R包
#install.packages("qqman")
library(qqman)
library(tidyverse)
#查看原始數(shù)據(jù)
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 我們知道Manhattan圖實際就是點圖，橫坐標(biāo)是chr，縱坐標(biāo)是-log(Pvalue) ，原始P值越小，-log轉(zhuǎn)化后的值越大，在圖中就越高。 原始數(shù)據(jù)中重要的“元素”都有了 ，我們自己的數(shù)據(jù)也是只需要這四列就可以了。注意繪制前需要轉(zhuǎn)化一下： 2）處理原始數(shù)據(jù)---計算SNP的累計位置# 1)計算chr長度
chr_len <- gwasResults %>%
group_by(CHR) %>%
summarise(chr_len=max(BP))
# 2） 計算每條chr的初始位置
chr_pos <- chr_len %>%
mutate(total = cumsum(chr_len) - chr_len) %>%
select(-chr_len)
#3)計算累計SNP的位置
Snp_pos <- chr_pos %>%
left_join(gwasResults, ., by="CHR") %>%
arrange(CHR, BP) %>%
mutate( BPcum = BP + total)
#查看轉(zhuǎn)化后的數(shù)據(jù)
head(Snp_pos,2)
SNP CHR BP P total BPcum
1 rs1 1 1 0.9148060 0 1
2 rs2 1 2 0.9370754 0 2 數(shù)據(jù)準(zhǔn)備完成，開始繪圖。 二 ggplot2繪制Manhattan圖1 縱坐標(biāo)為P值轉(zhuǎn)-log10（）ggplot(Snp_pos, aes(x=BPcum, y=-log10(P))) +
geom_point( aes(color=as.factor(CHR))) 
2 繪制加強版Manhattan圖1) 準(zhǔn)備X軸標(biāo)簽位置--在每條chr的中間 X_axis <- Snp_pos %>% group_by(CHR) %>% summarize(center=( max(BPcum) + min(BPcum) ) / 2 ) 2）繪制“改良版”Manhattan圖 p <- ggplot(Snp_pos, aes(x=BPcum, y=-log10(P))) +
#設(shè)置點的大小，透明度
geom_point( aes(color=as.factor(CHR)), alpha=0.8, size=1.3) +
#設(shè)置顏色
scale_color_manual(values = rep(c("grey", "skyblue"), 22 )) +
#設(shè)定X軸
scale_x_continuous( label = X_axis$CHR, breaks= X_axis$center ) +
#去除繪圖區(qū)和X軸之間的gap
scale_y_continuous(expand = c(0, 0) ) +
#添加閾值線
geom_hline(yintercept = c(6, -log10(0.05/nrow(Snp_pos))), color = c('green', 'red'), size = 1.2, linetype = c("dotted", "twodash")) +
#設(shè)置主題
theme_bw() +
theme(
legend.position="none",
panel.border = element_blank(),
axis.line.y = element_line(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()
) 
這時候是不是就可以了，??。 當(dāng)然了既然是ggplot2繪制的Manhattan圖（點圖），那么關(guān)于點，線，坐標(biāo)，主題的設(shè)置當(dāng)然都可以設(shè)置了，看這里 ggplot2|詳解八大基本繪圖要素 ggplot2|theme主題設(shè)置，詳解繪圖優(yōu)化-“精雕細琢” 3 玩轉(zhuǎn)Manhattan圖1) 利用數(shù)據(jù)集自帶的snpsOfInterest標(biāo)示顯著的位點，展示重要的基因信息library(ggrepel)
#準(zhǔn)備數(shù)據(jù)
data <- Snp_pos %>%
# 添加高亮和注釋信息：snpsOfInterest中的rs編號和P值大于6的點
mutate( is_highlight=ifelse(SNP %in% snpsOfInterest, "yes", "no")) %>%
mutate( is_annotate=ifelse(-log10(P)>6, "yes", "no"))
# 繪圖
p1 <- ggplot(data, aes(x=BPcum, y=-log10(P))) +
geom_point( aes(color=as.factor(CHR)), alpha=0.8, size=1.3) +
scale_color_manual(values = rep(c("grey", "skyblue"), 22 )) +
scale_x_continuous( label = X_axis$CHR, breaks= X_axis$center ) +
scale_y_continuous(expand = c(0, 0) ) +
# 添加高亮點
geom_point(data=subset(data, is_highlight=="yes"), color="orange", size=2) +
# 添加高亮label，且防止重疊
geom_label_repel( data=subset(data, is_annotate=="yes"), aes(label=SNP), size=2) +
theme_bw() +
theme(
legend.position="none",
panel.border = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()
) 
如果我們自己的gwas結(jié)果數(shù)據(jù)是Gene的話，label更改即可標(biāo)示基因。 2) 自定義重要的基因，標(biāo)示如果有某些“目的基因”，想查看這些基因的P值呢？ 新加gene和gene_annotate列即可！ #準(zhǔn)備數(shù)據(jù),使用基礎(chǔ)函數(shù)
data <- Snp_pos
#根據(jù)目的基因的位置，新加gene和gene_annotate列
data$gene[data$CHR == 3 & data$BP == 366] <- "geneA"
data$gene_annotate[data$CHR == 3 & data$BP == 366] <- "yes"
data$gene[data$SNP == "rs4064"] <- "geneB"
data$gene_annotate[data$SNP == "rs4064"] <- "yes"
# 繪圖
p2 <- ggplot(data, aes(x=BPcum, y=-log10(P))) +
geom_point( aes(color=as.factor(CHR)), alpha=0.8, size=1.3) +
scale_color_manual(values = rep(c("grey", "skyblue"), 22 )) +
scale_x_continuous( label = X_axis$CHR, breaks= X_axis$center ) +
scale_y_continuous(expand = c(0, 0) ) +
geom_label_repel( data=subset(data, gene_annotate=="yes"), aes(label=gene), size=4, col = "red") +
theme_bw() +
theme(
legend.position="none",
panel.border = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()
) 
3）區(qū)域放大展示重點展示某一區(qū)域的P值情況 library(ggforce)
data <- Snp_pos %>%
# 添加高亮和注釋信息：snpsOfInterest中的rs編號和P值大于6的點
mutate( is_highlight=ifelse(SNP %in% snpsOfInterest, "yes", "no")) %>%
mutate( is_annotate=ifelse(-log10(P)>6, "yes", "no"))
p3 <- ggplot(data, aes(x=BPcum, y=-log10(P))) +
geom_point( aes(color=as.factor(CHR)), alpha=0.8, size=1.3) +
scale_color_manual(values = rep(c("grey", "skyblue"), 22 )) +
scale_x_continuous( label = X_axis$CHR, breaks= X_axis$center ) +
scale_y_continuous(expand = c(0, 0) ) +
geom_point(data=subset(data, is_highlight=="yes"), color="orange", size=2) +facet_zoom(x = BPcum >= 3000 & BPcum <=3500)+
theme_bw() +
theme(
legend.position="none",
panel.border = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()
) 
可參考ggforce|繪制區(qū)域輪廓-區(qū)域放大-尋找你的“onepiece” 4）plotly 交互展示 library(plotly)
data <- Snp_pos %>%
mutate( is_highlight=ifelse(SNP %in% snpsOfInterest, "yes", "no")) %>% filter(-log10(P)>0.5) #過濾一些點，交互式壓力小
# 準(zhǔn)備SNP展示的text信息
data$text <- paste("SNP: ", data$SNP, "\nPosition: ", data$BP, "\nChromosome: ", data$CHR, "\nLOD score:", -log10(data$P) %>% round(2), "\nWhat else do you wanna know", sep="")
p4 <- ggplot(data, aes(x=BPcum, y=-log10(P), text=text)) +
geom_point( aes(color=as.factor(CHR)), alpha=0.8, size=1.3) +
scale_color_manual(values = rep(c("grey", "skyblue"), 22 )) +
scale_x_continuous( label = X_axis$CHR, breaks= X_axis$center ) +
scale_y_continuous(expand = c(0, 0) ) +
ylim(0,9) +
geom_point(data=subset(data, is_highlight=="yes"), color="orange", size=2) +
theme_bw() +
theme(
legend.position="none",
panel.border = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()
)
ggplotly(p4, tooltip="text") 
好吧，其實這個用處不太大，，， 以上就是ggplot2繪制一些常見的Manhattan圖，好處當(dāng)然就是兼容ggplot2的參數(shù)，也就可以根據(jù)需要自行設(shè)置。 ◆ ◆ ◆ ◆ ◆ 精心整理（含圖版）|你要的全拿走！有備無患 （R統(tǒng)計，ggplot2繪圖，生信圖形可視化匯總） 【覺得不錯，右下角點個“在看”，期待您的轉(zhuǎn)發(fā)，謝謝！】
|
