ggplot2|玩轉Manhattan圖-你有被要求這么畫嗎？

本文轉載自查看原文 2019-12-11 10:21 518

本文首發於“生信補給站”，ggplot2|玩轉Manhattan圖-你有被要求這么畫嗎？更多關於R語言，ggplot2繪圖，生信分析的內容，敬請關注小號。

Manhattan圖算是GWAS分析的標配圖了，可參考Bio|manhattan圖進行繪制。

由於Manhattan點太多，后期AI/PS修改的話難度有點大，如果可以“個性化”繪制的話那是極好的！

一載入R包，數據

1）載入數據處理的tidyverse包，使用qqman中gwasResults示例數據集

#載入R包
#install.packages("qqman")
library(qqman)
library(tidyverse)
#查看原始數據
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圖實際就是點圖，橫坐標是chr，縱坐標是-log(Pvalue) ，原始P值越小，-log轉化后的值越大，在圖中就越高。

原始數據中重要的“元素”都有了，我們自己的數據也是只需要這四列就可以了。注意繪制前需要轉化一下：

2）處理原始數據---計算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)

#查看轉化后的數據
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

數據准備完成，開始繪圖。

二 ggplot2繪制Manhattan圖

1 縱坐標為P值轉-log10（）

ggplot(Snp_pos, aes(x=BPcum, y=-log10(P))) +
geom_point( aes(color=as.factor(CHR)))

基本圖形出來了，但是有點怪；不急，一點點改進：
- 橫坐標標簽設置在每個chr中間位置；
- 背景色去掉，線去掉等
- 去掉點和X軸之間的 “gap” （很多地方可用）
- 添加閾值線

2 繪制加強版Manhattan圖

1) 准備X軸標簽位置--在每條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))) +
#設置點的大小，透明度
   geom_point( aes(color=as.factor(CHR)), alpha=0.8, size=1.3) +
#設置顏色
   scale_color_manual(values = rep(c("grey", "skyblue"), 22 )) +
#設定X軸
   scale_x_continuous( label = X_axis$CHR, breaks= X_axis$center ) +
#去除繪圖區和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")) + 
#設置主題
   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()
  )

這時候是不是就可以了，🤭。

當然了既然是ggplot2繪制的Manhattan圖（點圖），那么關於點，線，坐標，主題的設置當然都可以設置了，看這里

ggplot2|詳解八大基本繪圖要素

ggplot2|theme主題設置，詳解繪圖優化-“精雕細琢”

3 玩轉Manhattan圖

1) 利用數據集自帶的`snpsOfInterest`標示顯著的位點，展示重要的基因信息

library(ggrepel)
#准備數據
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結果數據是Gene的話，label更改即可標示基因。

2) 自定義重要的基因，標示

如果有某些“目的基因”，想查看這些基因的P值呢？

新加gene和gene_annotate列即可！

#准備數據,使用基礎函數
data <- Snp_pos
#根據目的基因的位置，新加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）區域放大展示

重點展示某一區域的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|繪制區域輪廓-區域放大-尋找你的“onepiece”

4）plotly 交互展示

library(plotly)
 data <- Snp_pos %>%
 mutate( is_highlight=ifelse(SNP %in% snpsOfInterest, "yes", "no")) %>% filter(-log10(P)>0.5) #過濾一些點，交互式壓力小
# 准備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")