数据挖掘—KEGG/GO分析中的富集因子、P值等概念理解

my_KEGG <- function(target_genes, species = "human", pvalue_cutoff = 1, qvalue_cutoff = 1) {
  # 根据物种选择相应的参数
  if (species == "human") {
    organism <- "hsa"
    OrgDb <- org.Hs.eg.db
  } else if (species == "mouse") {
    organism <- "mmu"
    OrgDb <- org.Mm.eg.db
  } else if (species == "rat") {
    organism <- "rno"
    OrgDb <- org.Rn.eg.db
  } else {
    stop("Unsupported species. Please choose 'human', 'mouse', or 'rat'.")
  }
  
  # 将基因符号转换为ENTREZ ID，并去除缺失值
  target_entrez <- as.character(na.omit(AnnotationDbi::select(
    OrgDb,
    keys = target_genes,
    columns = "ENTREZID",
    keytype = "SYMBOL"
  )[, 2]))
  
  # 打印转换后的基因数量
  message("转换后的基因数量: ", length(target_entrez))
  message("pvalue_cutoff: ", pvalue_cutoff)
  message("qvalue_cutoff: ", qvalue_cutoff)
  
  # 进行KEGG富集分析
  kk <- enrichKEGG(
    gene          = target_entrez,
    organism      = organism,
    pvalueCutoff  = pvalue_cutoff,
    qvalueCutoff  = qvalue_cutoff
  )
  
  # 将结果设置为可读形式
  kk <- setReadable(kk, OrgDb = OrgDb, keyType = "ENTREZID")
  
  # 返回结果
  return(kk)
}


kk <- my_KEGG(Target_all)
kk_result <- kk@result

# 解析 GeneRatio 和 BgRatio 字段为数值
kk_result$GeneRatio_num <- sapply(kk_result$GeneRatio, function(x) {
    eval(parse(text = x))  # 比如 "15/150" -> 0.1
})
kk_result$BgRatio_num <- sapply(kk_result$BgRatio, function(x) {
    eval(parse(text = x))  # 比如 "100/2000" -> 0.05
  })
  
# 计算富集因子并添加为新列
kk_result$enrichmentFactor <- kk_result$GeneRatio_num / kk_result$BgRatio_num
kk_result$GeneRatio_num <- NULL
kk_result$BgRatio_num <- NULL

table(kk_result$enrichmentFactor>1)
 #TRUE 
 #160 
table(kk_result$p.adjust < 0.05)
 #FALSE  TRUE 
 # 57   103 

#可视化
library(ggplot2)
# 添加 -log10(p.adjust) 列
kk_result$logP <- -log10(kk_result$p.adjust)
  
# 绘制散点图
ggplot(kk_result, aes(x = enrichmentFactor,y = logP)) +
    geom_point(aes(size = Count, color = logP)) +
    scale_color_gradient(low = "blue", high = "red") +
    labs(
      title = "Enrichment Factor vs -log10(p.adjust)",
      x = "Enrichment Factor",
      y = "-log10(Adjusted p-value)",
      color = "-log10(p.adjust)",
      size = "Gene Count"
    ) +
    theme_minimal()