水稻如何做KEGG富集分析
本文介绍了如何对水稻LOC4334374或LOC_Os01g01010.1或Os04t0485300命名方式的基因作KEGG富集分析。本文适用于任何能够在http://www.genome.jp/kegg/catalog/org_list.html 搜索到的物种。
Y叔的clusterProfiler是做富集分析的神器，但是会用并且用好它的人却不多。在群众的要求下，Y叔写了非模式基因GO富集分析：以玉米为例+使用OrgDb，旨在教大家如何学会用Bioconductor的AnnotationHub在线检索获取物种注释数据库，然后利用clusterProfiler搞定GO富集分析。
这一篇在Y叔基础上介绍如何使用clusterProfiler做水稻的KEGG富集分析
enrichKEGG函数介绍
Y叔的clusterProfiler可能是目前KEGG富集分析最好的软件，因为能爬取最新的KEGG在线版数据库，而不是用不再更新的KEGG.db。
本部分内容参照KEGG enrichment analysis with latest online data using clusterProfiler（阅读原文直达）
函数是enrichKEGG
enrichKEGG(gene, organism = "hsa", keyType = "kegg", pvalueCutoff = 0.05,
pAdjustMethod = "BH", universe, minGSSize = 10, maxGSSize = 500,
qvalueCutoff = 0.2, use_internal_data = FALSE)介绍一下主要参数的使用方法：
gene： 基因名，要和keyType对应
organism: 需要参考 http://www.genome.jp/kegg/catalog/org_list.html
keyType: 基因的命名方式， “kegg”, ‘ncbi-geneid’, ‘ncib-proteinid’ and ‘uniprot’选择其一
大部分不会用KEGG做水稻富集分析的人，主要就卡在了上面这三个参数上，不知道应该填什么参数。我以水稻作为案例分析如何确定这三个参数。
确定organisam
orgnisam参数是最好确定的， 而且是能为后面两个参数提供帮助。 打开http://www.genome.jp/kegg/catalog/org_list.html 搜索rice
上图左边的osa和dosa就是要填到organisam里参数，选择其中一个就行，比如说我选择dosa。
确定keytypes
学会如何使用enrichKEGG的一个关键在于理解它是如何获取数据的。Y叔说他是在线爬取数据库，相当于我们在KEGG上手动输入基因名，然后查询。
所以，让我们把enrichKEGG当做浏览器，试出合适的keytypes。
# kegg
> enrichKEGG(gene = "abc", organism = "dosa", keyType = "kegg")
--> No gene can be mapped....
--> Expected input gene ID: Os06t0664200-01,Os02t0534400-01,Os06t0185100-00,Os07t0691100-01,Os04t0485300-01,Os03t0330200-00
--> return NULL...
# ncbi-geneid
> enrichKEGG(gene = "abc", organism = "dosa", keyType = "ncbi-geneid")
Error in KEGG_convert("kegg", keyType, species) :
ncbi-geneid is not supported for dosa ...
# ncbi-proteinid
> enrichKEGG(gene = "abc", organism = "dosa", keyType = "ncbi-proteinid")
Error in KEGG_convert("kegg", keyType, species) :
ncbi-proteinid is not supported for dosa ...
# uniprot
> enrichKEGG(gene = "abc", organism = "dosa", keyType = "uniprot")
Error in KEGG_convert("kegg", keyType, species) :
uniprot is not supported for dosa ...因此，我们需要提供的keyType是kegg，基因名应该类似于Os06t0664200-01.按照同样的套路，对于osa的kegg命名方式应该3131385能够用于在NCBI上查找的gene ID。
转换基因名
一般而言，我们无法直接拿到类似于Os06t0664200-01或者3131385的命名数据，通常拿到的是类似于 LOC4334374或LOC_Os01g01010.1的结果。
对于LOC4334374，我们可以利用Bioconductor上的Annotation查找水稻的Org物种数据库， 然后转换。
library(AnnotationHub)
ah <- AnnotationHub()
query(ah, 'oryza')
oryza <- ah[['AH55775']]
# LOC4334374 对应的keytpe是SYMBOL
> head(keys(oryza, keytype = "SYMBOL"))
[1] "LOC107275246" "LOC107275247" "LOC107275248" "LOC107275249"
[5] "LOC107275250" "LOC107275251"你需要把SYMBOL转成GID，见我写的教程。
对于RGAP水稻的基因编号，如LOC_Os01g01010.1 我们要把它变成Os06t0664200-01 RAP-ID的命名方式，符合dosa的要求。相关数据库到https://shigen.nig.ac.jp/rice/oryzabase/ 下载, 速度比较慢。
$ grep LOC_Os01g01010.1 rice_id_20140620174522.txt
Os01g0100100 Os01t0100100-01 J075199P03 301700 J075199P03 LOC_Os01g01010.1 AK242339写一行脚本把LOC_Os01g01010.1 转为 Os06t0664200-01。（同样适用于把Os06t0664200转为Os06t0664200-01）
cat LOC.txt | xargs -i awk 'BEGIN{FS="\t"} $0 ~/{}/ { print $2}' rice_id_20140620174522.txt > RAP_id.txtKEGG富集分析
这一步其实是最简单的一步，因为我们上面几步已经把参数都准备好了。
library(clusterProfiler)
# 对于GID
kk <- enrichKEGG(gid_list, organism="osa",
keyType = "kegg",
pvalueCutoff=0.05, pAdjustMethod="BH",
qvalueCutoff=0.1)
# 对于RAP ID
kk <- enrichKEGG(rap_list, organism="dosa",
keyType = "kegg",
pvalueCutoff=0.05, pAdjustMethod="BH",
qvalueCutoff=0.1)根据具体情况，把gid_list，rap_list, organism替换成你自己的数据。
由于我不做水稻，所以对于水稻的命名方式实在搞不定， 但是本篇文章的思路在于
明确enrichKEGG需要的输入是什么， 然后提供对应的内容就行了。 同样本文适用于任何能够在http://www.genome.jp/kegg/catalog/org_list.html 搜索到的物种。
至于如何批量转ID，找到对应的命名数据库，下载，然后写个脚本
如果有打赏的话，再写一个续集吧