最近,有粉丝运行了我以前的数据挖掘成套代码里面的 run_kegg 函数,如下所示:
library(clusterProfiler)
run_kegg(gene_up,gene_down,pro='comp1')
出现了如下所示的报错:
Reading KEGG annotation online:
fail to download KEGG data...
Error in download.KEGG.Path(species) :
'species' should be one of organisms listed in 'http://www.genome.jp/kegg/catalog/org_list.html'...
In addition: Warning message:
In utils::download.file(url, quiet = TRUE, method = method, ...) :
URL 'https://rest.kegg.jp/link/mmu/pathway': status was 'Failure when receiving data from the peer'
Called from: download.KEGG.Path(species)
然后就找我,以为是我们的标准代码有问题,实际上我的 run_kegg 函数仅仅是包装了 Y叔的 clusterProfiler包而已 ,实际上里面没有啥玄机,如下所示:
## KEGG pathway analysis
### 做KEGG数据集超几何分布检验分析,重点在结果的可视化及生物学意义的理解。
run_kegg <- function(gene_up,gene_down,pro='test'){
library(ggplot2)
gene_up=unique(gene_up)
gene_down=unique(gene_down)
gene_diff=unique(c(gene_up,gene_down))
### over-representation test
# 下面把3个基因集分开做超几何分布检验
# 首先是上调基因集。
kk.up <- enrichKEGG(gene = gene_up,
organism = 'mmu',
#universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.up)[,1:6]
kk=kk.up
dotplot(kk)
kk=DOSE::setReadable(kk, OrgDb='org.Mm.eg.db',keyType='ENTREZID')
write.csv(kk@result,paste0(pro,'_kk.up.csv'))
# 首先是下调基因集。
kk.down <- enrichKEGG(gene = gene_down,
organism = 'mmu',
#universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.down)[,1:6]
kk=kk.down
dotplot(kk)
kk=DOSE::setReadable(kk, OrgDb='org.Mm.eg.db',keyType='ENTREZID')
write.csv(kk@result,paste0(pro,'_kk.down.csv'))
# 最后是上下调合并后的基因集。
kk.diff <- enrichKEGG(gene = gene_diff,
organism = 'mmu',
pvalueCutoff = 0.05)
head(kk.diff)[,1:6]
kk=kk.diff
dotplot(kk)
kk=DOSE::setReadable(kk, OrgDb='org.Mm.eg.db',keyType='ENTREZID')
write.csv(kk@result,paste0(pro,'_kk.diff.csv'))
}
就是针对上下调基因,以及其合并的差异基因,分别独立走 enrichKEGG 函数而已,它来自于 Y叔的 clusterProfiler包。
实际上,我们根据前面的报错,进入 https://www.genome.jp/kegg/catalog/org_list.html 可以很清楚的看到,我们的人和鼠,都是在里面的:
hsa Homo sapiens (human)
mmu Mus musculus (house mouse)
所以不可能是物种问题,而且这个kegg数据库的官网也可以访问,那么合理的推测,是不是Y叔的 clusterProfiler包有问题呢?其实仔细看报错提示是 download.KEGG.Path 函数问题,也是Y叔的包里面的。这个download.KEGG.Path 函数的代码是:
keggpathid2extid.df <- kegg_link(species, "pathway")
if (is.null(keggpathid2extid.df))
stop("'species' should be one of organisms listed in 'http://www.genome.jp/kegg/catalog/org_list.html'...")
keggpathid2extid.df[, 1] %<>% gsub("[^:]+:", "",
.)
keggpathid2extid.df[, 2] %<>% gsub("[^:]+:", "",
.)
keggpathid2name.df <- kegg_list("pathway")
keggpathid2name.df[, 1] %<>% gsub("path:map", species,
.)
keggpathid2extid.df <- keggpathid2extid.df[keggpathid2extid.df[,
1] %in% keggpathid2name.df[, 1], ]
return(list(KEGGPATHID2EXTID = keggpathid2extid.df, KEGGPATHID2NAME = keggpathid2name.df))
也就是说在kegg_link函数就开始报错了,继续看kegg_link的代码:
url <- paste0("http://rest.kegg.jp/link/", target_db,
"/", source_db, collapse = "")
kegg_rest(url)
实际上, 它就根据我们的物种,拼凑了一个 网页,比如:https://rest.kegg.jp/link/hsa/pathway,可以使用浏览器打开,能看到里面的通路的id和基因id的对应关系。所以真正的问题是 kegg_rest 这个函数,没有办法像浏览器那样直接获取网页里面的内容。我们继续看kegg_rest 这个函数的源代码:
dl <- mydownload(rest_url, destfile = f)
if (is.null(dl)) {
message("fail to download KEGG data...")
return(NULL)
}
真气人,就是一个套娃, 里面又是mydownload函数,那我们继续看mydownload函数的源代码 :
if (is.null(method))
method <- getOption("clusterProfiler.download.method")
if (!is.null(method) && method != "auto") {
dl <- tryCatch(utils::download.file(url, quiet = TRUE,
method = method, ...), error = function(e) NULL)
}
else {
dl <- tryCatch(downloader::download(url, quiet = TRUE,
...), error = function(e) NULL)
}
return(dl)
有意思的是,我看了看,它里面的两个下载函数(utils::download.file和 downloader::download)我都测试了,是可以独立下载的:
utils::download.file( "http://rest.kegg.jp/link/hsa/pathway" ,
tempfile() )
# downloader::download 也可以
trying URL 'http://rest.kegg.jp/link/hsa/pathway'
Content type 'text/plain; charset=utf-8' length unknown
downloaded 807 KB
但是这两个下载函数(utils::download.file和 downloader::download),我是使用默认参数进行下载,都没有设置协议,我看了看函数里面的协议是:
getOption("clusterProfiler.download.method")
[1] "libcurl"
也就是说,它们两个下载函数(utils::download.file和 downloader::download)使用默认的方法,也就是 'auto' 是可以去访问我们的浏览器可以访问的:https://rest.kegg.jp/link/hsa/pathway,的内容。但是如果给它们 "libcurl"的协议,就会失败。
所以,最简单的解决方案就是修改这个默认的协议即可,我给一个解决方案,就是:
install.packages('R.utils')
R.utils::setOption( "clusterProfiler.download.method",'auto' )
果然,接下来的富集分析就非常丝滑了。
出图如下所示:
KEGG数据库没有倒闭, Y叔的 clusterProfiler包也问题不大,我的一个 run_kegg 函数更不可能有问题。仅仅是因为R语言里面的下载文件的函数的协议需要注意,这两个函数两个下载函数(utils::download.file和 downloader::download),都太底层了。初学者确实不容易找到问题所在。