ggplot2（r包）绘制基因棒棒糖图

生信技能树

发布于 2025-02-05 05:42:03

18201

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：1

代码可运行

前面我们已经学习了四个包来绘制展示基因突变信息的棒棒图，其实，ggplot2也可以绘制，见资源：https://stackoverflow.com/questions/77473777/adding-branches-to-ggplot-mutation-lollipop-plot

前面已经介绍的四个软件：

数据准备

这里制作了四个位点突变新的示例数据：

rm(list=ls())
library(ggplot2)
library(ggrepel)

mut.df <- data.frame("AA" = c(201, 203, 500, 601), 
                     "Mut" = c("V201L", "R203H", "Q500*", "P601fs"), 
                     "Type" = c("Missense", "Missense", "Nonsense", "Frameshift"), 
                     "Freq" = c(2,3,4,1))

mut.df
# AA    Mut       Type Freq
# 1 201  V201L   Missense    2
# 2 203  R203H   Missense    3
# 3 500  Q500*   Nonsense    4
# 4 601 P601fs Frameshift    1

domain.df <- data.frame("Feature" = c("Start", "Dom1", "Dom2", "End"), 
                        "Type" = c("str", "dom", "dom", "str"),
                        "Start" = c(1, 180, 480, 650), 
                        "End" = c(1, 230, 630, 650))
                        
domain.df

# Feature Type Start End
# 1   Start  str     1   1
# 2    Dom1  dom   180 230
# 3    Dom2  dom   480 630
# 4     End  str   650 650

str.fill <- "#E1E1E1"
str.col <- "#16161D"

dom.fill <- c("Dom2" = "#FF0000", "Dom1" = "#AA4BAB")
dom.col <- c("#16161D")

ggplot2 绘制

1、使用`geom_rect`函数绘制边框

## 绘图
# 绘制边框
gp <- ggplot() +
  geom_rect(data = subset(domain.df, Type == "str"),
            mapping = aes(xmin = min(Start), xmax = max(End), ymin = 0.3, ymax = 0.7),
            fill = str.fill,
            colour = str.col)
gp

2、将上边绘制的边框压缩成一个长条形

# 添加y轴范围，刻度，将上边绘制的边框压缩成一个长条形
gp <- gp + scale_y_continuous(limits = c(0, 10), breaks = 0:10)
gp

3、添加棒棒图

使用geom_segment添加棒棒图的棒子，geom_point添加棒棒图上面的圈圈，geom_text_repel添加对应的文字

# 添加棒棒图
gp <- gp + geom_segment(data = mut.df, 
                        mapping = aes(x = AA, xend = AA, y = 0.7, yend = Freq)) +
  geom_point(data = mut.df,
             mapping = aes(x = AA, y = Freq, fill = Type),
             shape = 21,
             size = 2) +
  geom_text_repel(data = mut.df,
                  mapping = aes(x = AA, y = Freq, label = Mut),
                  bg.colour = "white",
                  seed = 12345,
                  nudge_y = 0.25)

gp

4、添加结构区域

再使用geom_rect添加突变区域：

# 添加结构区域
gp <- gp + geom_rect(data = subset(domain.df, Type == "dom"),
                     mapping = aes(xmin = Start, xmax = End, ymin = 0.2, ymax = 0.8, fill = Feature, group = Feature),
                     fill = dom.fill[subset(domain.df, Type == "dom")$Feature],
                     colour = dom.col)
gp

5、修改主题

优化一下配色，主题：

# 修改主题
gp <- gp +
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(linetype = "dotted")) +
  labs(x = "AA", y = "Freq", fill = "Mutation")

gp

6、再优化：将两个重叠的棒子分开不重叠

整理一下数据变成适合的数据格式：

###############################
library(tidyverse)
library(tidygraph)
library(ggraph)

data <- mut.df %>%
  select(Mut, AA, Type, Freq) %>%
  mutate(Base = paste0(Mut, '_0'), .after = 'Mut') %>%
  as_tbl_graph() %>%
  mutate(AA = rep(mut.df$AA, 2),
         Freq = c(mut.df$Freq, rep(0.5, nrow(mut.df))),
         Type = c(mut.df$Type, rep(NA, nrow(mut.df)))) %>%
  mutate(dist = sapply(AA, \(x) min(abs(x - mut.df$AA[!mut.df$AA %in% x])))) %>%
  mutate(AA = ifelse(!is.na(Type) & dist < 20, 
                     runif(n(), -50, 50), 0) + AA) 

data

绘图：不重叠的棒棒使用geom_edge_elbow函数

data %>%
  ggraph(layout = 'manual', x = AA, y = Freq) + # 绘制一个空白图
  # 绘制边框
  geom_rect(data = subset(domain.df, Type == "str"),
            mapping = aes(xmin = min(Start), xmax = max(End), 
                          ymin = 0.3, ymax = 0.7),
            fill = str.fill,
            colour = str.col) + 
  # 绘制突变结构域
  geom_rect(data = subset(domain.df, Type == "dom"),
            mapping = aes(xmin = Start, xmax = End, ymin = 0.2, 
                          ymax = 0.8, fill = Feature, group = Feature),
            fill = dom.fill[subset(domain.df, Type == "dom")$Feature],
            colour = dom.col) +
  # 添加不重叠的棒棒
  geom_edge_elbow(aes(direction = 1), strength = 0.5) +
  # 添加棒棒上面的圈圈
  geom_node_point(shape = 21, aes(fill = Type, size = Type)) +
  # 添加棒棒上的突变信息
  geom_node_text(aes(label = ifelse(is.na(Type), '', name)), 
                 angle = 90, hjust = -0.3) +
  # 优化圈圈的大小
  scale_size_manual(values = rep(3, 3), breaks = unique(mut.df$Type),
                    guide = 'none') +
  scale_fill_discrete(breaks = unique(mut.df$Type)) +
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(linetype = "dotted")) +
  labs(x = "AA", y = "Freq", fill = "Mutation") +
  ylim(c(0, 5))

结果如下：

7、再优化：两个棒棒之间添加一下空格

首先，作者写了一个函数增加两个重叠的棒棒间的空格：

# function to shift the x-axis coordinates when points are too close
shift.lollipop.x <- function(mut.pos = NULL, total.length = NULL, shift.factor = 0.05){
  
  pos.dif <- 0
  for (i in 1:length(mut.pos)){
    pos.dif <- c(pos.dif, mut.pos[i+1] - mut.pos[i])
  }
  
  idx <- which(pos.dif < shift.factor*total.length)
  
  ## deal with odd and even sets of points
  if (median(idx) %% 1==0){
    mut.pos[idx[idx < median(idx)]] <- mut.pos[idx[idx < median(idx)]] - shift.factor*total.length
    mut.pos[idx[idx > median(idx)]] <- mut.pos[idx[idx > median(idx)]] + shift.factor*total.length
  } else {
    mut.pos[idx[idx == median(idx)-0.5]] <- mut.pos[idx[idx == median(idx)-0.5]] - 0.5*shift.factor*total.length
    mut.pos[idx[idx == median(idx)+0.5]] <- mut.pos[idx[idx == median(idx)+0.5]] + 0.5*shift.factor*total.length
    
    mut.pos[idx[idx < median(idx)-0.5]] <- mut.pos[idx[idx < median(idx)-0.5]] - shift.factor*total.length
    mut.pos[idx[idx > median(idx)+0.5]] <- mut.pos[idx[idx > median(idx)+0.5]] + shift.factor*total.length
  }
  
  mut.pos
  
}

# function to split the segment into 3 parts
shift.lollipop.y <- function(x, start.y = 0.7){
  mod.start <- x - start.y
  
  set1 <- start.y + mod.start/3
  set2 <- set1 + mod.start/3
  
  as.data.frame(cbind(set1,set2))
}

使用上面定义的函数修改数据并绘图：

# 修改数据
mut.df$Shift.AA <- shift.lollipop.x(mut.df$AA, 650)
mut.df <- cbind(mut.df, shift.lollipop.y(mut.df$Freq, 0.7))

str.fill <- "#E1E1E1"
str.col <- "#16161D"
dom.fill <- c("Dom2" = "#FF0000", "Dom1" = "#AA4BAB")
dom.col <- c("#16161D")

# 绘图
gp <- ggplot() +
  geom_rect(data = subset(domain.df, Type == "str"),
            mapping = aes(xmin = min(Start), xmax = max(End), ymin = 0.3, ymax = 0.7),
            fill = str.fill,
            colour = str.col) + 
  scale_y_continuous(limits = c(0, 10), breaks = 0:10) + 
  geom_segment(data = mut.df, 
               mapping = aes(x = AA, xend = AA, y = 0.7, yend = set1)) +
  geom_segment(data = mut.df, 
               mapping = aes(x = AA, xend = Shift.AA, y = set1, yend = set2)) +
  geom_segment(data = mut.df, 
               mapping = aes(x = Shift.AA, xend = Shift.AA, y = set2, yend = Freq)) +
  geom_point(data = mut.df,
             mapping = aes(x = Shift.AA, y = Freq, fill = Type),
             shape = 21,
             size = 2) +
  geom_text_repel(data = mut.df,
                  mapping = aes(x = Shift.AA, y = Freq, label = Mut),
                  bg.colour = "white",
                  seed = 12345,
                  nudge_y = 0.25,
                  angle = 90) + 
  geom_rect(data = subset(domain.df, Type == "dom"),
                     mapping = aes(xmin = Start, xmax = End, ymin = 0.2, ymax = 0.8, fill = Feature, group = Feature),
                     fill = dom.fill[subset(domain.df, Type == "dom")$Feature],
                     colour = dom.col)

gp <- gp +
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(linetype = "dotted")) +
  labs(x = "AA", y = "Freq", fill = "Mutation")

gp