如何优雅的可视化基因集富集分析

library(ggplot2)
library(readr)
library(dplyr)
library(RColorBrewer)

# ============================================================
# 0. 手动设置区域
# ============================================================

pathway_file <- "example_pathway.tsv"      # 单个 pathway 的 GSEA 结果文件

pathway_name <- "Necroptosis"
NES_value <- -1.353
p_value <- 0

# pathway 文件中用于绘图的列
# 你之前用的是 select(4:6)，分别对应 rank, metric, score
pathway_rank_col <- 4
pathway_metric_col <- 5
pathway_score_col <- 6

# 完整 ranked gene list 的总基因数
# 只用 pathway 文件时，必须手动提供这个值
# 如果不知道，可以先用 max(pathway rank) 近似，但不推荐
n_total_gene <- 40727

# ============================================================
# 1. 读取 pathway 文件：用于 ES 曲线和 hit 竖线
# pathway_df: 包含rank,metric,score的三列数据，核心数据框
# es_df: 用于绘制折线图
# hit_df: 用于绘制竖线图
# color_block_df: 用于绘制底部colorbar
# ============================================================

pathway_df <- read_tsv(pathway_file, show_col_types = FALSE)

pathway_df <- pathway_df %>%
  select(
    rank = all_of(pathway_rank_col),
    metric_hit = all_of(pathway_metric_col),
    score = all_of(pathway_score_col)
  ) %>%
  mutate(
    rank = as.numeric(rank),
    metric_hit = as.numeric(metric_hit),
    score = as.numeric(score)
  ) %>%
  filter(
    !is.na(rank),
    !is.na(score)
  ) %>%
  arrange(rank)

es_df <- pathway_df %>%
  transmute(
    x = rank,
    ES = score
  )

hit_df <- pathway_df %>%
  transmute(
    x = rank
  )

# ============================================================
# 2. enrichplot::gseaplot2() 风格的底部色块逻辑
#    不使用所有基因 metric，只使用 pathway hit rank
# ============================================================

n <- n_total_gene

position <- rep(0, n)
hit_rank <- pathway_df$rank
hit_rank <- hit_rank[hit_rank >= 1 & hit_rank <= n]
position[hit_rank] <- 1

# enrichplot::gseaplot2() 风格色块逻辑
n_col <- 10

v <- seq(1, sum(position), length.out = n_col)

inv <- findInterval(
  rev(cumsum(position)),
  v,
  all.inside = TRUE
)

col <- c(
  rev(RColorBrewer::brewer.pal(5, "Blues")),
  RColorBrewer::brewer.pal(5, "Reds")
)

block_id <- unique(inv)

xmin <- which(!duplicated(inv))
xmax <- xmin + as.numeric(table(inv)[as.character(block_id)])

color_block_df <- tibble(
  xmin = xmin,
  xmax = xmax,
  col = col[block_id]
) %>%
  mutate(
    xmin = pmax(xmin, 1),
    xmax = pmin(xmax, n)
  )

# 防止 xmax 超出 n
color_block_df <- color_block_df %>%
  mutate(
    xmin = pmax(xmin, 1),
    xmax = pmin(xmax, n)
  )

# ============================================================
# 3. 自动设置画图范围和标注位置
# ============================================================

es_min <- min(es_df$ES, na.rm = TRUE)
es_max <- max(es_df$ES, na.rm = TRUE)

# ES 区域下方额外留空间，避免 ES 曲线和 hit 竖线重叠
extra_bottom <- 0.13

y_min <- es_min - extra_bottom
y_max <- es_max + 0.05

# 底部 colorbar：加厚
bar_ymin <- y_min + 0.010
bar_ymax <- y_min + 0.035

# hit 竖线：放在 colorbar 上方，但不要碰到 ES 曲线
hit_ymin <- bar_ymax + 0.008
hit_ymax <- hit_ymin + 0.070

# 在 ES 区域和 hit 区域之间加一条横线，模仿第二张图
sep_y <- hit_ymax + 0.008

# NES 和 P value 的默认位置
nes_x <- n * 0.07
nes_y <- es_min + (es_max - es_min) * 0.18

p_x <- n * 0.60
p_y <- es_max - (es_max - es_min) * 0.18

# y 轴刻度只显示 ES 曲线附近的刻度，不显示底部装饰区域
y_breaks <- pretty(c(es_min, es_max), n = 4)
y_breaks <- y_breaks[y_breaks >= es_min & y_breaks <= es_max]

# ============================================================
# 4. 作图
# ============================================================

ggplot() +
  
  # enrichplot::gseaplot2() 风格底部色块，加厚
  geom_rect(
    data = color_block_df,
    aes(
      xmin = xmin,
      xmax = xmax,
      ymin = bar_ymin,
      ymax = bar_ymax
    ),
    fill = color_block_df$col,
    color = NA
  ) +
  
  # hit 竖线：独立放在 ES 曲线下方
  geom_segment(
    data = hit_df,
    aes(
      x = x,
      xend = x,
      y = hit_ymin,
      yend = hit_ymax
    ),
    linewidth = 0.35,
    color = "black"
  ) +
  
  # ES 区域和 hit 区域之间的分隔线
  geom_hline(
    yintercept = sep_y,
    linewidth = 1.6,
    color = "black"
  ) +
  
  # ES running score 曲线
  geom_line(
    data = es_df,
    aes(x = x, y = ES),
    color = "#8a2be2",
    linewidth = 0.75
  ) +
  
  annotate(
    "text",
    x = nes_x,
    y = nes_y,
    label = paste0("NES = ", NES_value),
    hjust = 0,
    vjust = 0,
    size = 7.5,
    color = "black"
  ) +
  
  annotate(
    "text",
    x = p_x,
    y = p_y,
    label = paste0("italic(P) == ", p_value),
    hjust = 0,
    size = 6.6,
    color = "black",
    parse = TRUE
  ) +
  
  scale_x_continuous(
    limits = c(1, n),
    expand = c(0, 0),
    breaks = pretty_breaks(n = 8)
  ) +
  
  scale_y_continuous(
    limits = c(y_min, y_max),
    breaks = y_breaks,
    labels = function(x) sprintf("%.1f", x),
    expand = c(0, 0)
  ) +
  
  labs(
    title = pathway_name,
    x = NULL,
    y = NULL
  ) +
  
  theme_classic(base_size = 24) +
  theme(
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    
    plot.title = element_text(
      hjust = 0.5,
      size = 34,
      face = "plain",
      margin = margin(b = 2)
    ),
    
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    
    axis.text.y = element_text(size = 28, color = "black"),
    axis.ticks.y = element_line(linewidth = 1.1, color = "black"),
    axis.ticks.length = unit(0.18, "cm"),
    
    panel.border = element_rect(
      color = "black",
      fill = NA,
      linewidth = 2.2
    ),
    axis.line = element_blank(),
    
    panel.grid.major.x = element_line(
      color = "grey88",
      linewidth = 0.5
    ),
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank(),
    
    plot.margin = margin(5, 8, 5, 8)
  )

ggsave("GSEA_Necroptosis.pdf",width=6.61,height = 2.85)