KEGG enrichment analysis with latest online data using clusterProfiler

[This article was first published on YGC » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

KEGG.db is not updated since 2012. The data is now pretty old, but many of the Bioconductor packages still using it for KEGG annotation and enrichment analysis.

As pointed out in ‘Are there too many biological databases‘, there is a problem that many out of date biological databases often don’t get offline. This issue also exists in web-server or software that using out-of-date data. For example, the WEGO web-server stopped updating GO annotation data since 2009, and WEGO still online with many people using it. The biological story may changed totally if using a recently updated data. Seriously, We should keep an eye on this issue.

Now enrichKEGG function is reloaded with a new parameter use.KEGG.db. This parameter is by default setting to FALSE, and enrichKEGG function will download the latest KEGG data for enrichment analysis. If the parameter use.KEGG.db is explicitly setting to TRUE, it will use the KEGG.db which is still supported but not recommended.

With this new feature, supported species is unlimited if only there are KEGG annotations available in KEGG database. You can access the full list of species supported by KEGG via: http://www.genome.jp/kegg/catalog/org_list.html

Now the organism parameter in enrichKEGG should be abbreviation of academic name, for example ‘hsa’ for human and ‘mmu’ for mouse. It accepts any species listed in http://www.genome.jp/kegg/catalog/org_list.html.

In the current release version of clusterProfiler (in Bioconductor 3.0), enrichKEGG supports about 20 species, and the organism parameter accept common name of species, for instance “human” and “mouse”. For these previously supported species, common name is also supported. So that you script is still working with new version of clusterProfiler. For other species, common name is not supported, since I don’t want to maintain such a long mapping list with many species have no common name available and it may also introduce unexpected bugs.

Example 1: Using online KEGG annotation

?View Code RSPLUS
1
2
3
4
5
6
7
8
library(DOSE)
data(geneList)
de <- names(geneList)[geneList > 1]
 
library(clusterProfiler)
kk <- enrichKEGG(de, organism="hsa", pvalueCutoff=0.05, pAdjustMethod="BH", 
                 qvalueCutoff=0.1, readable=TRUE)
head(summary(kk))
> head(summary(kk))
               ID                            Description GeneRatio  BgRatio
hsa04110 hsa04110                             Cell cycle    31/247 124/6861
hsa03030 hsa03030                        DNA replication     9/247  36/6861
hsa04060 hsa04060 Cytokine-cytokine receptor interaction    25/247 265/6861
hsa04114 hsa04114                         Oocyte meiosis    14/247 113/6861
hsa04115 hsa04115                  p53 signaling pathway    10/247  68/6861
hsa04062 hsa04062            Chemokine signaling pathway    18/247 189/6861
               pvalue     p.adjust       qvalue
hsa04110 2.280256e-18 9.349050e-17 5.280593e-17
hsa03030 3.527197e-06 7.230753e-05 4.084123e-05
hsa04060 8.404037e-06 1.148552e-04 6.487326e-05
hsa04114 4.827484e-05 4.948171e-04 2.794859e-04
hsa04115 1.406946e-04 9.801620e-04 5.536217e-04
hsa04062 1.434383e-04 9.801620e-04 5.536217e-04
                                                                                                                                                                              geneID
hsa04110 CDC45/CDC20/CCNB2/CCNA2/CDK1/MAD2L1/TTK/CHEK1/CCNB1/MCM5/PTTG1/MCM2/CDC25A/CDC6/PLK1/BUB1B/ESPL1/CCNE1/ORC6/ORC1/CCNE2/MCM6/MCM4/DBF4/SKP2/CDC25B/BUB1/MYC/PCNA/E2F1/CDKN2A
hsa03030                                                                                                                            MCM5/MCM2/MCM6/MCM4/FEN1/RFC4/PCNA/RNASEH2A/DNA2
hsa04060                           CXCL10/CXCL13/CXCL11/CXCL9/CCL18/IL1R2/CCL8/CXCL3/CCL20/IL12RB2/CXCL8/TNFRSF11A/CCL5/CXCR6/IL2RA/CCR1/CCL2/IL2RG/CCL4/CCR8/CCR7/GDF5/IL24/LTB/IL7
hsa04114                                                                                          CDC20/CCNB2/CDK1/MAD2L1/CALML5/AURKA/CCNB1/PTTG1/PLK1/ESPL1/CCNE1/CCNE2/BUB1/FBXO5
hsa04115                                                                                                               CCNB2/RRM2/CDK1/CHEK1/CCNB1/GTSE1/CCNE1/CCNE2/SERPINB5/CDKN2A
hsa04062                                                                       CXCL10/CXCL13/CXCL11/CXCL9/CCL18/CCL8/CXCL3/CCL20/CXCL8/CCL5/CXCR6/CCR1/STAT1/CCL2/CCL4/HCK/CCR8/CCR7
         Count
hsa04110    31
hsa03030     9
hsa04060    25
hsa04114    14
hsa04115    10
hsa04062    18
>

In the KEGG.db, there are only 5894 human genes annotated. With current online data, the number of annotated gene increase to 6861 as shown above and of course, p-values changed.

User should pay attention to another issue that readable parameter is only available for those species that has an annotation db. For example, for human we use org.Hs.eg.db for mapping gene ID to Symbol.

Example 2: enrichment analysis of species which are not previously supported

Here, I use a gene list of Streptococcus pneumoniae D39 to demonstrate using enrichKEGG with species that are not supported previously.

?View Code RSPLUS
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
> gene
 [1] "SPD_0065" "SPD_0071" "SPD_0293" "SPD_0295" "SPD_0296" "SPD_0297"
 [7] "SPD_0327" "SPD_0426" "SPD_0427" "SPD_0428" "SPD_0559" "SPD_0560"
[13] "SPD_0561" "SPD_0562" "SPD_0580" "SPD_0789" "SPD_1046" "SPD_1047"
[19] "SPD_1048" "SPD_1050" "SPD_1051" "SPD_1052" "SPD_1053" "SPD_1057"
[25] "SPD_1326" "SPD_1432" "SPD_1534" "SPD_1582" "SPD_1612" "SPD_1613"
[31] "SPD_1633" "SPD_1634" "SPD_1648" "SPD_1678" "SPD_1919"
> spdKEGG = enrichKEGG(gene, organism="spd")
> summary(spdKEGG)
               ID                                 Description GeneRatio BgRatio
spd00052 spd00052                        Galactose metabolism     35/35  35/752
spd02060 spd02060             Phosphotransferase system (PTS)     12/35  47/752
spd01100 spd01100                          Metabolic pathways     28/35 341/752
spd00520 spd00520 Amino sugar and nucleotide sugar metabolism      9/35  43/752
               pvalue     p.adjust       qvalue
spd00052 4.961477e-61 2.480739e-60 5.222608e-61
spd02060 2.470177e-07 6.175443e-07 1.300093e-07
spd01100 1.958319e-05 3.263866e-05 6.871296e-06
spd00520 6.534975e-05 8.168718e-05 1.719730e-05
                                                                                                                                                                                                                                                                                                                             geneID
spd00052 SPD_0065/SPD_0071/SPD_0293/SPD_0295/SPD_0296/SPD_0297/SPD_0327/SPD_0426/SPD_0427/SPD_0428/SPD_0559/SPD_0560/SPD_0561/SPD_0562/SPD_0580/SPD_0789/SPD_1046/SPD_1047/SPD_1048/SPD_1050/SPD_1051/SPD_1052/SPD_1053/SPD_1057/SPD_1326/SPD_1432/SPD_1534/SPD_1582/SPD_1612/SPD_1613/SPD_1633/SPD_1634/SPD_1648/SPD_1678/SPD_1919
spd02060                                                                                                                                                                                                                SPD_0293/SPD_0295/SPD_0296/SPD_0297/SPD_0426/SPD_0428/SPD_0559/SPD_0560/SPD_0561/SPD_1047/SPD_1048/SPD_1057
spd01100                                                                SPD_0071/SPD_0426/SPD_0427/SPD_0428/SPD_0559/SPD_0560/SPD_0561/SPD_0562/SPD_0580/SPD_0789/SPD_1046/SPD_1047/SPD_1048/SPD_1050/SPD_1051/SPD_1052/SPD_1053/SPD_1057/SPD_1326/SPD_1432/SPD_1534/SPD_1582/SPD_1612/SPD_1613/SPD_1633/SPD_1634/SPD_1648/SPD_1919
spd00520                                                                                                                                                                                                                                           SPD_0580/SPD_1326/SPD_1432/SPD_1612/SPD_1613/SPD_1633/SPD_1634/SPD_1648/SPD_1919
         Count
spd00052    35
spd02060    12
spd01100    28
spd00520     9

Summary

To summarize, clusterProfiler supports downloading the latest KEGG annotation for enrichment analysis and it supports all species that have KEGG annotation available in KEGG database.

To install the devel version of clusterProfiler, start R and enter the following command:

?View Code RSPLUS
1
2
install.packages(c("DOSE", "clusterProfiler"),
                repos="http://www.bioconductor.org/packages/devel/bioc")

Related Posts

To leave a comment for the author, please follow the link and comment on their blog: YGC » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)