clusterProfiler — 基因功能富集分析与可视化 R 包¶

# 安装 clusterProfiler 及常用依赖
if (!require("BiocManager")) install.packages("BiocManager")

BiocManager::install(c(
    "clusterProfiler",    # 主包（当前 v4.12+）
    "org.Hs.eg.db",       # 人类基因注释数据库
    "org.Mm.eg.db",       # 小鼠基因注释数据库
    "DOSE",               # 疾病本体富集分析
    "enrichplot"          # 富集分析可视化
))
install.packages("ggplot2")

# 验证安装
library(clusterProfiler)
packageVersion("clusterProfiler")    # 应显示 4.12 或更高

library(clusterProfiler)
library(org.Hs.eg.db)   # 人类基因注释
library(enrichplot)

# 准备输入：差异基因列表（基因 Symbol）
deg_genes <- c("GLUT1","HK2","LDHA","PKM","ENO1",  # 上调基因示例
               "FOXO1","PPARA","ACOX1","CPT1A")

# ORA GO 富集（三个本体：BP/CC/MF）
ego <- enrichGO(
    gene = deg_genes,                  # 输入基因列表
    OrgDb = org.Hs.eg.db,             # 物种注释数据库
    keyType = "SYMBOL",                # 基因 ID 类型（SYMBOL/ENTREZID/ENSEMBL）
    ont = "BP",                        # 本体类型（BP/CC/MF/ALL）
    pAdjustMethod = "BH",              # 多重校正方法（BH=Benjamini-Hochberg）
    pvalueCutoff = 0.05,               # p 值阈值
    qvalueCutoff = 0.2,                # q 值阈值
    readable = TRUE                    # 将基因 ID 转为可读 Symbol
)

# 查看结果
head(ego@result, 10)                   # 前 10 个 GO term
dim(ego@result)                        # 富集到的 GO term 总数

# KEGG 分析需要 Entrez ID（不是 Symbol）
# 基因 ID 转换（Symbol → Entrez ID）
gene_entrez <- bitr(deg_genes,
                    fromType = "SYMBOL",    # 输入类型
                    toType = "ENTREZID",    # 目标类型
                    OrgDb = org.Hs.eg.db)  # 使用人类注释数据库

# KEGG 通路富集
kk <- enrichKEGG(
    gene = gene_entrez$ENTREZID,       # Entrez ID
    organism = "hsa",                  # 物种代码（hsa=人类）
    pvalueCutoff = 0.05,
    pAdjustMethod = "BH"
)
head(kk@result, 10)

# 准备有序基因列表（Entrez ID + 排序统计量）
deseq_res <- read.csv("DESeq2_results.csv", row.names=1)

# 转换基因 ID
gene_list_df <- bitr(rownames(deseq_res),
                     fromType="SYMBOL",
                     toType="ENTREZID",
                     OrgDb=org.Hs.eg.db)

# 合并统计量
deseq_merged <- merge(gene_list_df, deseq_res,
                       by.x="SYMBOL", by.y=0)

# 创建有序数值向量
geneList <- setNames(deseq_merged$stat,          # DESeq2 Wald 统计量
                      deseq_merged$ENTREZID)      # Entrez ID 为名称
geneList <- sort(geneList, decreasing=TRUE)       # 降序排列

# GSEA + KEGG
gse_kegg <- gseKEGG(
    geneList = geneList,
    organism = "hsa",
    minGSSize = 10,                               # 最小基因集大小
    maxGSSize = 500,
    pvalueCutoff = 0.05,
    verbose = FALSE
)

# GSEA + GO
gse_go <- gseGO(
    geneList = geneList,
    OrgDb = org.Hs.eg.db,
    ont = "BP",
    minGSSize = 10,
    maxGSSize = 500,
    pvalueCutoff = 0.05,
    verbose = FALSE
)

library(enrichplot)

# 1. 条形图（barplot）— 最直观
barplot(ego, showCategory=20)          # 展示前 20 个 GO term

# 2. 气泡图（dotplot）— 展示更多信息
dotplot(ego, showCategory=20)          # 气泡大小=基因数，颜色=p.adjust

# 3. 富集图（emapplot）— 通路关联网络
ego_pair <- pairwise_termsim(ego)      # 计算 GO term 相似性
emapplot(ego_pair, showCategory=30)    # 网络图（相似 term 聚在一起）

# 4. 基因概念网络图（cnetplot）— 基因和通路的关联
cnetplot(ego,
         categorySize="pvalue",        # 节点大小=p值
         foldChange=setNames(deseq_res$log2FoldChange, rownames(deseq_res)))

# 5. GSEA 富集图（ridgeplot）
ridgeplot(gse_kegg) +
    labs(x="Enrichment Distribution")

# 6. GSEA 走势图（gseaplot2）
gseaplot2(gse_kegg,
          geneSetID=1:3,               # 展示前 3 个通路
          title="KEGG GSEA Results")

# GO ORA 分析（最简）
library(clusterProfiler); library(org.Hs.eg.db)
ego <- enrichGO(gene=deg_genes, OrgDb=org.Hs.eg.db, keyType="SYMBOL", ont="BP", pvalueCutoff=0.05)
dotplot(ego, showCategory=20)

# KEGG ORA 分析
entrez <- bitr(deg_genes, fromType="SYMBOL", toType="ENTREZID", OrgDb=org.Hs.eg.db)
kk <- enrichKEGG(gene=entrez$ENTREZID, organism="hsa", pvalueCutoff=0.05)

# GSEA-GO
gse <- gseGO(geneList=geneList, OrgDb=org.Hs.eg.db, ont="BP", minGSSize=10, pvalueCutoff=0.05)