通路富集GSEA分析¶

步骤1：对所有基因按某个统计量排序
        → 比如按 log2FC × (-log10 pvalue) 从大到小排

步骤2：从排名第1个基因开始"走楼梯"
        → 遇到基因集成员，分数+1步（上楼）
        → 遇到非成员，分数-1步（下楼）
        → 最大偏离零点就是ES（富集分数）

步骤3：通过排列检验（permutation）估算p值
        → 随机打乱基因顺序1000次
        → 看真实ES在随机分布中有多极端

# 安装fgsea包（2026最新版）
if (!require("BiocManager")) install.packages("BiocManager")  # 安装BiocManager
BiocManager::install("fgsea")         # 安装fgsea
BiocManager::install("msigdbr")       # MSigDB基因集的R接口
BiocManager::install("org.Hs.eg.db")  # 人类注释

library(fgsea)       # 加载fgsea
library(msigdbr)     # 加载MSigDB
library(ggplot2)     # 画图
library(data.table)  # 高效数据处理

# 假设你已有DESeq2的差异分析结果
# deseq_result 包含 log2FoldChange, pvalue, padj 等列

# 方法1：用 -log10(pvalue) × sign(log2FC) 作为排名统计量
deseq_result$rank_metric <- -log10(deseq_result$pvalue) * sign(deseq_result$log2FoldChange)

# 方法2：直接用 log2FC（更简单，但信息量少）
# deseq_result$rank_metric <- deseq_result$log2FoldChange

# 方法3：用t统计量（如果有的话）
# deseq_result$rank_metric <- deseq_result$stat

# 构建命名向量（名字=基因名，值=排名统计量）
gene_ranks <- deseq_result$rank_metric          # 取排名值
names(gene_ranks) <- deseq_result$gene_symbol   # 用基因名命名
gene_ranks <- sort(gene_ranks, decreasing = TRUE)  # 从大到小排序
gene_ranks <- gene_ranks[!is.na(gene_ranks)]    # 去除NA值
gene_ranks <- gene_ranks[!duplicated(names(gene_ranks))]  # 去除重复基因

head(gene_ranks)  # 查看排名前几的基因
tail(gene_ranks)  # 查看排名后几的基因

# 获取Hallmark基因集（最常用，50个标志性通路）
hallmark_sets <- msigdbr(
  species = "Homo sapiens",   # 物种
  category = "H"              # H = Hallmark基因集
)

# 转换为fgsea需要的list格式
hallmark_list <- split(hallmark_sets$gene_symbol, hallmark_sets$gs_name)  # 按基因集名分组
length(hallmark_list)  # 查看有多少个基因集（应该是50个）
head(names(hallmark_list))  # 查看基因集名称

# 其他常用基因集类别
# C2: 策划基因集（包含KEGG、Reactome等）
c2_sets <- msigdbr(species = "Homo sapiens", category = "C2")
c2_list <- split(c2_sets$gene_symbol, c2_sets$gs_name)

# C5: GO基因集
c5_sets <- msigdbr(species = "Homo sapiens", category = "C5")
c5_list <- split(c5_sets$gene_symbol, c5_sets$gs_name)

# 运行fgsea分析（使用多层级蒙特卡洛方法）
gsea_result <- fgsea(
  pathways = hallmark_list,   # 基因集列表
  stats = gene_ranks,         # 排名统计量
  minSize = 15,               # 最小基因集大小（过滤太小的）
  maxSize = 500,              # 最大基因集大小（过滤太大的）
  eps = 0                     # 设为0可以计算任意小的p值
)

# 查看结果（按NES排序）
gsea_result <- gsea_result[order(gsea_result$NES, decreasing = TRUE), ]  # NES降序排列
head(gsea_result, 10)  # 看上调最显著的10个通路

# 结果列说明
# pathway: 基因集名称
# pval: 原始p值
# padj: BH校正后p值
# log2err: 对数误差估计
# ES: 富集分数
# NES: 标准化富集分数（跨基因集可比）
# size: 基因集中被检测到的基因数
# leadingEdge: 核心基因列表

# 筛选显著富集的通路（padj < 0.05）
sig_pathways <- gsea_result[gsea_result$padj < 0.05, ]  # 取显著的

# 分别看上调和下调的通路
up_pathways <- sig_pathways[sig_pathways$NES > 0, ]    # NES>0 = 上调富集
down_pathways <- sig_pathways[sig_pathways$NES < 0, ]  # NES<0 = 下调富集

cat("上调通路数:", nrow(up_pathways), "\n")   # 打印上调通路数
cat("下调通路数:", nrow(down_pathways), "\n")  # 打印下调通路数

# 单个通路的GSEA图（最经典的可视化）
plotEnrichment(
  hallmark_list[["HALLMARK_TNFA_SIGNALING_VIA_NFKB"]],  # 指定通路
  gene_ranks                                               # 排名统计量
) +
  labs(title = "TNFa signaling via NFkB") +  # 添加标题
  theme_classic()                              # 经典主题

# 同时展示多个通路的GSEA图（表格格式）
topPathways <- gsea_result[head(order(gsea_result$pval), 10), ]$pathway  # 取Top10通路

plotGseaTable(
  hallmark_list[topPathways],  # 通路列表
  gene_ranks,                   # 排名
  gsea_result,                  # 结果
  gseaParam = 0.5               # ES计算参数
)

# 自定义气泡图
library(ggplot2)

# 准备绘图数据
plot_data <- gsea_result[gsea_result$padj < 0.05, ]  # 取显著通路
plot_data$pathway <- gsub("HALLMARK_", "", plot_data$pathway)  # 去掉前缀
plot_data$pathway <- gsub("_", " ", plot_data$pathway)         # 下划线换空格

# 画气泡图
ggplot(plot_data, aes(x = NES, y = reorder(pathway, NES))) +  # NES为X轴
  geom_point(aes(size = size, color = padj)) +  # 大小=基因集大小，颜色=padj
  scale_color_gradient(low = "red", high = "blue") +  # 颜色梯度
  labs(x = "Normalized Enrichment Score",  # X轴标签
       y = NULL,                            # Y轴不加标签
       title = "GSEA Hallmark Pathways") +  # 标题
  theme_minimal() +                         # 简洁主题
  geom_vline(xintercept = 0, linetype = "dashed")  # 添加零线

library(clusterProfiler)  # 加载clusterProfiler

# clusterProfiler也支持GSEA（底层调用fgsea）
gsea_go <- gseGO(
  geneList = gene_ranks,       # 排名基因列表
  OrgDb = org.Hs.eg.db,        # 注释数据库
  keyType = "SYMBOL",          # 基因ID类型
  ont = "BP",                   # GO子类
  minGSSize = 15,              # 最小基因集
  maxGSSize = 500,             # 最大基因集
  pvalueCutoff = 0.05,         # p值阈值
  eps = 0                       # 精确p值
)

# GSEA-KEGG
gsea_kegg <- gseKEGG(
  geneList = gene_ranks_entrez,  # 需要Entrez ID的排名列表
  organism = "hsa",               # 物种
  minGSSize = 15,
  maxGSSize = 500,
  pvalueCutoff = 0.05,
  eps = 0
)

# clusterProfiler的GSEA可视化
gseaplot2(gsea_go, geneSetID = 1:3, pvalue_table = TRUE)  # 前3个通路
dotplot(gsea_go, showCategory = 15)  # 气泡图
ridgeplot(gsea_go, showCategory = 15)  # 山脊图

# 推荐方案（按推荐程度排序）：
# 1. DESeq2的stat列（Wald统计量）— 最推荐，信息量最大
gene_ranks <- setNames(res$stat, res$gene)

# 2. -log10(pvalue) × sign(log2FC) — 通用性强
gene_ranks <- setNames(-log10(res$pvalue) * sign(res$log2FoldChange), res$gene)

# 3. log2FC — 最简单，但忽略了显著性信息
gene_ranks <- setNames(res$log2FoldChange, res$gene)

# 注意：必须去除NA和Inf值！
gene_ranks <- gene_ranks[!is.na(gene_ranks) & is.finite(gene_ranks)]

# Leading Edge是GSEA结果的精华——贡献最大的核心基因
# 提取某个通路的leading edge基因
pathway_of_interest <- "HALLMARK_TNFA_SIGNALING_VIA_NFKB"
le_genes <- gsea_result[gsea_result$pathway == pathway_of_interest, ]$leadingEdge[[1]]
print(le_genes)  # 打印核心基因

# 批量提取所有显著通路的leading edge
le_all <- gsea_result[gsea_result$padj < 0.05, c("pathway", "leadingEdge")]

# fgsea快速流程
准备排名 → msigdbr()获取基因集 → fgsea() → plotEnrichment()

# clusterProfiler GSEA快速流程
准备排名 → gseGO()/gseKEGG() → gseaplot2()/dotplot()

# 排名统计量推荐
首选: DESeq2 stat列
次选: -log10(p) × sign(logFC)
简单: log2FC

# MSigDB常用类别
H: Hallmark（50个，首选）
C2: KEGG+Reactome
C5: GO基因集

# 显著性判断
fgsea原版: padj < 0.25（宽松）或 < 0.05（严格）
发文章推荐: padj < 0.05 + |NES| > 1