GSEA 结果展示

# ============================================================
# GSEA 结果展示新函数：gene rank 泳道 + p.adjust / NES 展示
# ------------------------------------------------------------
# 来源：KS 原创函数合集（2024-7-12 教程）
#   KS代码/KS-account-codes/5-KS原创函数合集/
#   "2024-7-12 GSEA结果展示新函数-gene rank-p.adjust-NES值展示"
# 整理：2026-08-25，仅移除原教程的 setwd() 硬编码路径，
#       ks_plot_gsea 函数体保持原样未改动。
#
# 输入：clusterProfiler::GSEA() 返回的 gseaResult 对象
# 输出：三联图 = gene rank 排名泳道 + p.adjust 面板 + NES 面板
#   labels_size = T : p.adjust 与 NES 以气泡形式展示
#                     （颜色映射 -log10(p.adjust) / NES，大小映射 setSize）
#   labels_size = F : p.adjust（科学计数法）与 NES（4 位小数）以数字文本展示
# ============================================================

# ---- 使用示例（取消注释即可运行）----
# 需要本机的 res.RData（gseaResult 对象，仅存于原始素材目录，不入库）：
#   load("res.RData")
#
# labels_size = T：气泡模式（效果见 1.png）
# ks_plot_gsea(data = res,
#              term = res@result$Description[1:10],
#              rankColor = "#90191B",
#              labels_size = T)
#
# labels_size = F：数字模式（效果见 2.png）
# ks_plot_gsea(data = res,
#              term = res@result$Description[1:10],
#              rankColor = "#90191B",
#              labels_size = F)
#
# term 也可以手动指定通路名称子集，例如：
# ks_plot_gsea(data = res,
#              term = c("immunoglobulin production",
#                       "B cell receptor signaling pathway",
#                       "T cell cytokine production"),
#              rankColor = "#90191B",
#              labels_size = T)

ks_plot_gsea <- function(data,#clusterProfiler分析GSEA得到的结果文件，Large gseaResult
                         term,#需要展示的terms
                         rankColor,#线的颜色
                         labels_size=T#选择为T，展示NES和pvalue使用气泡形式。否则使用数字展示
                         ){

  require(ggplot2)
  require(patchwork)
  require(clusterProfiler)
  require(viridis)
  require(grDevices)

  #clusterprofile
  gseaScores <- getFromNamespace("gseaScores", "DOSE")

  # define function
  gsInfo <- function(object, terms) {
    geneList <- object@geneList

    gsea_result_table <- object@result
    site=which(gsea_result_table$Description==terms, arr.ind = TRUE)
    genesetid=gsea_result_table$ID[site]

    geneSetID <- object@result[genesetid, "ID"]

    geneSet <- object@geneSets[[geneSetID]]
    exponent <- object@params[["exponent"]]
    df <- gseaScores(geneList, geneSet, exponent, fortify=TRUE)
    df$ymin <- 0
    df$ymax <- 0
    pos <- df$position == 1
    h <- diff(range(df$runningScore))/20
    df$ymin[pos] <- -h
    df$ymax[pos] <- h
    df$geneList <- geneList

    df$Description <- object@result[geneSetID, "Description"]
    return(df)
  }




  if (length(term) == 1) {

    gsea_plotData <- gsInfo(data, terms = term)

  } else {

    gsea_plotData <- do.call(rbind, lapply(term, gsInfo, object = data))

  }

  #get gsea rank data

  gsea_plotData <- gsea_plotData %>%
    mutate("gene_name" = rep(data@gene2Symbol, length(term))) %>%
    filter(position == 1)

  test = gsea_plotData


  #process gsea rank data for plot
  test$xend = test$x+1
  test$xend2 = ''
  test$xend3 = ''

  for (i in 1:nrow(test)) {
    test$xend2[i] = test$xend[i+1] - test$xend[i]
  }

  for (i in 1:ncol(test)){
    test[,i][is.na(test[,i])] <- 5
  }

  test$xend2 <- as.numeric(test$xend2)

  for (i in 1:nrow(test)) {
    test$xend3[i] = (test$x[i] + test$xend2[i])
  }

  test$xend3 <- as.numeric(test$xend3)


  test$y = rep(1:length(term), as.numeric(table(test$Description)))


  label_breaks <- seq(1:length(term))+0.5

  p=ggplot(test) +
    # geom_rect(aes(xmin = x-1,xmax = length(data@gene2Symbol), ymin = y-0.01 , ymax = y+0.8, fill=x), lwd=4)+
    # scale_fill_gradientn(colours = colorRampPalette(c("#90191B","white","#3F90C9"))(100))+
    geom_rect(aes(xmin = x,xmax = xend, ymin = y-0.01, ymax = y+0.8, color=x), color=rankColor, size=0.5)+
    scale_y_continuous(expand = c(0,0),labels=term, breaks = label_breaks)+
    scale_x_continuous(expand = c(0,0))+
    theme_bw() +
    theme(panel.grid.major=element_blank(),
          panel.grid.minor=element_blank(),
          axis.text.y = element_text(size = 12, color = "black"),
          axis.title.y = element_text(size = 15,colour = 'black'),
          axis.ticks.x = element_blank(),
          legend.position = "none")




  #NES和adjust Pvalue
  gesa_table <- data@result
  gsea_table = gesa_table[gesa_table$Description %in% term, ]


  if(labels_size==T){
    p1 = ggplot(gsea_table, aes(x=1, y=label_breaks))+
      geom_point(aes(fill=-log(p.adjust),size=setSize), shape=21)+
      scale_fill_gradientn(colors=viridis(256))+
      theme(panel.grid.major=element_blank(),
            panel.grid.minor=element_blank(),
            panel.background = element_blank(),
            axis.text = element_blank(),
            axis.title= element_blank(),
            axis.ticks = element_blank(),
            axis.line = element_blank(),
            legend.key.height = unit(0.5, "cm"))+
      scale_size(limits = c(0,400), breaks = c(100, 200, 300))




    p2 = ggplot(gsea_table, aes(x=1, y=label_breaks))+
      geom_point(aes(fill=NES,size=setSize), shape=21)+
      scale_fill_gradientn(colours = colorRampPalette(c("#3F90C9","white","#90191B"))(100))+
      theme(panel.grid.major=element_blank(),
            panel.grid.minor=element_blank(),
            panel.background = element_blank(),
            axis.text = element_blank(),
            axis.title= element_blank(),
            axis.ticks = element_blank(),
            axis.line = element_blank(),
            legend.key.height = unit(0.5, "cm"))+
      scale_size(limits = c(0,400), breaks = c(100, 200, 300))

  }else{


    format_number_in_scientific_notation <- function(number) {
      # 使用sprintf函数格式化数字为科学计数法
      formatted_number <- sprintf("%e", number)
      return(formatted_number)
    }


    p1 = ggplot(gsea_table, aes(x=1, y=label_breaks))+
      geom_text(aes(label=format_number_in_scientific_notation(p.adjust)))+
      theme(panel.grid.major=element_blank(),
            panel.grid.minor=element_blank(),
            panel.background = element_blank(),
            axis.text = element_blank(),
            axis.title= element_blank(),
            axis.ticks = element_blank(),
            axis.line = element_blank())+
      ggtitle("p.adjust")+
      theme(plot.title = element_text(hjust = 0.5))



    p2 = ggplot(gsea_table, aes(x=1, y=label_breaks))+
      geom_text(aes(label=round(NES,4)))+
      theme(panel.grid.major=element_blank(),
            panel.grid.minor=element_blank(),
            panel.background = element_blank(),
            axis.text = element_blank(),
            axis.title= element_blank(),
            axis.ticks = element_blank(),
            axis.line = element_blank())+
      ggtitle("NES")+
      theme(plot.title = element_text(hjust = 0.5))



  }


  # library(patchwork)

  if(labels_size==T){

    plot_fig <- p+p1+p2+plot_layout(ncol = 3,width=c(2,0.2,0.2),guides = 'collect')&
      theme(legend.position='bottom')

  }else{

    plot_fig <- p+p1+p2+plot_layout(ncol = 3,width=c(2,0.5,0.5))

  }



  return(plot_fig)

}

# ============================================================
# GSEA 渐变峰峦图 + NES 气泡图（"日照金山"配色）
# ------------------------------------------------------------
# 来源：微信公众号「科易璟」教程《R绘图：日照金山，峰峦图升级版》
#       https://mp.weixin.qq.com/s/tfy8T2pz1PQ_C-HDFMkgVQ
# 整理：2026-09-15。相对原教程做了三处调整（见文末"与原教程的差异"）：
#       1) 数据源改为从本目录的 res.RData（gseaResult 对象）自动重建，
#          去掉原教程的 E:/ 硬编码路径；
#       2) 补上 patchwork 自动拼图，替代原教程用 AI 手工拼图那一步；
#       3) 全程输出 PDF（原教程也是 ggsave 出 PDF）。
#
# 图的形式（对应 参考图_文章原图/04_成品_日照金山.png）：
#   左：峰峦图 —— x = Log2FoldChange，y = 通路，填充映射 -log10(pvalue)
#       渐变 紫→白→黄→金，即所谓"日照金山"
#   右：NES 气泡图 —— 点的颜色与大小都映射 NES，蓝→白→红发散配色
#   两侧并排后，一张图同时交代基因表达方向（左）与富集方向/强度（右）
#
# 为什么要"分别绘制再拼"：原教程明确说"在一张图上绘制时，总是出现配色覆盖"，
#   因为峰峦图要按 -log10(p) 填色、气泡图要按 NES 配色，两者色标不同，
#   共用一个 scale_fill / scale_color 会互相覆盖。拆开画各自配色，最后拼起来即可。
# ============================================================

# ---------------------------------------------------------------- 0. 依赖
suppressMessages({
  library(tidyr)       # separate_rows()
  library(dplyr)       # 管道与数据处理
  library(ggplot2)
  library(ggridges)    # geom_density_ridges()
  library(patchwork)   # 拼图（替代 AI 手工拼）
  library(clusterProfiler)  # 读取 gseaResult 对象
})
# 可选：library(openxlsx)  # 若要把结果另存为 xlsx

# ---------------------------------------------------------------- 0b. 定位脚本自身所在目录
# 这样无论你从哪个工作目录调用，输入数据和输出图都落在脚本旁边，
# 不会因为 cwd 不同而把 result/ 建到别处。
#   - Rscript 调用：能从 commandArgs 拿到脚本路径，最可靠
#   - RStudio 里 source()：拿不到，退回当前工作目录（此时请先 setwd 到本目录）
script_dir <- local({
  a <- grep("^--file=", commandArgs(trailingOnly = FALSE), value = TRUE)
  if (length(a)) dirname(normalizePath(sub("^--file=", "", a[1]))) else getwd()
})
out_dir <- file.path(script_dir, "result")
if (!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE)

# ---------------------------------------------------------------- 1. 读入数据
# 本目录（或上一级 06_GSEA/）的 res.RData 是一个 clusterProfiler::GSEA() 的
# gseaResult 对象。没有的话见文末"用你自己的数据"。
find_res <- function() {
  rel <- c("res.RData", "../res.RData", "../../res.RData")
  # 先按脚本所在目录找，再按当前工作目录找
  cands <- c(file.path(script_dir, rel), rel)
  hit <- cands[file.exists(cands)]
  if (!length(hit)) {
    stop("没找到 res.RData（已在脚本目录与当前目录下查找）。\n",
         "  脚本目录: ", script_dir, "\n",
         "  当前目录: ", getwd(), "\n",
         "请把 GSEA 结果放到上述任一目录，或参考文末『用你自己的数据』改用 txt 输入。")
  }
  normalizePath(hit[1])
}
res_file <- find_res()
# 注意：load() 会把变量 res 注入当前环境（source 时即全局环境）。
# 这是刻意的——原教程和本脚本后续都直接用 `res`。若不想污染全局，
# 可改为 e <- new.env(); load(res_file, envir = e); res <- e$res
load(res_file)
cat("[数据] ", res_file, " —— ", nrow(res@result), " 条通路\n", sep = "")

# ---------------------------------------------------------------- 2. 构造绘图用的两张表
# 原教程需要两份输入：
#   genelist     : 两列 gene / logfc（每个基因的 log2FoldChange）
#   KEGG_gsearesult : GSEA 结果表（含 Description / pvalue / NES / core_enrichment）
# 这两份都能从 gseaResult 对象里还原：
#   logfc      <- res@geneList   （names 是 ENTREZ ID，值是 logFC）
#   gene symbol<- res@gene2Symbol（ENTREZ ID -> Symbol）
genelist <- data.frame(
  gene  = unname(res@gene2Symbol[names(res@geneList)]),
  logfc = as.numeric(res@geneList),
  stringsAsFactors = FALSE
) %>%
  dplyr::filter(!is.na(gene)) %>%
  dplyr::distinct(gene, .keep_all = TRUE)   # 一个 symbol 对应多个 ENTREZ 时只留第一个
cat("[genelist] ", nrow(genelist), " 个基因\n", sep = "")

# gsearesult1 用于气泡图：取最显著的 20 条通路
gsearesult1 <- res@result %>%
  dplyr::arrange(pvalue) %>%
  head(20) %>%
  dplyr::mutate(Description = factor(Description, levels = rev(unique(Description))))

# gsearesult 用于峰峦图：把 core_enrichment（"GENE/GENE/..." 长字符串）
# 拆成一行一个基因，再 join 上各自的 logfc，才能画密度曲线
gsearesult <- gsearesult1 %>%
  dplyr::mutate(log10P = -log10(pvalue)) %>%
  separate_rows(core_enrichment, sep = "/") %>%
  dplyr::left_join(genelist, by = c("core_enrichment" = "gene")) %>%
  dplyr::mutate(Description = factor(Description, levels = rev(unique(Description))))

na_n <- sum(is.na(gsearesult$logfc))
cat("[gsearesult] 长表 ", nrow(gsearesult), " 行", sep = "")
if (na_n > 0) cat("（其中 ", na_n, " 行没匹配到 logfc，会被 ggplot 自动丢弃）", sep = "")
cat("\n")

# ---------------------------------------------------------------- 3. 峰峦图
# 配色解析：#8075ad 紫 → #f5edf0 近白 → #f5da73 黄 → #ffb800 金
custom_colors <- colorRampPalette(c("#8075ad", "#f5edf0", "#f5da73", "#ffb800"))(100)

p <- ggplot(gsearesult, aes(x = logfc, y = Description, fill = log10P)) +
  geom_density_ridges(alpha = 0.7, scale = 1.5) +
  labs(x = "Log2FoldChange", y = "") +
  scale_fill_gradientn(colors = custom_colors, name = "-Log10(pvalue)") +
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA, linewidth = 1),
    axis.text = element_text(color = "black", size = 12),
    axis.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    legend.title = element_text(size = 12)
  )
ggsave(file.path(out_dir, "ridge_kegg.pdf"), plot = p, height = 6.5, width = 7.8)

# ---------------------------------------------------------------- 4. NES 气泡图
# 配色解析：#008bd0 蓝 → #eeeeee 白 → #fc4e00 红（发散色，0 附近为白）
custom_colors1 <- colorRampPalette(c("#008bd0", "#eeeeee", "#fc4e00"))(100)

p1 <- ggplot(gsearesult1, aes(x = 1, y = Description)) +
  geom_point(aes(size = abs(NES), color = NES)) +        # 大小和颜色都映射 NES
  scale_color_gradientn(colors = custom_colors1, name = "NES", limits = c(-3, 3)) +
  scale_size_continuous(range = c(2, 6), guide = "none") +  # 见下方"关于尺寸图例"
  labs(x = NULL, y = NULL) +
  theme_minimal() +
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    panel.border = element_blank(),
    axis.ticks = element_blank(),
    axis.text.x = element_blank()
  )
ggsave(file.path(out_dir, "ridge_point.pdf"), p1, width = 5, height = 6)

# ---------------------------------------------------------------- 5. 拼图（原教程用 AI 手工做）
# 原教程最后一步是在 Adobe Illustrator 里手动对齐两张图；这里用 patchwork 自动拼。
#
# 关键：拼图顺序与宽度比例。原文成品的构图是
#       [ 通路名 | NES 气泡列 | 峰峦密度曲线 | 图例 ]
#   通路名只出现一次（它是峰峦图的 y 轴标签，位于最左侧），气泡列紧贴在峰峦图左边。
#   所以这里要把**气泡图放在左边并去掉它的 y 轴文字**，
#   否则通路名会被打印两遍，且 patchwork 的 widths 是"面板比例"而非绝对宽度——
#   若把带长标签的峰峦图放左边再随手给个 4:1，峰峦面板会被挤成一条窄缝
#   （实测：绘图面板从单图的 129pt 缩到 32pt，山脊全糊在一起）。
p_compact <- p +
  theme(
    axis.text.y = element_blank(),                     # 通路名交给左侧的气泡图显示
    plot.margin = ggplot2::margin(5.5, 0, 5.5, 0)      # 收掉右侧留白，让山脊贴紧气泡
  )

# 顺序：气泡图在左（带通路名）、峰峦图在右（不带），即 [通路名 | 气泡 | 山脊]
combined <- p1 + p_compact + plot_layout(widths = c(1, 12))
ggsave(file.path(out_dir, "ridge_combined.pdf"), combined, width = 13, height = 6.5)

cat("\n[产物] ", out_dir, " 下：ridge_kegg.pdf / ridge_point.pdf / ridge_combined.pdf\n", sep = "")
cat("       成品效果对照：参考图_文章原图/04_成品_日照金山.png\n")
cat("       （combined 是 patchwork 自动拼；构图与原文 AI 手工成品一致，但细节可能有微小差异）\n")

# ============================================================
# 与原教程的差异
# ------------------------------------------------------------
# 1. 数据源：
#    原教程  genelist     <- read.table("E:/.../genelist.txt", ...)
#            gsearesult   <- read.table("E:/.../KEGG_gsearesult.txt", ...)
#    本脚本  从 res.RData（gseaResult）自动重建这两份表，不需要外部 txt。
#    若你手上正好有那两份 txt，直接改成：
#      genelist   <- read.table("genelist.txt", sep="\t", header=TRUE)
#      gsearesult <- read.table("KEGG_gsearesult.txt", sep="\t", header=TRUE)
#    其余处理步骤完全一致（genelist 的 gene 列需为字符型，见原教程的 mutate）。
#
# 2. 拼图：原教程用 AI 手工拼；本脚本用 patchwork 自动拼。
#    想完全复刻手工版，就只跑第 3、4 节导出两张图，再自行拼。
#
# 3. 配色与主题：与原教程逐字一致，未做改动。
#
# 注意：原教程第 3 节注释里写"# 使用 iris 数据集，根据不同物种绘制峰峦图"——
#       那是从别处复制代码时残留的注释，与本文数据无关，已删除。
# ============================================================