今天在做富集分析的时候遇到下面的一个报错：
wrong orderBy parameter; set to default
orderBy = "x"
Error in DOSE::setReadable(kk, OrgDb = "org.Hs.eg.db", keytype = "ENTREZID") :
参数没有用(keytype = "ENTREZID")
先理解一下这个参数没有用，到底是没有用这个keytype = "ENTREZID"参数呢，还是用了这个keytype = "ENTREZID"这个参数但是用的不对所以没有用呢？感叹中国语言的博大精深！然后就犯了一个大家经常都会时不时犯的错误，就是没经过思索就去问别人。额，我就问了老大。老大解释这个没有用是你这个参数用的不对而不是你没有用这个参数。那我问好❓一下这个setReadable，得到的是如下图
但是我依然没有看出明道，还在将之前的kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID') 改成kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='auto')。又错了！最后还是在老大的提醒下，知道其实就是包被作者更改了，改之前是keytype，改之后是keyType，就是小写的t变成了大写的T，顿感的我根本没看出来，并且实际上auto就是ENTREZID，如果不知道这个的话就还要去查一下背景了。
果然，我把keytype改成keyType后，就不报错了，不过出现下面的一个wrong的提示，当然不用管，一样出图，也可以自己理解一下这个提示的意思，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
接下来，我是想看看这个我是想知道setReadable这个代码的意思，毕竟之前总是拿来老大的代码直接做富集分析了，之前也没遇见这次的报错，所以还是想找到之前注释的代码，研究一下看看，之前的代码如下
## KEGG pathway analysis
### 做KEGG数据集超几何分布检验分析，重点在结果的可视化及生物学意义的理解。
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数据集超几何分布检验分析，重点在结果的可视化及生物学意义的理解。
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个基因集分开做超几何分布检验
# 首先是上调基因集。
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'))
# 首先是下调基因集。
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'))
# 最后是上下调合并后的基因集。
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
#画图设置, 这个图很丑，大家可以自行修改。
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算法跟上面的使用差异基因集做超几何分布检验不一样。
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'))
# 这里找不到显著下调的通路，可以选择调整阈值，或者其它。
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'))
#
}
}
上面代码最主要的代码是什么呢？没错，就是每次注释以后多了一个kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')，还是研究一下这个setReadable函数。
参考https://yulab-smu.github.io/clusterProfiler-book/chapter14.html#setReadable
定位到
setReadable函数，如下图
啥意思呢？看下面两个截图就懂了
用setReadable进行转换，注意下面的图哦，这里面其实局势keyType了！同时geneID与上图的差别就是，是我们可识别的gene symbol了，而不是上图中在做富集分析前需要转换的ENTREZID了，就是那些我不认识他们，他们也不认识我的数字了！
好！那我就赶紧试试吧，在现在的代码中，去掉那句setReadable代码看看，然后再加上看看，得到的富集分析的结果有什么不同吧！下面上图是没有用setReadable函数进行转换为基因名，下图是用了setReadable函数转换为基因名。
通过上面尝试，我突然想起来，我之前写过的一个小推文，就是之前我也是得到了富集分析后的结果，然后想要获得基因名，见：KEGG和GO富集分析获取通路中的基因，当时可是费了九牛二虎之力，才搞出来，其实原来一个函数setReadable就能搞定了！
总结：多看，多做，多思考，多尝试，一定会有收获！重要的是，一定要有个关键时刻能给你指点迷津的好老大！
