今天在做富集分析的時(shí)候遇到下面的一個(gè)報(bào)錯(cuò)：
wrong orderBy parameter; set to default
orderBy = "x"
Error in DOSE::setReadable(kk, OrgDb = "org.Hs.eg.db", keytype = "ENTREZID") :
參數(shù)沒有用(keytype = "ENTREZID")
先理解一下這個(gè)參數(shù)沒有用份乒，到底是沒有用這個(gè)keytype = "ENTREZID"參數(shù)呢管呵，還是用了這個(gè)keytype = "ENTREZID"這個(gè)參數(shù)但是用的不對(duì)所以沒有用呢兴垦？感嘆中國語言的博大精深日矫！然后就犯了一個(gè)大家經(jīng)常都會(huì)時(shí)不時(shí)犯的錯(cuò)誤茵汰，就是沒經(jīng)過思索就去問別人枢里。額，我就問了老大蹂午。老大解釋這個(gè)沒有用是你這個(gè)參數(shù)用的不對(duì)而不是你沒有用這個(gè)參數(shù)栏豺。那我問好?一下這個(gè)setReadable，得到的是如下圖
但是我依然沒有看出明道画侣，還在將之前的kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID') 改成kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='auto')冰悠。又錯(cuò)了！最后還是在老大的提醒下配乱，知道其實(shí)就是包被作者更改了溉卓，改之前是keytype，改之后是keyType搬泥，就是小寫的t變成了大寫的T桑寨，頓感的我根本沒看出來，并且實(shí)際上auto就是ENTREZID忿檩，如果不知道這個(gè)的話就還要去查一下背景了尉尾。
果然，我把keytype改成keyType后燥透，就不報(bào)錯(cuò)了沙咏，不過出現(xiàn)下面的一個(gè)wrong的提示，當(dāng)然不用管班套，一樣出圖肢藐，也可以自己理解一下這個(gè)提示的意思，wrong的提示如下：
wrong orderBy parameter; set to default
orderBy = "x"
wrong orderBy parameter; set to defaultorderBy = "x"
wrong orderBy parameter; set to defaultorderBy = "x"
Saving 7 x 7 in image
接下來吱韭，我是想看看這個(gè)我是想知道setReadable這個(gè)代碼的意思吆豹，畢竟之前總是拿來老大的代碼直接做富集分析了，之前也沒遇見這次的報(bào)錯(cuò)，所以還是想找到之前注釋的代碼痘煤，研究一下看看凑阶，之前的代碼如下
## KEGG pathway analysis
### 做KEGG數(shù)據(jù)集超幾何分布檢驗(yàn)分析，重點(diǎn)在結(jié)果的可視化及生物學(xué)意義的理解衷快。
if(T){
### over-representation test
kk.up <- enrichKEGG(gene = gene_up,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.up)[,1:6]
dotplot(kk.up );ggsave('kk.up.dotplot.png')
kk.down <- enrichKEGG(gene = gene_down,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.down)[,1:6]
dotplot(kk.down );ggsave('kk.down.dotplot.png')
kk.diff <- enrichKEGG(gene = gene_diff,
organism = 'hsa',
pvalueCutoff = 0.05)
head(kk.diff)[,1:6]
dotplot(kk.diff );ggsave('kk.diff.dotplot.png')
kegg_diff_dt <- as.data.frame(kk.diff)
kegg_down_dt <- as.data.frame(kk.down)
kegg_up_dt <- as.data.frame(kk.up)
down_kegg<-kegg_down_dt[kegg_down_dt$pvalue<0.05,];down_kegg$group=-1
up_kegg<-kegg_up_dt[kegg_up_dt$pvalue<0.05,];up_kegg$group=1
source('functions.R')
g_kegg=kegg_plot(up_kegg,down_kegg)
print(g_kegg)
ggsave(g_kegg,filename = 'kegg_up_down.png')
### GSEA
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 120,
pvalueCutoff = 0.9,
verbose = FALSE)
head(kk_gse)[,1:6]
gseaplot(kk_gse, geneSetID = rownames(kk_gse[1,]))
down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
g_kegg=kegg_plot(up_kegg,down_kegg)
print(g_kegg)
ggsave(g_kegg,filename = 'kegg_up_down_gsea.png')
}
但是這次有setReadable的代碼宙橱，是如下
## KEGG pathway analysis
### 做KEGG數(shù)據(jù)集超幾何分布檢驗(yàn)分析，重點(diǎn)在結(jié)果的可視化及生物學(xué)意義的理解蘸拔。
run_kegg <- function(gene_up,gene_down,geneList=F,pro='test'){
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 = 'hsa',
#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.Hs.eg.db',keyType='ENTREZID')
write.csv(kk@result,paste0(pro,'_kk.up.csv'))
# 首先是下調(diào)基因集。
kk.down <- enrichKEGG(gene = gene_down,
organism = 'hsa',
#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.Hs.eg.db',keyType='ENTREZID')
write.csv(kk@result,paste0(pro,'_kk.down.csv'))
# 最后是上下調(diào)合并后的基因集都伪。
kk.diff <- enrichKEGG(gene = gene_diff,
organism = 'hsa',
pvalueCutoff = 0.05)
head(kk.diff)[,1:6]
kk=kk.diff
dotplot(kk)
kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
write.csv(kk@result,paste0(pro,'_kk.diff.csv'))
kegg_diff_dt <- as.data.frame(kk.diff)
kegg_down_dt <- as.data.frame(kk.down)
kegg_up_dt <- as.data.frame(kk.up)
down_kegg<-kegg_down_dt[kegg_down_dt$pvalue<0.01,];down_kegg$group=-1
up_kegg<-kegg_up_dt[kegg_up_dt$pvalue<0.01,];up_kegg$group=1
#畫圖設(shè)置, 這個(gè)圖很丑呕乎，大家可以自行修改。
g_kegg=kegg_plot(up_kegg,down_kegg)
print(g_kegg)
ggsave(g_kegg,filename = paste0(pro,'_kegg_up_down.png') )
if(geneList){
### GSEA
## GSEA算法跟上面的使用差異基因集做超幾何分布檢驗(yàn)不一樣陨晶。
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 20,
pvalueCutoff = 0.9,
verbose = FALSE)
head(kk_gse)[,1:6]
gseaplot(kk_gse, geneSetID = rownames(kk_gse[1,]))
gseaplot(kk_gse, 'hsa04110',title = 'Cell cycle')
kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
tmp=kk@result
write.csv(kk@result,paste0(pro,'_kegg.gsea.csv'))
# 這里找不到顯著下調(diào)的通路猬仁，可以選擇調(diào)整閾值，或者其它先誉。
down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
g_kegg=kegg_plot(up_kegg,down_kegg)
print(g_kegg)
ggsave(g_kegg,filename = paste0(pro,'_kegg_gsea.png'))
#
}
}
上面代碼最主要的代碼是什么呢湿刽？沒錯(cuò)，就是每次注釋以后多了一個(gè)kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')褐耳，還是研究一下這個(gè)setReadable函數(shù)诈闺。
參考https://yulab-smu.github.io/clusterProfiler-book/chapter14.html#setReadable
定位到
setReadable函數(shù)，如下圖
啥意思呢铃芦？看下面兩個(gè)截圖就懂了
用setReadable進(jìn)行轉(zhuǎn)換雅镊，注意下面的圖哦，這里面其實(shí)局勢(shì)keyType了刃滓！同時(shí)geneID與上圖的差別就是仁烹，是我們可識(shí)別的gene symbol了，而不是上圖中在做富集分析前需要轉(zhuǎn)換的ENTREZID了咧虎，就是那些我不認(rèn)識(shí)他們卓缰，他們也不認(rèn)識(shí)我的數(shù)字了！
好砰诵！那我就趕緊試試吧征唬，在現(xiàn)在的代碼中，去掉那句setReadable代碼看看茁彭，然后再加上看看总寒，得到的富集分析的結(jié)果有什么不同吧！下面上圖是沒有用setReadable函數(shù)進(jìn)行轉(zhuǎn)換為基因名尉间，下圖是用了setReadable函數(shù)轉(zhuǎn)換為基因名偿乖。
通過上面嘗試，我突然想起來哲嘲，我之前寫過的一個(gè)小推文贪薪，就是之前我也是得到了富集分析后的結(jié)果，然后想要獲得基因名眠副，見：KEGG和GO富集分析獲取通路中的基因画切，當(dāng)時(shí)可是費(fèi)了九牛二虎之力，才搞出來囱怕，其實(shí)原來一個(gè)函數(shù)setReadable就能搞定了霍弹！
總結(jié)：多看，多做娃弓，多思考典格，多嘗試，一定會(huì)有收獲台丛！重要的是耍缴，一定要有個(gè)關(guān)鍵時(shí)刻能給你指點(diǎn)迷津的好老大！
