社区首页 >专栏 >数据挖掘—KEGG/GO分析中的富集因子、P值等概念理解

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

原创

sheldor没耳朵

发布于 2025-05-07 03:47:42

7990

文章被收录于专栏：数据挖掘数据挖掘

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

基于这样的问题探究

1.概念理解

clusterProfiler 的结果表中没有直接显示富集因子enrichmentFactor，但是可以自己计算。首先还是辨析下，bgRatio（背景比例）、geneRatio（基因比例）和富集因子三个不同的指标

bgRatio（背景比例）：
- 这个指标表示在背景基因集中（即整个基因组或研究中考虑的所有基因），属于特定GO条目的基因所占的比例。
- 计算公式为：bgRatio=背景基因集的总基因数背景基因集中属于该GO条目的基因数
- 它反映了在没有特定条件下（如差异表达），某个GO条目在基因集中的普遍性。
geneRatio（基因比例）：
- 这个指标表示在差异表达基因集中，属于特定GO条目的基因所占的比例。
- 计算公式为：geneRatio=差异表达基因集的总基因数差异表达基因集中属于该GO条目的基因数
- 它反映了在特定条件下（如差异表达），某个GO条目在基因集中的富集程度。
富集因子（Enrichment Factor）：
- 富集因子是基因比例与背景比例的比值，用于衡量特定GO条目在差异表达基因集中相对于背景基因集的富集程度。
- 计算公式为：Enrichment Factor=bgRatio/geneRatio
- 如果富集因子大于1，表示该GO条目在差异表达基因集中富集；如果小于1，表示该GO条目在差异表达基因集中不富集。

2.问：为什么可以直接根据P值选top terms？

以KEGG为例，额外计算下enrichmentFactor，并作图展示enrichmentFactor与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

计算enrichmentFactor，可见enrichmentFactor全是>1的

# 解析 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()