KEGG.db的包自2011年后不再更新，clusterprofiler做KEGG分析的時候，可以使用KEGG提供的API對所需物種信息進(jìn)行獲取，但是KEGG的服務(wù)器有時候會斷網(wǎng)，這對我們批量進(jìn)行數(shù)據(jù)分析的時候，會發(fā)生錯誤而分析中斷。另外，目前我們在云上使用docker是無法聯(lián)網(wǎng)的（不是云平臺的原因）。因此我需要自己自備數(shù)據(jù)分析來獲得更為全面正確的分析結(jié)果。
> library('KEGG.db')
> ls('package:KEGG.db') #
[1] 'KEGG' 'KEGG_dbconn' 'KEGG_dbfile' 'KEGG_dbInfo' 'KEGG_dbschema'
[6] 'KEGGENZYMEID2GO' 'KEGGEXTID2PATHID' 'KEGGGO2ENZYMEID' 'KEGGMAPCOUNTS' 'KEGGPATHID2EXTID'
[11] 'KEGGPATHID2NAME' 'KEGGPATHNAME2ID'
> columns(KEGG.db)
Error in columns(KEGG.db) : object 'KEGG.db' not found
> # columns(org.Hs.eg.db) # show colums in .db
> # columns(KEGG.db) # Error in columns(KEGG.db) : object 'KEGG.db' not found
> # available.dbschemas()
> frame = toTable(KEGGPATHID2EXTID)
> head(frame)
pathway_id gene_or_orf_id
1 hsa00232 10
2 hsa00983 10
3 hsa01100 10
4 hsa00230 100
5 hsa01100 100
6 hsa05340 100
> frame = toTable(KEGGPATHID2NAME)
> head(frame)
path_id path_name
1 hsa00010 Glycolysis / Gluconeogenesis
2 hsa00020 Citrate cycle (TCA cycle)
3 hsa00030 Pentose phosphate pathway
4 hsa00040 Pentose and glucuronate interconversions
5 hsa00051 Fructose and mannose metabolism
6 hsa00052 Galactose metabolism
> class(KEGGPATHID2EXTID)
[1] 'AnnDbBimap'
attr(,'package')
[1] 'AnnotationDbi'
ls('package:KEGG.db')返回的結(jié)果是一些AnnObObj (Bimap objects)。是一種老的AnnotationDbi的interface，現(xiàn)在已經(jīng)不怎么建議使用了?，F(xiàn)在推薦使用select方法（columns()、cols()、keytypes()等）。
KEGG（目錄，手動創(chuàng)建）
├── DESCRIPTION（文件，手動創(chuàng)建）
├── inst（目錄）
│ └── extdata（目錄）
│ └── KEGG.sqlite（文件，通過代碼生成）
├── LICENSE（文件，手動創(chuàng)建）
├── NAMESPACE（文件，手動創(chuàng)建）
└── R（目錄，手動創(chuàng)建）
└── zzz.R（文件，手動創(chuàng)建）
參考2011版本的KEGG.db來寫：
Package: KEGG.db
Title: KEGG.db for KEGG enrichment analysis.
Description: KEGG.db for KEGG enrichment analysis.
Version: 1.0
Author: xxx <xxx@xxx.xxx.com>
Maintainer: xxx <xxx@xxx.xxx.com>
Depends: R (>= 2.7.0), methods, AnnotationDbi (>= 1.44.0)
Imports: methods, AnnotationDbi
Suggests: DBI
License: BSD
License_restricts_use: yes
biocViews: AnnotationData, FunctionalAnnotation
參考2011版本的KEGG.db來寫：
Free for academic use. Non-academic users are requested to obtain a license agreement with KEGG.
參考2011版本的KEGG.db來寫：
import(methods)
import(AnnotationDbi)
### Only put what is statically exported here. All the AnnObj instances
### created at load time are dynamically exported (refer to R/zzz.R for
### the details).
export(
KEGG,
KEGG_dbconn,
KEGG_dbfile,
KEGG_dbschema,
KEGG_dbInfo
)
查看clusterprofiler的enrichKEGG函數(shù)，當(dāng)使用KEGG.db的數(shù)據(jù)時，使用了get_data_from_KEGG_db()函數(shù)進(jìn)行數(shù)據(jù)獲取。
get_data_from_KEGG_db <- function(species) {
PATHID2EXTID <- as.list(get_KEGG_db('KEGGPATHID2EXTID'))
if (!any(grepl(species, names(PATHID2EXTID)))) {
stop('input species is not supported by KEGG.db...')
}
idx <- grep(species, names(PATHID2EXTID))
PATHID2EXTID <- PATHID2EXTID[idx]
PATHID2EXTID.df <- stack(PATHID2EXTID)
PATHID2EXTID.df <- PATHID2EXTID.df[, c(2,1)]
PATHID2NAME <- as.list(get_KEGG_db('KEGGPATHID2NAME'))
PATHID2NAME.df <- data.frame(path=paste0(species, names(PATHID2NAME)),
name=unlist(PATHID2NAME))
build_Anno(PATHID2EXTID.df, PATHID2NAME.df)
}
get_KEGG_db <- function(kw) {
annoDb <- 'KEGG.db'
suppressMessages(requireNamespace(annoDb))
eval(parse(text=paste0(annoDb, '::', kw)))
}
get_data_from_KEGG_db()使用get_KEGG_db()函數(shù)從'KEGG.db'中獲取了'KEGGPATHID2EXTID'和'KEGGPATHID2NAME'這兩個AnnDbBimap。
具體來說，對于get_KEGG_db('KEGGPATHID2EXTID')其實就是執(zhí)行KEGG.db::KEGGPATHID2EXTID這一條命令。as.list(get_KEGG_db('KEGGPATHID2EXTID'))就是將KEGG.db::KEGGPATHID2EXTID這一條命令返回的結(jié)果轉(zhuǎn)成list。
> x <- as.list(KEGG.db::KEGGPATHID2EXTID)
> head(x,2)
$hsa00232
[1] '10' '1544' '1548' '1549' '1553' '7498' '9'
$hsa00983
[1] '10' '1066' '10720' '10941' '151531' '1548' '1549' '1551' '1553'
[10] '1576' '1577' '1806' '1807' '1890' '221223' '2990' '3251' '3614'
[19] '3615' '3704' '51733' '54490' '54575' '54576' '54577' '54578' '54579'
[28] '54600' '54657' '54658' '54659' '54963' '574537' '64816' '7083' '7084'
[37] '7172' '7363' '7364' '7365' '7366' '7367' '7371' '7372' '7378'
[46] '7498' '79799' '83549' '8824' '8833' '9' '978'
那AnnDbBimap是怎么生成的？在AnnotationDbi的createAnnObjs.KEGG_DB.R文件里有生成bimap的代碼：
KEGG_DB_AnnDbBimap_seeds <- list(
list(
objName='PATHID2NAME',
Class='AnnDbBimap',
L2Rchain=list(
list(
tablename='pathway2name',
Lcolname='path_id',
Rcolname='path_name'
)
)
),
list(
objName='PATHID2EXTID',
Class='AnnDbBimap',
L2Rchain=list(
list(
tablename='pathway2gene',
Lcolname='pathway_id',
Rcolname='gene_or_orf_id'
)
)
),
list(
objName='ENZYMEID2GO',
Class='AnnDbBimap',
L2Rchain=list(
list(
tablename='ec2go',
Lcolname='ec_no',
Rcolname='go_id'
)
)
)
)
createAnnObjs.KEGG_DB <- function(prefix, objTarget, dbconn, datacache)
{
checkDBSCHEMA(dbconn, 'KEGG_DB')
## AnnDbBimap objects
seed0 <- list(
objTarget=objTarget,
datacache=datacache
)
ann_objs <- createAnnDbBimaps(KEGG_DB_AnnDbBimap_seeds, seed0)
## Reverse maps
ann_objs$PATHNAME2ID <- revmap(ann_objs$PATHID2NAME, objName='PATHNAME2ID')
ann_objs$EXTID2PATHID <- revmap(ann_objs$PATHID2EXTID, objName='EXTID2PATHID')
ann_objs$GO2ENZYMEID <- revmap(ann_objs$ENZYMEID2GO, objName='GO2ENZYMEID')
## 1 special map that is not an AnnDbBimap object (just a named integer vector)
ann_objs$MAPCOUNTS <- createMAPCOUNTS(dbconn, prefix)
prefixAnnObjNames(ann_objs, prefix)
}
具體參考2011版本的KEGG.db。
datacache <- new.env(hash=TRUE, parent=emptyenv())
KEGG <- function() showQCData('KEGG', datacache)
KEGG_dbconn <- function() dbconn(datacache)
KEGG_dbfile <- function() dbfile(datacache)
KEGG_dbschema <- function(file='', show.indices=FALSE) dbschema(datacache, file=file, show.indices=show.indices)
KEGG_dbInfo <- function() dbInfo(datacache)
.onLoad <- function(libname, pkgname)
{
## Connect to the SQLite DB
dbfile <- system.file('extdata', 'KEGG.sqlite', package=pkgname, lib.loc=libname)
assign('dbfile', dbfile, envir=datacache)
dbconn <- dbFileConnect(dbfile)
assign('dbconn', dbconn, envir=datacache)
## Create the AnnObj instances
ann_objs <- createAnnObjs.SchemaChoice('KEGG_DB', 'KEGG', 'KEGG', dbconn, datacache)
mergeToNamespaceAndExport(ann_objs, pkgname)
packageStartupMessage(AnnotationDbi:::annoStartupMessages('KEGG.db'))
}
.onUnload <- function(libpath)
{
dbFileDisconnect(KEGG_dbconn())
}
這里，為了簡化過程，我們通過AnnotationDbi包的createAnnObjs.SchemaChoice()函數(shù)調(diào)用了AnnotationDbi里的createAnnObjs.KEGG_DB.R來產(chǎn)生AnnDbBimap。你也可以自己自定義一個加到zzz.R里。具體操作可以參考AnnotationForge,version 2.11 of Bioconductor的Creating an annotation package with a new database schema
packagedir <- '你創(chuàng)建的KEGG目錄所在的路徑'
sqlite_path <- paste(packagedir, sep=.Platform$file.sep, 'inst', 'extdata')
if(!dir.exists(sqlite_path)){dir.create(sqlite_path,recursive = TRUE)}
dbfile <- file.path(sqlite_path, 'KEGG.sqlite')
unlink(dbfile)
###################################################
### download data
###################################################
species <- 'hsa'
# KEGG的物種列表查看：
# https://www.genome.jp/kegg/catalog/org_list.html
# 或者
# http://rest.kegg.jp/list/organism
kegg <- clusterProfiler::download_KEGG(species)
KEGGPATHID2NAME <- kegg$KEGGPATHID2NAME
colnames(KEGGPATHID2NAME) <- c('path_id','path_name')
KEGGPATHID2NAME$path_id <- sub(species,'',KEGGPATHID2NAME$path_id)
KEGGPATHID2EXTID <- kegg$KEGGPATHID2EXTID
colnames(KEGGPATHID2EXTID) <- c('pathway_id','gene_or_orf_id')
###################################################
### create database
###################################################
## Create the database file
library(RSQLite)
drv <- dbDriver('SQLite')
db <- dbConnect(drv, dbname=dbfile)
## dbDisconnect(db)
###################################################
### put the data into the tables
###################################################
dbWriteTable(conn = db, 'pathway2name', KEGGPATHID2NAME, row.names=FALSE)
dbWriteTable(conn = db, 'pathway2gene', KEGGPATHID2EXTID, row.names=FALSE)
dbListTables(db)
test <- dbReadTable(conn = db,'pathway2name')
head(test)
test <- dbReadTable(conn = db,'pathway2gene')
head(test)
###################################################
### append the metadata
###################################################
#dbSendQuery(conn = db,'drop table if exists metadata')
metadata <- rbind(c('PATHNAMESOURCENAME', 'KEGG PATHWAY'),
c('PATHNAMESOURCEURL', 'ftp://ftp.genome.jp/pub/kegg/pathway'),
c('PATHNAMESOURCEDATE', format(Sys.Date(), '%Y%m%d')),
c('KEGGSOURCENAME', 'KEGG GENOME'),
c('KEGGSOURCEURL', 'ftp://ftp.genome.jp/pub/kegg/genomes'),
c('KEGGSOURCEDATE', format(Sys.Date(), '%Y%m%d')),
c('GOEXTSOURCEDATE', '2015-Sepec2go27'),
c('GOEXTSOURCENAME', 'Gene Ontology External Link'),
c('GOEXTSOURCEURL', 'http://www.geneontology.org/external2go'),
c('Db type', 'KEGGDB'),
c('DBSCHEMA', 'KEGG_DB'),
c('DBSCHEMAVERSION', '2.1'))
metadata <- as.data.frame(metadata)
colnames(metadata) <- c('name', 'value') #makeAnnDbPkg規(guī)定的
dbWriteTable(conn = db, 'metadata', metadata, row.names=FALSE)
map.counts <- rbind(c('pathway2name', nrow(KEGGPATHID2NAME)),
c('pathway2gene', nrow(KEGGPATHID2EXTID)))
map.counts <- as.data.frame(map.counts)
colnames(map.counts) <- c('map_name','count')
dbWriteTable(conn = db, 'map_counts', map.counts, row.names=FALSE)
dbListTables(db)
dbListFields(conn = db, 'metadata')
dbReadTable(conn = db,'metadata')
map.metadata <- rbind(c('ENZYMEID2GO','Gene Ontology External Link','http://www.geneontology.org/external2go','2015-Sepec2go27'),
c('GO2ENZYMEID','Gene Ontology External Link','http://www.geneontology.org/external2go','2015-Sepec2go27'),
c('EXTID2PATHID','KEGG GENOME','ftp://ftp.genome.jp/pub/kegg/genomes','2011-Mar15'),
c('PATHID2EXTID','KEGG GENOME','ftp://ftp.genome.jp/pub/kegg/genomes','2011-Mar15'),
c('PATHNAME2ID','KEGG PATHWAY','ftp://ftp.genome.jp/pub/kegg/pathway',format(Sys.Date(),'%Y%m%d')),
c('PATHID2NAME','KEGG PATHWAY','ftp://ftp.genome.jp/pub/kegg/pathway',format(Sys.Date(),'%Y%m%d')))
map.metadata <- as.data.frame(map.metadata)
colnames(map.metadata) <- c('map_name','source_name','source_url','source_date')
dbWriteTable(conn = db, 'map_metadata', map.metadata, row.names=FALSE)
dbListTables(db)
dbListFields(conn = db, 'map_metadata')
dbReadTable(conn = db,'map_metadata')
dbDisconnect(db)
新開一個RStudio ：File → New Project
install.packages('D:/test/KEGG.db_1.0.tar.gz', repos=NULL,type = 'source')
library(KEGG.db)
library(clusterProfiler)
data(geneList, package='DOSE')
head(geneList)
gene <- names(geneList)[abs(geneList) > 2]
head(gene)
kk <- enrichKEGG(gene = gene,
organism = 'hsa',
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
use_internal_data =T)
# 使用線上拉取的數(shù)據(jù)，則use_internal_data = F
head(kk)
nrow(kk)
結(jié)果對比：
好了，打完收工。
注：最簡單粗暴的方式其實是在安裝了KEGG.db后，使用本文添加數(shù)據(jù)的方法，將sqlite里的數(shù)據(jù)直接替換掉。
上面的內(nèi)容是「生信媛」公眾號的小伙伴投的稿。 寫得詳細(xì)，雖然說follow下來的人都可以創(chuàng)建屬于自己的KEGG.db，但對于小白來說，感覺看到這么多代碼，會蒙圈。當(dāng)然打好包供大家下載也不是辦法，一來可能KEGG會找麻煩，他自己是開放給人家去爬數(shù)據(jù)，但你打包給別人用可能是不允許的。第二那么多的物種，提供給你的可能不是你想要。第三打包還得隔段時間更新。 在我看完這篇文章之后，我當(dāng)場就提議，不如一起寫個R包，這個R包的功能，就是你給個物種名，它負(fù)責(zé)去抓數(shù)據(jù)，然后生成KEGG.db，也就是一鍵到最后一步，用戶只要安裝，就可以在clusterProfiler中離線用最新的KEGG數(shù)據(jù)了。完美～ 寫個R包，讓它去生成另一個R包，這操作我也覺得挺6的開始Y叔叔的表演
remotes::install_github('YuLab-SMU/createKEGGdb')
只有一條命令，create_kegg_db，參數(shù)是你給定的物種，比如人的是hsa，函數(shù)的輸出就是在當(dāng)前目錄下生成一個KEGG.db_1.0.tar.gz的R包。
也是一條指令：
install.packages('KEGG.db_1.0.tar.gz', type='source')
然后你就可以使用前面文章中的第8步用clusterProfiler去分析了。
在線抓最新數(shù)據(jù)固然好，但不好重復(fù)啊，你分析了之后，過了幾個月了，準(zhǔn)備投稿，發(fā)現(xiàn)之前選定的通路重復(fù)不出來了，欲哭無淚怎么辦？現(xiàn)在自己打KEGG.db包，也就是你把時間點給固定了，什么時候重新跑，只要你用的還是原來的KEGG.db，你必須就可以重復(fù)出來。也就是它給了你可以回到過去的「月光寶盒」。
聯(lián)系客服
