如何计算基序富集p值？请求验证R代码正确性

问题背景 #

用户提供如下基因集合与基序统计数据：

Category,Total_Genes,UUAGGG_motif
Background,22591,18190
SetA,122,102
SetB,198,182
SetC,90,82

需要计算SetA、SetB、SetC分别与Background对比的p值，判断哪个集合中UUAGGG基序存在富集，并验证自己编写的R代码是否正确。

用户编写的R代码：

library(hypeR)

# Number of background genes
N <- 22591
# Number of background genes with motif
K <- 18190

# Set A
n_A <- 122
k_A <- 102

# Set B
n_B <- 198
k_B <- 182

# Set C
n_C <- 90
k_C <- 82

# Perform hypergeometric test for Set A
p_value_A <- 1 - phyper(k_A - 1, K, N - K, n_A, lower.tail = TRUE)

# Perform hypergeometric test for Set B
p_value_B <- 1 - phyper(k_B - 1, K, N - K, n_B, lower.tail = TRUE)

# Perform hypergeometric test for Set C
p_value_C <- 1 - phyper(k_C - 1, K, N - K, n_C, lower.tail = TRUE)

代码正确性分析 #

你的代码是正确的，理由如下：

• 超几何检验是基因集富集分析的标准方法，完全适配这种"从背景中抽取子集，检验子集内目标特征（带基序基因）比例是否显著高于背景"的场景。
• phyper函数的参数对应逻辑准确：
  • k_A - 1：要计算"至少观测到k_A个带基序基因"的概率，需用1减去累积到k_A-1的左侧概率
  • K：背景中带基序的基因总数（总体成功数）
  • N - K：背景中不带基序的基因总数（总体失败数）
  • n_A：当前集合的基因总数（抽取的样本量）
• 用1 - phyper(..., lower.tail=TRUE)计算右侧尾概率，正好对应"集合中带基序基因数显著多于随机预期"的富集检验逻辑。

额外建议 #

1. 批量处理优化：可以把数据整理成数据框，用循环或向量化操作批量计算，避免重复代码：

df <- read.csv(text = "Category,Total_Genes,UUAGGG_motif
Background,22591,18190
SetA,122,102
SetB,198,182
SetC,90,82")

# 提取背景数据
bg_total <- df$Total_Genes[df$Category == "Background"]
bg_motif <- df$UUAGGG_motif[df$Category == "Background"]

# 处理目标集合
test_sets <- df[df$Category != "Background", ]
test_sets$p_value <- sapply(1:nrow(test_sets), function(i) {
  n <- test_sets$Total_Genes[i]
  k <- test_sets$UUAGGG_motif[i]
  1 - phyper(k - 1, bg_motif, bg_total - bg_motif, n, lower.tail = TRUE)
})

2. 多重检验校正：由于同时对3个集合做检验，建议对p值进行多重校正（比如Bonferroni或FDR校正），避免假阳性结果：

test_sets$adj_p_value <- p.adjust(test_sets$p_value, method = "fdr")

3. 利用hypeR封装函数：你加载了hypeR包，其实可以直接用它的hyper_enrichment函数，封装更完善，还能输出可视化结果：

# 构造基因列表示例（假设你有具体基因ID）
setA_genes <- paste0("gene_", 1:122)
bg_genes <- paste0("gene_", 1:22591)
motif_genes <- paste0("gene_", 1:18190)

# 运行富集分析
hyp_obj <- hyper_enrichment(setA_genes, bg_genes, genesets = list(UUAGGG = motif_genes))
print(hyp_obj)

内容的提问来源于stack exchange，提问作者Apex