还在使用同源基因转换获取非人类物种基因集吗，那你就out啦！来看看这个方法

生信技能树

发布于 2025-06-26 09:06:43

1560

今天分享业内大佬顾祖光的一个包，这个软件包旨在提供一个关于基因集富集分析的全面介绍，超棒的一个包，分享给大家~

https://jokergoo.github.io/GSEAtopics/index.html

首先是安装，非常简单：

## 使用西湖大学的 Bioconductor镜像
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
options("repos"=c(CRAN="https://mirrors.westlake.edu.cn/CRAN/"))

library(devtools)
install_github("jokergoo/GSEAtopics")

1.基因集来源

基因集顾名思义就是一堆基因的集合，在R语言中有三种格式的基因集：

1.1 list对象

格式如下：

lt = list(
    geneset1 = c("gene1", "gene2", "gene3"),
    geneset2 = c("gene2", "gene4"),
    geneset3 = c("gene1", "gene3", "gene5", "gene6")
)
lt

1.2 数据框格式：固定的两列

格式如下：有一些软件要求第一列是通路名第二列是对应的基因，有一些软件正好相反，这里需要看情况

df = data.frame(
    geneset = c(rep("geneset1", 3), rep("geneset2", 2), rep("geneset3", 4)),
    gene = c("gene1", "gene2", "gene3", "gene2", "gene4", "gene1", "gene3", "gene5", "gene6")
)
df
# df = df[, 2:1]

以上两种格式转换

list to data frame：

data.frame(
    geneset = rep(names(lt), times = sapply(lt, length)),
    gene = unlist(lt)
)

data frame to list：

split(df$gene, df$geneset)

GSAEtopics 包中有很方便的函数可以进行转换

# GSEAtopics 包中的函数进行转换
list_to_data_frame(lt)
data_frame_to_list(df)

1.3 二进制矩阵

这个格式用的不多，基因与基因集之间的关系可以表示为一个二进制矩阵：

m = matrix(0, nrow = 3, ncol = 6)
rownames(m) = unique(df$geneset)
colnames(m) = unique(df$gene)

for(i in seq_len(nrow(df))) {
    m[df[i, 1], df[i, 2]] = 1
}
m

如果基因有权重信息，这种格式会非常有用：

for(i in seq_len(nrow(df))) {
    m[df[i, 1], df[i, 2]] = runif(1)
}
m

2.获取不同数据库的基因集

我们常见的手段是去数据库下载各种通路与基因的对应关系，这里来看看如何通过代码获取常见数据库的基因集。

2.1 GO数据库

使用 org.*.db 包获取 GO terms和 gene的对应关系，以人类的 org.*.db 包为例

#### 
library(org.Hs.eg.db)
org.Hs.eg.db
keytypes(org.Hs.eg.db)
## 使用select
# 获取所有的基因
all_genes = keys(org.Hs.eg.db, keytype = "ENTREZID")
tb = select(org.Hs.eg.db, keys = all_genes, keytype = "ENTREZID", columns = c("GOALL", "ONTOLOGYALL"))
head(tb)
tail(tb)

# 提取bp,并去重
tb = tb[tb$ONTOLOGYALL == "BP", 1:2]
tb = unique(tb)
head(tb)
tail(tb)

此外，Biocondutor core team 包括了18个 org.*.db 物种包

2.2 KEGG数据库 gene sets

获取基因和通路对应关系：

#### KEGG
# 获取基因和通路对应关系：
df1 = read.table(url("https://rest.kegg.jp/link/pathway/hsa"), sep = "\t")
head(df1)
# 去掉基因和通路的前缀
df1[, 1] = gsub("hsa:", "", df1[, 1])
df1[, 2] = gsub("path:",  "", df1[, 2])
head(df1)

# 获取通路名称和通路ID对应关系
df2 = read.table(url("https://rest.kegg.jp/list/pathway/hsa"), sep = "\t")
head(df2)
# 去掉通路后面的物种名
df2[, 2] = gsub(" - Homo.*$", "", df2[, 2])
head(df2)

2.3 MSigDB数据库 gene sets

MSigDB定义了一种简单的.gmt格式用于存储基因集。

#### MSigDB
library(GSEAtopics)
list_msigdb_versions()
list_msigdb_collections("2025.1.Hs")
list_msigdb_collections("2025.1.Mm")

lt = get_msigdb(version = "2025.1.Mm", collection = "mh.all", gene_id_type = c("symbols"))
lt[1:2]

2.4 Reactome数据库 gene sets

Reactome是另一个流行的通路数据库。它以层次化的方式组织通路，其中包含通路和子通路或通路组件。最新的通路数据可以直接在以下网址找到：https://reactome.org/download-data（映射文件 -> 标识符映射文件 -> NCBI到所有通路）。

reactomePATHID2EXTID 包含了Reactome通路ID与基因EntreZ ID之间的映射关系。

reactomePATHID2NAME 包含了通路名称。

Reactome包含多个生物体的通路，在Reactome ID中，第二部分包含了生物体的信息。

# BiocManager::install("reactome.db")
library(reactome.db)
reactome.db
tb = toTable(reactomePATHID2EXTID)
head(tb)

p2n = toTable(reactomePATHID2NAME)
head(p2n)

# 统计物种
sort(table( gsub("^R-(\\w+)-\\d+$", "\\1", p2n[, 1]) ))

## 从数据库下载
download.file("https://reactome.org/download/current/NCBI2Reactome_All_Levels.txt", destfile = "NCBI2Reactome_All_Levels.txt")
tb = read.table("NCBI2Reactome_All_Levels.txt", sep = "\t", comment.char = "")
head(tb)

# 选取人的
tb2 = tb[grepl("HSA", tb[, 2]), ]
gs = split(tb2[, 1], tb2[, 2])
gs = lapply(gs, function(x) unique(as.character(x)))
gs[[1]]
head(gs)

2.5 UniProt 数据库 gene sets

UniProt数据库提供了一个以关键词形式表示的受控词汇列表，用于基因或蛋白质（https://www.uniprot.org/keywords/）。这对于以简洁的方式总结基因功能非常有用。

################ UniProt
# BiocManager::install("UniProtKeywords")
library(UniProtKeywords)
UniProtKeywords
# 获取可能得到数据的物种列表
UniProtKeywords:::ORGANISM

# 人: 可以使用三种关键词物种ID 获取
gl = load_keyword_genesets(9606)
gl[sample(length(gl), 2)]

gl = load_keyword_genesets("human")
gl[sample(length(gl), 2)]

gl = load_keyword_genesets("Homo sapiens")
gl[sample(length(gl), 2)]

# 返回表格格式
tb = load_keyword_genesets(9606, as_table = TRUE)
head(tb)
tail(tb)

2.6 GO/KEGG其他物种的基因集

KEGG pathway 支持大量的物种，可以通过一下方式获取：

read.table(url("https://rest.kegg.jp/list/pathway/hsa"), sep = "\t")
read.table(url("https://rest.kegg.jp/list/pathway/ptr"), sep = "\t")
read.table(url("https://rest.kegg.jp/list/pathway/pps"), sep = "\t")
read.table(url("https://rest.kegg.jp/list/pathway/ggo"), sep = "\t")
read.table(url("https://rest.kegg.jp/list/pathway/pon"), sep = "\t")
...

## 或者 GSEAtopics 这个包里面的获取方式
load_kegg_genesets(organism)

org.*.db 包，修改一下物种

select(org.*.db, keys = keys(org.*.db), columns = c("GOALL", "ONTOLOGYALL"))
select(org.*.db, keys = "BP", keytype = "ONTOLOGYALL", columns = c("ENTREZID", "GOALL"))
as.list(org.*.egGO2ALLEGS)

AnnotationHub 包支持 ~2000个物种

library(AnnotationHub)
ah = AnnotationHub()
query(ah, c("cat", "OrgDb"))
query(ah, c("Felis catus", "OrgDb"))
# It is annoying that ID changes between different package versions
org_db = ah[["AH117537"]]  # using `[[` downloads the dataset
columns(org_db)

## GO database
all_genes = keys(org_db, keytype = "ENTREZID")
tb = select(org_db, keys = all_genes, keytype = "ENTREZID", 
    columns = c("GOALL", "ONTOLOGYALL"))
head(tb)

tb = tb[!is.na(tb$GOALL), ]
tb = unique(tb)

# 或者
library(GSEAtopics)
lt = load_go_genesets(org_db, "BP")

今天分享到这里~

我们后面还会介绍罕见物种，非模式物种如何构建基因集~

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2025-06-25，如有侵权请联系 cloudcommunity@tencent.com 删除