查看原文
其他

水稻如何做KEGG富集分析

徐洲更 biobabble 2020-02-05

本文转载自《生信媛》,已获得授权。

本文介绍了如何对水稻LOC4334374LOC_Os01g01010.1Os04t0485300命名方式的基因作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做水稻富集分析的人,主要就卡在了上面这三个参数上,不知道应该填什么参数。我以水稻作为案例分析如何确定这三个参数。

确定organism

organism参数是最好确定的, 而且是能为后面两个参数提供帮助。 打开http://www.genome.jp/kegg/catalog/org_list.html 搜索rice

上图左边的osadosa就是要填到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的命名数据,通常拿到的是类似于 LOC4334374LOC_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=" "}  $0 ~/{}/ { print $2}' rice_id_20140620174522.txt > RAP_id.txt

KEGG富集分析

这一步其实是最简单的一步,因为我们上面几步已经把参数都准备好了。

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,找到对应的命名数据库,下载,然后写个脚本

如果有打赏的话,再写一个续集吧

    您可能也对以下帖子感兴趣

    文章有问题?点此查看未经处理的缓存