# ============================================================
# 教程配套脚本 06：GO / KEGG 功能富集分析
# 输入：差异基因列表（gene symbol）
# 输出：富集结果表 + 气泡图/柱状图
# ============================================================

library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ggplot2)

deg <- read.csv("results/DEGs_tumor_vs_normal.csv")
gene_list <- deg$gene

# ---- 1. GO 富集（BP/CC/MF 全部）----
ego <- enrichGO(gene = gene_list,
                OrgDb = org.Hs.eg.db,
                keyType = "SYMBOL",          # 输入是基因 Symbol
                ont = "ALL",                 # BP+CC+MF
                pAdjustMethod = "BH",
                pvalueCutoff = 0.05,
                qvalueCutoff = 0.05,
                readable = TRUE)
head(as.data.frame(ego))

# ---- 2. KEGG 富集（人类 hsa）----
ekegg <- enrichKEGG(gene = gene_list,
                    organism = "hsa",
                    keyType = "kegg",
                    pvalueCutoff = 0.05)
# 若基因是 Symbol，先用 bitr 转 ENTREZID 再跑 KEGG 更稳妥：
# gene_entrez <- bitr(gene_list, fromType = "SYMBOL", toType = "ENTREZID",
#                     OrgDb = org.Hs.eg.db)
# ekegg <- enrichKEGG(gene = gene_entrez$ENTREZID, organism = "hsa")

# ---- 3. 保存结果 ----
write.csv(as.data.frame(ego), "results/GO_enrichment.csv", row.names = FALSE)
if (!is.null(ekegg)) write.csv(as.data.frame(ekegg), "results/KEGG_enrichment.csv", row.names = FALSE)

# ---- 4. 可视化：气泡图 ----
if (nrow(as.data.frame(ego)) > 0) {
  p_go <- dotplot(ego, showCategory = 15) +
    ggtitle("GO enrichment") +
    theme(axis.text.y = element_text(size = 8))
  ggsave("figures/R_go_dotplot.png", p_go, width = 9, height = 7, dpi = 300)
}

if (!is.null(ekegg) && nrow(as.data.frame(ekegg)) > 0) {
  p_kegg <- dotplot(ekegg, showCategory = 15) +
    ggtitle("KEGG enrichment") +
    theme(axis.text.y = element_text(size = 8))
  ggsave("figures/R_kegg_dotplot.png", p_kegg, width = 9, height = 6, dpi = 300)
}

# ---- 5. （进阶）GSEA：用全部基因的 log2FC 排序 ----
# gene_rank <- results %>% arrange(desc(logFC)) %>% select(gene, logFC) %>%
#   deframe()   # 命名向量：基因 -> log2FC
# gsea_res <- gseKEGG(geneList = gene_rank, organism = "hsa",
#                     pvalueCutoff = 0.25, seed = 123)
# gseaplot2(gsea_res, geneSetID = 1)   # 画富集曲线

cat("富集分析完成！\n")
