如何批量查询单个基因的功能通路？

生信技能树

发布于 2025-02-19 13:33:26

9300

代码可运行

运行总次数：0

代码可运行

我们的马拉松授课最后一周给大家讲解了KEGG和GO数据库的网页版本的单个基因与多个基因的功能注释，还讲解了功能注释与功能富集的区别。富集其实就是显著注释到的功能通路。然后上课有个学员跟着实际操作演练，发现他的基因在KEGG数据库网页中没有找到通路：

KEGG 网页对某个基因的注释操作为：

学员的查询结果：

然后我验证了一下他的操作，真的没有找到！下面来看看其他的方法！

批量代码版本注释一个基因

学员的基因我们就不放了，我们随便找一个基因比如 EPCAM，一个上皮细胞的经典基因。我们上面讲了，基因的注释与富集含义是一样的，只不过富集是显著注释的通路，通过卡一个pvalue或者其他的阈值挑选得到，那我们第一种方法就用功能富集的方式拿到这个基因注释到的通路：

这里需要注意参数的设置，我们的目的是拿到所有可能的通路，所以设置了：pvalueCutoff= 1,qvalueCutoff= 1, minGSSize = 1, maxGSSize = 100000，ont为ALL，包括了CC，MF，BP三个集合的：

rm(list = ls())
library(clusterProfiler)
library(org.Hs.eg.db)
library(GSEABase)
library(ggplot2)
library(tidyverse)

# 单个基因查询
# GO通路
ego_ALL <- enrichGO(gene="EPCAM", OrgDb= 'org.Hs.eg.db', keyType='SYMBOL', ont="ALL", 
                   pvalueCutoff= 1,qvalueCutoff= 1, minGSSize = 1, maxGSSize = 100000)


# KEGG通路
# 拿到这个基因的ENTREZID
bitr(gene="EPCAM", fromType="SYMBOL", toType="ENTREZID", OrgDb='org.Hs.eg.db')
ekegg_ALL <- enrichKEGG(gene = "4072", organism = 'hsa', pvalueCutoff = 1, qvalueCutoff = 1,minGSSize = 1, maxGSSize = 100000)

以 GO数据结果为例，我们拿到了 138条通路：

具体的result：

还没有其他的方法？更全的通路：MsigDB

我们也可以拿到 MsigDB 数据库中所有含有EPCAM基因的通路！

首先去 MsigDB 数据库下载所有通路的gmt文件：msigdb.v2024.1.Hs.symbols.gmt，只有20几M，很快就下载完了！

https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp

## 所有通路
geneset <- read.gmt("msigdb.v2024.1.Hs.symbols.gmt")
geneset_EPCAM <- geneset[geneset$gene == "EPCAM",]
head(geneset_EPCAM)
dim(geneset_EPCAM)
str(geneset_EPCAM)