# 安装必要的包（如果未安装）
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

required_packages <- c("clusterProfiler", "org.Hs.eg.db", "enrichplot", 
                       "ggplot2", "DOSE", "stringr", "dplyr", "tidyr")
for (pkg in required_packages) {
  if (!require(pkg, character.only = TRUE)) {
    BiocManager::install(pkg)
  }
  library(pkg, character.only = TRUE)
}

# 读取基因列表（从你的数据）
# 假设你的基因列表在一个名为"genes.txt"的文件中，每行一个基因
# 或者你可以直接使用向量

# 这里我创建了基因列表向量（从你提供的数据）
gene_list <- c(
  "TWF1", "KLHL20", "PPARGC1B", "MKRN3", "B3GNT5", "RORA", "FZD3", "BRWD3", "CELSR3",
  "XPO1", "PDE7A", "NT5E", "EED", "RFX6", "RFX7", "KLHL28", "TNRC6A", "SCN9A", "CTNNB1",
  "MTDH", "NEDD4", "CHD1", "LMBR1", "PLAGL2", "FRZB", "SETD5", "SOX9", "EML1", "REEP3",
  "RUNX1", "FOXG1", "STK39", "VIM", "KLF10", "NFIB", "SPEN", "SH2B3", "NECAP1", "STIM2",
  "ITGA6", "RUNX2", "OTUD6B", "SAMD8", "SOCS1", "GOLGA1", "AZIN1", "CNKSR2", "CFL2",
  "PFN2", "COL13A1", "SOCS3", "ANKRD17", "EEA1", "CLOCK", "TNRC6B", "SEC23A", "TMOD2",
  "RAB8A", "WNT3A", "PSMD7", "BNIP3L", "SH3RF1", "RAPGEF4", "RAP1B", "CAMK2D", "LRRC8C",
  "KMT2C", "PAWR", "PGM1", "LPP", "TAOK1", "SCN1A", "PRLR", "ERLIN1", "SACS", "OXR1",
  "OSBPL8", "VPS26B", "PICALM", "ZNF711", "PPARGC1A", "MCF2L", "RASA1", "CAND1", "GATM",
  "PER2", "SIX4", "LRRC8D", "ARHGEF6", "MYO5A", "CALU", "PRKAA2", "DNMT3A", "EDNRA",
  "TENM3", "CBX2", "STX16", "ACTR1A", "NRIP1", "RAB23", "PDE4D", "ABL1", "LOX", "EPG5",
  "VAT1", "NDEL1", "ZDHHC17", "RHOB", "DLGAP1", "MAP3K5", "OSTM1", "SLC25A36", "PPP3CA",
  "FLVCR1", "CHD7", "EDEM3", "TMEM87A", "SNX10", "ELAVL2", "HOXA1", "NAA25", "PPP1R12A",
  "GJA1", "BCL9", "IDH1", "KMT2A", "PIP4K2B", "HIPK2", "MAML1", "ACVR1", "BCOR", "DDIT4",
  "ARID1A", "SLC38A2", "GLI2", "GALNT2", "DOCK7", "ATXN1", "ATP2A2", "SNAI1", "CSNK1G1",
  "NF1", "DSTYK", "CADPS", "ELOVL5", "GOLGA4", "IRF4", "NR4A2", "JOSD1", "NLGN1", "EFNA3",
  "CHMP2B", "PGM3", "IFNAR2", "RNF220", "SLC9A8", "PGP", "ST8SIA4", "IRS1", "ITGB3",
  "LRP6", "ZBTB18", "REEP1", "EPB41", "LIFR", "DNAJC13", "PDCD10", "INPP4A", "AP4E1",
  "MIA3", "EDC3", "DPYSL2", "UBE3C", "GCLC", "PI4K2B", "SLC4A7", "NUS1", "CAMK4", "PDGFRB",
  "RAPGEF2", "RTN4R", "TBC1D2B", "STAG2", "MAPK8", "SP4", "MAP3K7", "TRIO", "FAM13A",
  "ATP2B1", "ATL2", "GRB10", "PRPF40A", "CRKL", "IGF2R", "STT3B", "SLC5A3", "ZEB2",
  "PIGA", "TMCC1", "EPC2", "SLC6A6", "FIGN", "KCTD5", "TRAF3", "SBF1", "CADM2", "RASGEF1B",
  "NR3C1", "NEDD4L", "NAPG", "B4GALT6", "KRAS", "PAX9", "TNPO3", "RPS6KA2", "USP45",
  "LEPR", "SIK3", "AGO1", "ASCC3", "SLC7A6", "HNRNPC", "ATP2B2", "GCNT2", "QKI", "ITSN1",
  "MTTP", "CLCF1", "UGT8", "GALNT3", "AP3S1", "FXR1", "PCDH10", "HNRNPA2B1", "USP2",
  "MBNL1", "SALL4", "CPEB4", "LRCH2", "KCTD20", "BTBD10", "SLC6A9", "PHF6", "PLXNA2",
  "FAM217B", "TET1", "WDR26", "PPFIA2", "AFF3", "HNRNPA3", "IRS2", "MBNL2", "SETD3",
  "SLC7A11", "PPIP5K2", "SOS1", "WDFY3", "SEC23IP", "KCTD3", "SLC35A5", "TRPS1", "FYCO1",
  "TNKS", "NFIA", "ARF4", "SNAPIN", "LCOR", "LARP4", "UBE2D3", "HERC2", "JARID2", "ALPK3",
  "ZNF746", "GLUD1", "ZFYVE26", "NID1", "PBRM1", "PLXNA1", "SLC38A1", "SATB1", "DGKZ",
  "USP22", "UNKL", "PGM2L1", "CEP170B", "SON", "ZNF148", "MAN1B1", "ZBTB7A", "SIRT1",
  "MTF2", "YES1", "NCOA3", "SOX4", "RAB10", "SLC29A3", "H6PD", "NOVA1", "RNMT", "ACTN1",
  "NSD1", "CCNY", "PITX1", "JAG2", "USP15", "BAZ1A", "AGO2", "ADAM10", "MINPP1", "MIB1",
  "TIMP2", "TXNDC5", "SNX27", "FUCA1", "CHKA", "FOXP4", "SUCLG2", "CHD9", "GNAQ", "TBL1X",
  "APC", "RAB11A", "TIA1", "EPB41L4B", "DOC2A", "FRMD4A", "MMP9", "Axin 2", "CTH", "MECP2",
  "SMARCD2", "Cyclin D1", "TUBGCP3", "Axin 1"
)

# 注意：基因符号需要标准化（去除空格等）
gene_list <- gsub(" ", "", gene_list)  # 去除空格
gene_list <- unique(gene_list)  # 去重

cat("基因列表数量:", length(gene_list), "\n")

# 将基因符号转换为ENTREZID
gene_entrez <- bitr(gene_list, 
                    fromType = "SYMBOL", 
                    toType = "ENTREZID", 
                    OrgDb = org.Hs.eg.db)

cat("成功转换的基因数量:", nrow(gene_entrez), "\n")

if (nrow(gene_entrez) == 0) {
  stop("没有基因被成功转换，请检查基因符号是否正确")
}

# 1. GO富集分析
go_bp <- enrichGO(gene = gene_entrez$ENTREZID,
                  OrgDb = org.Hs.eg.db,
                  keyType = "ENTREZID",
                  ont = "BP",  # BP:生物过程，MF:分子功能，CC:细胞组分
                  pvalueCutoff = 0.05,
                  qvalueCutoff = 0.2,
                  readable = TRUE)

go_mf <- enrichGO(gene = gene_entrez$ENTREZID,
                  OrgDb = org.Hs.eg.db,
                  keyType = "ENTREZID",
                  ont = "MF",
                  pvalueCutoff = 0.05,
                  qvalueCutoff = 0.2,
                  readable = TRUE)

go_cc <- enrichGO(gene = gene_entrez$ENTREZID,
                  OrgDb = org.Hs.eg.db,
                  keyType = "ENTREZID",
                  ont = "CC",
                  pvalueCutoff = 0.05,
                  qvalueCutoff = 0.2,
                  readable = TRUE)

# 2. KEGG富集分析
kegg <- enrichKEGG(gene = gene_entrez$ENTREZID,
                   organism = 'hsa',  # 人类
                   pvalueCutoff = 0.05,
                   qvalueCutoff = 0.2)

# 将KEGG结果转换为可读的基因符号
if (!is.null(kegg)) {
  kegg <- setReadable(kegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}

# 保存结果到文件
output_dir <- "enrichment_results"
if (!dir.exists(output_dir)) {
  dir.create(output_dir)
}

# 保存GO结果
if (nrow(go_bp) > 0) {
  write.csv(go_bp, file.path(output_dir, "GO_BP_enrichment.csv"), row.names = FALSE)
  cat("GO BP结果已保存到", file.path(output_dir, "GO_BP_enrichment.csv"), "\n")
}

if (nrow(go_mf) > 0) {
  write.csv(go_mf, file.path(output_dir, "GO_MF_enrichment.csv"), row.names = FALSE)
  cat("GO MF结果已保存到", file.path(output_dir, "GO_MF_enrichment.csv"), "\n")
}

if (nrow(go_cc) > 0) {
  write.csv(go_cc, file.path(output_dir, "GO_CC_enrichment.csv"), row.names = FALSE)
  cat("GO CC结果已保存到", file.path(output_dir, "GO_CC_enrichment.csv"), "\n")
}

# 保存KEGG结果
if (!is.null(kegg) && nrow(kegg) > 0) {
  write.csv(kegg, file.path(output_dir, "KEGG_enrichment.csv"), row.names = FALSE)
  cat("KEGG结果已保存到", file.path(output_dir, "KEGG_enrichment.csv"), "\n")
}

# 3. 可视化
# 创建可视化目录
vis_dir <- file.path(output_dir, "visualization")
if (!dir.exists(vis_dir)) {
  dir.create(vis_dir)
}

# GO富集条形图（BP为例）
if (nrow(go_bp) > 0) {
  pdf(file.path(vis_dir, "GO_BP_barplot.pdf"), width = 12, height = 8)
  print(barplot(go_bp, showCategory = 20, title = "GO Biological Process Enrichment"))
  dev.off()
  
  # 点图
  pdf(file.path(vis_dir, "GO_BP_dotplot.pdf"), width = 12, height = 8)
  print(dotplot(go_bp, showCategory = 20, title = "GO Biological Process Enrichment"))
  dev.off()
  
  # 网络图（显示基因与GO term的关系）
  if (nrow(go_bp) >= 5) {
    pdf(file.path(vis_dir, "GO_BP_cnetplot.pdf"), width = 14, height = 10)
    print(cnetplot(go_bp, categorySize = "pvalue", foldChange = NULL,
                   showCategory = 10, node_label = "all"))
    dev.off()
  }
}

# KEGG可视化
if (!is.null(kegg) && nrow(kegg) > 0) {
  pdf(file.path(vis_dir, "KEGG_barplot.pdf"), width = 12, height = 8)
  print(barplot(kegg, showCategory = 20, title = "KEGG Pathway Enrichment"))
  dev.off()
  
  pdf(file.path(vis_dir, "KEGG_dotplot.pdf"), width = 12, height = 8)
  print(dotplot(kegg, showCategory = 20, title = "KEGG Pathway Enrichment"))
  dev.off()
  
  # KEGG通路图（需要安装pathview）
  tryCatch({
    if (!require("pathview", quietly = TRUE)) {
      BiocManager::install("pathview")
    }
    library(pathview)
    
    # 为每个显著富集的KEGG通路生成通路图
    sig_kegg_ids <- kegg@result$ID[1:min(5, nrow(kegg))]
    for (pid in sig_kegg_ids) {
      tryCatch({
        pathview(gene.data = gene_entrez$ENTREZID,
                 pathway.id = pid,
                 species = "hsa",
                 limit = list(gene = 5, cpd = 1))
      }, error = function(e) {
        cat("无法生成通路图", pid, ":", e$message, "\n")
      })
    }
    # 移动生成的PNG文件到可视化目录
    png_files <- list.files(pattern = "hsa.*\\.png$")
    if (length(png_files) > 0) {
      file.copy(png_files, vis_dir)
      file.remove(png_files)
    }
  }, error = function(e) {
    cat("路径可视化失败:", e$message, "\n")
  })
}

# 4. 结果汇总
cat("\n========== 富集分析结果汇总 ==========\n")
cat("输入基因总数:", length(gene_list), "\n")
cat("成功转换的基因数:", nrow(gene_entrez), "\n")

cat("\nGO富集结果:\n")
cat("  Biological Process:", ifelse(nrow(go_bp) > 0, paste(nrow(go_bp), "个显著条目"), "无显著条目"), "\n")
cat("  Molecular Function:", ifelse(nrow(go_mf) > 0, paste(nrow(go_mf), "个显著条目"), "无显著条目"), "\n")
cat("  Cellular Component:", ifelse(nrow(go_cc) > 0, paste(nrow(go_cc), "个显著条目"), "无显著条目"), "\n")

cat("\nKEGG富集结果:\n")
cat("  显著通路:", ifelse(!is.null(kegg) && nrow(kegg) > 0, paste(nrow(kegg), "个"), "无"), "\n")

cat("\n结果文件保存位置:\n")
cat("  数据文件:", output_dir, "\n")
cat("  可视化文件:", vis_dir, "\n")

# 5. 显示部分结果
if (nrow(go_bp) > 0) {
  cat("\nTop 10 GO Biological Process:\n")
  print(head(go_bp@result[, c("Description", "pvalue", "qvalue", "Count")], 10))
}

if (!is.null(kegg) && nrow(kegg) > 0) {
  cat("\nTop 10 KEGG Pathways:\n")
  print(head(kegg@result[, c("Description", "pvalue", "qvalue", "Count")], 10))
}

# 6. 可选：生成HTML报告
tryCatch({
  if (!require("rmarkdown", quietly = TRUE)) {
    install.packages("rmarkdown")
  }
  
  # 创建简化的报告
  report_file <- file.path(output_dir, "enrichment_report.Rmd")
  
  report_content <- '---
title: "基因富集分析报告"
output: html_document
date: "`r Sys.Date()`"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ggplot2)