volcanoGSEA

Create publication-ready volcano plots with integrated GSEA-style enrichment barcodes, plus nested DMR donut summaries for hypo/hyper methylation contrasts.

创建带有GSEA富集条形码的出版级火山图。

Inspiration / 灵感来源

This package was inspired by the volcano plot visualization in:

Huang, H., et al. (2024). Vitamin C enhances TET activity and promotes DNA demethylation in mouse embryonic stem cells. Science. DOI: 10.1126/science.adk0775

本R包的火山图设计灵感来源于上述Science论文中的可视化方案。

Features / 功能特点

• Background shading for significant regions (up/down regulated) / 显著区域背景着色
• Significance boxes with up/down gene counts / 显著性框显示上调/下调基因数量
• Multiple pathway barcodes at bottom (each pathway in separate row) / 底部多通路条形码
• Shared gene detection (diamond markers for genes in multiple pathways) / 共享基因检测（菱形标记）
• ES and p-value display for each pathway / 每个通路显示ES和p值
• Group expression arrows indicating which group has higher expression / 组别表达箭头指示高表达方向
• Journal-specific color schemes (Nature, Cell, Science, NEJM, Lancet) / 期刊专用配色方案
• Nested donut plots for multi-comparison DMR summaries / 多比较DMR嵌套环形图
• Helpers to extract hypo/hyper gene sets and prioritize study-relevant terms / 提取hypo/hyper基因并筛选研究相关term

Installation / 安装

# From local directory / 本地安装
devtools::install_local("/path/to/volcanoGSEA")

# Or source directly / 或直接加载源码
source("/path/to/volcanoGSEA/R/volcano_gsea.R")
source("/path/to/volcanoGSEA/R/utils.R")

Quick Start / 快速开始

library(volcanoGSEA)

p <- volcano_gsea(

diff_data = my_diff_results,
gene_col = "gene_symbol",
logfc_col = "log2FoldChange", 
pval_col = "padj",
genesets = list(
  "Pathway_A" = pathway_a_genes,
  "Pathway_B" = pathway_b_genes
),
geneset_stats = list(
  "Pathway_A" = list(ES = "0.52", pval = "1.2e-5"),
  "Pathway_B" = list(ES = "-0.38", pval = "3.4e-3")
),
key_genes = c("TP53", "KRAS", "MYC"),
# Group expression arrows / 组别表达箭头
group_arrows = TRUE,
left_group = "Control",
right_group = "Treatment"
)

ggsave("volcano_gsea.png", p, width = 11, height = 9, dpi = 300)

Main Parameters / 主要参数

Parameter	Description	Default
diff_data	Differential expression data frame / 差异表达数据框	-
gene_col	Gene symbol column / 基因符号列名	"gene"
logfc_col	log2FC column / log2FC列名	"logFC"
pval_col	Adjusted p-value column / 校正p值列名	"FDR"
genesets	Named list of gene sets / 基因集列表	NULL
geneset_stats	GSEA stats list(ES=, pval=)	NULL
key_genes	Genes to highlight / 高亮基因	NULL
logfc_cutoff	log2FC threshold / log2FC阈值	0.5
pval_cutoff	P-value threshold / P值阈值	0.05
show_shade	Show background shading / 显示背景着色	TRUE
show_top_box	Show gene count boxes / 显示基因计数框	TRUE
colors	Color scheme / 配色方案	default_colors()
xlim	X-axis limits / X轴范围	c(-3, 3)
ylim	Y-axis limits / Y轴范围	NULL (auto)
barcode_ratio	Barcode area ratio / 条形码区域比例	0.25
group_arrows	Show group arrows / 显示组别箭头	FALSE
left_group	Left group label / 左侧组名	NULL
right_group	Right group label / 右侧组名	NULL

Complete Example / 完整示例

library(volcanoGSEA)

# Load DEG results / 加载差异表达结果
diff <- read.csv("DEG_results.csv")

# Define gene sets / 定义基因集
genesets <- list(
"Base Excision Repair" = c("Polb", "Parp1", "Ung", "Neil1", "Lig1", "Pcna"),
"Nucleotide Excision Repair" = c("Xpc", "Xpa", "Ercc1", "Lig1", "Pcna"),
"Mismatch Repair" = c("Msh2", "Msh6", "Mlh1", "Lig1", "Pcna")
)

# GSEA statistics / GSEA统计量
geneset_stats <- list(
"Base Excision Repair" = list(ES = "-0.42", pval = "2.3e-5"),
"Nucleotide Excision Repair" = list(ES = "0.31", pval = "8.7e-3"),
"Mismatch Repair" = list(ES = "0.18", pval = "0.042")
)

# Create plot / 创建图表
p <- volcano_gsea(
diff_data = diff,
gene_col = "gene_name",
logfc_col = "log2FoldChange",
pval_col = "padj",
genesets = genesets,
geneset_stats = geneset_stats,
key_genes = c("Polb", "Parp1", "Xpc", "Msh2"),
logfc_cutoff = 1,
pval_cutoff = 0.05,
xlim = c(-6, 6),
ylim = c(0, 80),
xlab = "Treatment vs Control (log2FC)",
barcode_ratio = 0.28,
# Group arrows / 组别箭头
group_arrows = TRUE,
left_group = "WT",
right_group = "KO"
)

ggsave("figure_volcano.png", p, width = 11, height = 9, dpi = 600)

Visual Layout / 可视化布局

┌─────────────────────────────────────────────────────┐
│  ┌──────┐ Down-regulated     Up-regulated ┌───────┐│
│  │Down:N│ (blue shade)       (red shade)  │ Up: N ││
│  │      │  ← WT high           KO high →  │       ││
│  │  ●   │    ● ○ ○ ○ ●                    │   ●   ││
│  └──────┘ ○ ○ ○ ○ ○ ○ ○ ○ ○               └───────┘│
│  ........α = 0.05.............................      │
│                                                     │
│  Pathway_A  ││││││││││││││         ES: 0.52        │
│  Pathway_B  │││││││││││││          ES: -0.38       │
│  Pathway_C  ││││││││││             ES: 0.18        │
└─────────────────────────────────────────────────────┘

Legend / 图例:
○ Background genes / 背景基因
● Pathway-specific genes / 通路特异基因
◇ Shared genes (in multiple pathways) / 共享基因

Helper Functions / 辅助函数

Function	Description
read_gmt()	Read MSigDB GMT file / 读取GMT文件
default_colors()	Get default color scheme / 默认配色
default_dmr_colors()	DMR donut color scheme / DMR环图配色
journal_colors(style)	Journal-specific colors / 期刊配色
volcano_gsea_deseq2()	Wrapper for DESeq2 / DESeq2封装
volcano_gsea_edger()	Wrapper for edgeR / edgeR封装
volcano_gsea_limma()	Wrapper for limma / limma封装
load_stj_dmr_inputs()	Load STJ DMR summary tables / 读取STJ DMR汇总表
extract_hypo_hyper_genes()	Build hypo/hyper gene sets from integrated tables / 从整合结果提取hypo/hyper基因
plot_dmr_nested_donuts()	Plot 4-comparison nested donut summary / 绘制4组嵌套环图
recommend_stj_gsea_terms()	Study-driven term keywords / 基于课题叙事推荐term关键词
select_stj_relevant_terms()	Filter enrichment tables to relevant terms / 从富集表筛选相关term

STJ Example / STJ项目示例

source("R/dmr_gsea.R")

dmr_inputs <- load_stj_dmr_inputs("/home/crafy/media/sdc/workspace/stj")

p <- plot_dmr_nested_donuts(
  summary_stats = dmr_inputs$summary_stats,
  regional_distribution = dmr_inputs$regional_distribution
)

ggsave("Fig1B_DMR_nested_donuts.png", p, width = 12, height = 9, dpi = 300)

gene_sets <- extract_hypo_hyper_genes(
  integrated_dir = "/home/crafy/media/sdc/workspace/stj/results",
  region = "promoter",
  annotation_file = "/home/crafy/media/sdc/workspace/stj/data/gene_annotations.csv",
  significant_only = FALSE,
  meth_diff_cutoff = 5
)

head(gene_sets$KO_vs_WT$hypo)

ko_terms <- read.csv("/home/crafy/media/sdc/workspace/stj/03.Downstream_Analysis/02_methylation_analysis/enrichment/GO_analysis/KO_vs_WT/GO_BP_hypo.csv")
select_stj_relevant_terms(ko_terms, comparison = "KO_vs_WT", direction = "Hypo")

Dependencies / 依赖

• ggplot2 (>= 3.4.0)
• ggrepel
• dplyr

License / 许可证

MIT