中文字幕理论片,69视频免费在线观看,亚洲成人app,国产1级毛片,刘涛最大尺度戏视频,欧美亚洲美女视频,2021韩国美女仙女屋vip视频

最近，有粉絲運(yùn)行了我以前的數(shù)據(jù)挖掘成套代碼里面的 run_kegg 函數(shù)，如下所示：

library(clusterProfiler) 
run_kegg(gene_up,gene_down,pro='comp1')

出現(xiàn)了如下所示的報(bào)錯(cuò)：

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)

然后就找我，以為是我們的標(biāo)準(zhǔn)代碼有問題，實(shí)際上我的 run_kegg 函數(shù)僅僅是包裝了 Y叔的 clusterProfiler包而已 ，實(shí)際上里面沒有啥玄機(jī)，如下所示：

## KEGG pathway analysis
### 做KEGG數(shù)據(jù)集超幾何分布檢驗(yàn)分析，重點(diǎn)在結(jié)果的可視化及生物學(xué)意義的理解。
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個(gè)基因集分開做超幾何分布檢驗(yàn)
  # 首先是上調(diào)基因集。
  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'))
  
  # 首先是下調(diào)基因集。
  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'))
  
  # 最后是上下調(diào)合并后的基因集。
  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')) 
}

就是針對(duì)上下調(diào)基因，以及其合并的差異基因，分別獨(dú)立走 enrichKEGG 函數(shù)而已，它來(lái)自于 Y叔的 clusterProfiler包。

實(shí)際上，我們根據(jù)前面的報(bào)錯(cuò)，進(jìn)入  https://www.genome.jp/kegg/catalog/org_list.html 可以很清楚的看到，我們的人和鼠，都是在里面的：

hsa   Homo sapiens (human)
mmu   Mus musculus (house mouse)

所以不可能是物種問題，而且這個(gè)kegg數(shù)據(jù)庫(kù)的官網(wǎng)也可以訪問，那么合理的推測(cè)，是不是Y叔的 clusterProfiler包有問題呢？其實(shí)仔細(xì)看報(bào)錯(cuò)提示是 download.KEGG.Path 函數(shù)問題，也是Y叔的包里面的。這個(gè)download.KEGG.Path 函數(shù)的代碼是：

 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))

也就是說(shuō)在kegg_link函數(shù)就開始報(bào)錯(cuò)了，繼續(xù)看kegg_link的代碼：

url <- paste0("http://rest.kegg.jp/link/", target_db, 
  "/", source_db, collapse = "")
 kegg_rest(url)

實(shí)際上， 它就根據(jù)我們的物種，拼湊了一個(gè) 網(wǎng)頁(yè)，比如：https://rest.kegg.jp/link/hsa/pathway，可以使用瀏覽器打開，能看到里面的通路的id和基因id的對(duì)應(yīng)關(guān)系。所以真正的問題是 kegg_rest 這個(gè)函數(shù)，沒有辦法像瀏覽器那樣直接獲取網(wǎng)頁(yè)里面的內(nèi)容。我們繼續(xù)看kegg_rest 這個(gè)函數(shù)的源代碼：

 dl <- mydownload(rest_url, destfile = f)
 if (is.null(dl)) {
  message("fail to download KEGG data...")
  return(NULL)
 }

真氣人，就是一個(gè)套娃， 里面又是mydownload函數(shù)，那我們繼續(xù)看mydownload函數(shù)的源代碼 ：

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)

有意思的是，我看了看，它里面的兩個(gè)下載函數(shù)（utils::download.file和  downloader::download）我都測(cè)試了，是可以獨(dú)立下載的：

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

但是這兩個(gè)下載函數(shù)（utils::download.file和  downloader::download），我是使用默認(rèn)參數(shù)進(jìn)行下載，都沒有設(shè)置協(xié)議，我看了看函數(shù)里面的協(xié)議是：

 getOption("clusterProfiler.download.method")
[1] "libcurl"

也就是說(shuō)，它們兩個(gè)下載函數(shù)（utils::download.file和  downloader::download）使用默認(rèn)的方法，也就是 'auto'  是可以去訪問我們的瀏覽器可以訪問的：https://rest.kegg.jp/link/hsa/pathway，的內(nèi)容。但是如果給它們 "libcurl"的協(xié)議，就會(huì)失敗。

所以，最簡(jiǎn)單的解決方案就是修改這個(gè)默認(rèn)的協(xié)議即可，我給一個(gè)解決方案，就是：

install.packages('R.utils')
R.utils::setOption( "clusterProfiler.download.method",'auto' )

果然，接下來(lái)的富集分析就非常絲滑了。

出圖如下所示：

KEGG數(shù)據(jù)庫(kù)沒有倒閉， Y叔的 clusterProfiler包也問題不大，我的一個(gè) run_kegg 函數(shù)更不可能有問題。僅僅是因?yàn)镽語(yǔ)言里面的下載文件的函數(shù)的協(xié)議需要注意，這兩個(gè)函數(shù)兩個(gè)下載函數(shù)（utils::download.file和  downloader::download），都太底層了。初學(xué)者確實(shí)不容易找到問題所在。