ORA

最后发布时间 : 2025-05-22 15:12:46 浏览量 :

Biomedical Knowledge Mining using GOSemSim and clusterProfiler

超几何检验(Hypergeometric Test)是一种常用于基因富集分析(如 GO/KEGG)的统计方法。下面是它的核心统计理念:

🎯 零假设(Null Hypothesis, H_0

零假设 H_0
你观察到的基因富集(例如某通路中的命中基因数)是完全随机的,不存在富集现象
换句话说,你的目标基因集在注释类(如某 GO term)中的出现频率,与整个背景基因集中的出现频率没有显著差异

🧮 统计量(Test Statistic)

超几何检验的统计量是一个离散随机变量 X,它表示:

从总共 N 个基因中,随机抽取 n 个基因,在一个已知包含 K 个目标基因的功能集合中,命中了 k 个目标基因的概率。

X \sim \text{Hypergeometric}(N, K, n)

  • N:总背景基因数
  • K:该功能类(如某 GO term)包含的基因数
  • n:你输入的感兴趣基因集合大小
  • k:你观察到输入基因集中命中该功能类的基因数
  • X统计量,你观测到的 k

📐 计算 p 值(右尾)

你实际计算的是:

P(X \geq k) = \sum_{i = k}^{\min(n, K)} \frac{\binom{K}{i} \binom{N-K}{n-i}}{\binom{N}{n}}

代表:在零假设成立的前提下,观察到至少 k 个命中的概率。

✅ 举例:

假设你有:

  • 背景基因集 N = 20,000
  • 某 GO term 中有 K = 500 个基因
  • 你的 DEGs 有 n = 300
  • 实际命中 GO term 的是 k = 20

就可以用超几何检验求 P(X \geq 20)

📊 R语言中超几何检验函数

# 计算 P(X >= k)
phyper(q = k - 1, m = K, n = N - K, k = n, lower.tail = FALSE)

如果你需要 Python 实现或对这部分代码更感兴趣,我也可以提供对应代码。是否需要?

KEGG

keggapi

devtools::install_github("BioinfoFungi/enrichplot2")

library(clusterProfiler)
library(org.Mm.eg.db)
library(KEGG.db)

annotations_orgDb <- AnnotationDbi::select(org.Mm.eg.db, 
keys = rownames(wt_control_deg),  
columns = c("ENTREZID","GENENAME"), 
keytype = "SYMBOL")
kegg<- clusterProfiler::enrichKEGG(gene =annotations_orgDb$ENTREZID,  organism = "mmu",
keyType = "SYMBOL",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
minGSSize = 10,
maxGSSize = 500,
qvalueCutoff = 10000,
use_internal_data = TRUE)
kegg <- setReadable(kegg, 'org.Mm.eg.db', 'ENTREZID')
png(filename = paste0(wt_control_vs_ko_control,"/KEGG.png"),width = 7,height = 6,res = 300,units = "in")
dotplot(kegg,showCategory=20,color="pvalue")
dev.off()
png(filename = paste0(wt_control_vs_ko_control,"/heat.png"),width = 10,height =5,res = 300,units = "in")
heatplot(kegg,showCategory=20,color="pvalue")
dev.off()
write_tsv(arrange(kegg@result,pvalue), file=paste0(wt_control_vs_ko_control,"/KEGG.tsv"))

注意pvalueCutoff=10000dotplot默认使用在enrichKEGG函数中设置的pvalueCutoff = 0.05筛选显著富集的通路。及使用0.05的阈值筛选矫正后的pvalue,进行可视化的。因此如果基因太少则不能富集到通路,此时可以将改值调大,或者修改dotplot函数的源码实现根据pvalue进行矫正。有关dotplot源码相关内容见enrichplot源码阅读

查看感兴趣基因与通路基因的交集

suppressMessages(library(KEGG.db))
suppressMessages(library(clusterProfiler))
suppressMessages(library(org.Mm.eg.db))
ls("package:KEGG.db")
frame = toTable(KEGGPATHID2EXTID)
x <- as.list(KEGG.db::KEGGPATHID2EXTID)
x["mmu04020"][[1]] |> length()
annotations <- AnnotationDbi::select(org.Mm.eg.db, 
keys = x["mmu04020"][[1]] ,  
columns = c("ENTREZID","SYMBOL"), 
keytype = "ENTREZID") |> 
na.omit()
annotations$SYMBOL |> length()
intersect(rownames(deg),annotations_foldchanges$SYMBOL)
ggvenn(
list("Calcium-signaling-pathway" = annotations_foldchanges$SYMBOL ,"mcherry-PKH vs mcherry deg" = rownames(deg) ),
stroke_size = 0.6, set_name_size = 2.8, text_size = 5,show_percentage=F
)

GO KEGG

library(clusterProfiler)
library(org.Hs.eg.db)
entrez_id <- bitr(top20_gene, fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db)%>%
pull("ENTREZID")%>%unique()
kegg<- clusterProfiler::enrichKEGG(gene =entrez_id,organism = 'hsa',pAdjustMethod="none")
kegg <- setReadable(kegg, 'org.Hs.eg.db', 'ENTREZID')
png(filename = "figure/KEGG.png",width = 7,height = 6,res = 300,units = "in")
dotplot(kegg,showCategory=20)
dev.off()
ego<-enrichGO(gene = entrez_id,OrgDb = org.Hs.eg.db,keyType = 'ENTREZID',ont = 'all',readable = T ,pAdjustMethod="BH")
png(filename = "figure/GO.png",width = 10,height = 9,res = 300,units = "in")
colorSel="qvalue"
dotplot(ego,showCategory = 10,split="ONTOLOGY",orderBy = "GeneRatio", color = colorSel) + facet_grid(ONTOLOGY~., scale='free')
dev.off()

KEGG本地化

https://github.com/YuLab-SMU/createKEGGdb

remotes::install_github("YuLab-SMU/createKEGGdb")
library(createKEGGdb)
species <-c("hsa","rno","mmu")
create_kegg_db(species)
install.packages("./KEGG.db_1.0.tar.gz", repos=NULL)
library(KEGG.db)
library(clusterProfiler)

data(geneList, package="DOSE")
gene <- names(geneList)[abs(geneList) > 2]

kk <- enrichKEGG(gene         = gene,
organism     = 'hsa',
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
use_internal_data =T)

合并GO结果

ALL

ONTOLOGYIDDescriptionGeneRatioBgRatiopvaluep.adjustqvaluegeneIDCount
BPGO:0006641triglyceride metabolic process2/5112/181123.744105884395701e-43.744105884395701e-40.009409082606148922Scd/Pnpla22
BPGO:0006639acylglycerol metabolic process2/5141/181125.92595747195959e-45.92595747195959e-40.009409082606148922Scd/Pnpla22
BPGO:0006638neutral lipid metabolic process2/5143/181126.09451941534646e-46.09451941534646e-40.009409082606148922Scd/Pnpla22
BPGO:0055088lipid homeostasis2/5171/181128.697839126288131e-48.697839126288131e-40.010071182146228362Scd/Pnpla22
BPGO:0010898positive regulation of triglyceride catabolic process1/512/181120.00330869900355446060.00330869900355446060.016408231859088875Pnpla21
BPGO:0048733sebaceous gland development1/513/181120.00358402809222746480.00358402809222746480.016408231859088875Scd1

BP

IDDescriptionGeneRatioBgRatiopvaluep.adjustqvaluegeneIDCount
GO:0006641triglyceride metabolic process2/5112/181123.744105884395701e-43.744105884395701e-40.009409082606148922Scd/Pnpla22
GO:0006639acylglycerol metabolic process2/5141/181125.92595747195959e-45.92595747195959e-40.009409082606148922Scd/Pnpla22
GO:0006638neutral lipid metabolic process2/5143/181126.09451941534646e-46.09451941534646e-40.009409082606148922Scd/Pnpla22
GO:0055088lipid homeostasis2/5171/181128.697839126288131e-48.697839126288131e-40.010071182146228362Scd/Pnpla22
GO:0010898positive regulation of triglyceride catabolic process1/512/181120.00330869900355446060.00330869900355446060.016408231859088875Pnpla21
GO:0048733sebaceous gland development1/513/181120.00358402809222746480.00358402809222746480.016408231859088875Scd1

Reference