有小伙伴说他要用gage这个包,因为可以选择sigmet这个index,然后得到的结果只有signaling and metabolic pathways,而不会有他不关心的disease pathways。然而也有各种不爽,他最喜欢的还是clusterProfiler,但没办法只做某些pathways。

我发现大家对clusterProfiler有各种误解,各种觉得没办法,我也很无语啊,明明我写了大量的文档,你们偏不看。clusterProfiler啥都可以做,你想做COG,domain这些没有内置支持的富集分析都可以的,因为clusterProfiler是通用的分析工具,啥都能做。

说到gage的pathway index,这其实是他们对pathway有个分类,这个数据就在https://pathview.uncc.edu/data/khier.tsv可以下载到,要支持他还不容易,但我不喜欢把别人的东西打包在自己的包里,所谓走别人的路,让别人无路可走,这可不是什么好主意。所以呢,我不会内置支持的,你们自己玩。

要玩这个,也很容易。无非是拿它的gene sets来做嘛。我们先来看看它分了什么类:

> khier$category %>% table
.
                  Cellular Processes                     Drug Development
                                  33                                   75
Environmental Information Processing       Genetic Information Processing
                                  40                                   22
                      Human Diseases                           Metabolism
                                  83                                  181
                  Organismal Systems
                                  81

总共分了7个,信号通路类的,包括:

  • Genetic Information Processing
  • Environmental Information Processing
  • Cellular Processes
  • Organismal Systems

代谢通路类的,只有自己Metabolism,这也是KEGG最大的一个类别。而所谓的sigmet,就是信号通路+代谢通路。

疾病类的是Human Diseases,而drug development这一类别,gage也有拿出来给大家用。这样大家就知道大概怎么回事了。

那么怎么来和clusterProfiler衔接?gage的分析步骤首先就是准备gene sets,准备完之后,你就可以直接拿来用了,所以说我啥都不用干,本身就是支持的。

> require(gage)
> data(kegg.gs)
> kg.mouse<- kegg.gsets("mouse")
> kegg.gs<- kg.mouse$kg.sets[kg.mouse$sigmet.idx]

下面这个代码,把这个gene sets从list变成data.frame,然后我随机拿100个基因来做演示:

> kegg.gs.df = data.frame(TERM=rep(names(kegg.gs), sapply(kegg.gs, length)), GENE= unlist(kegg.gs))
> de = sample(kegg.gs.df$GENE, 100)

clusterProfiler提供了通用的enricher用于做ORA,和GSEA用于做GSEA,啥都可以搞,不要再以为只有GO和KEGG了,当然你也可以用enricher和GSEA来搞GO和KEGG。

> require(clusterProfiler)
> x = enricher(de, pvalueCutoff=1, qvalueCutoff=1, TERM2GENE= kegg.gs.df)
> head(x[,-8], 2)
                                                                                             ID
mmu04921 Oxytocin signaling pathway                         mmu04921 Oxytocin signaling pathway
mmu04261 Adrenergic signaling in cardiomyocytes mmu04261 Adrenergic signaling in cardiomyocytes
                                                                                    Description
mmu04921 Oxytocin signaling pathway                         mmu04921 Oxytocin signaling pathway
mmu04261 Adrenergic signaling in cardiomyocytes mmu04261 Adrenergic signaling in cardiomyocytes
                                                GeneRatio  BgRatio       pvalue
mmu04921 Oxytocin signaling pathway                 14/99 153/7523 9.784148e-09
mmu04261 Adrenergic signaling in cardiomyocytes     12/99 147/7523 4.245576e-07
                                                    p.adjust       qvalue Count
mmu04921 Oxytocin signaling pathway             1.770931e-06 1.019611e-06    14
mmu04261 Adrenergic signaling in cardiomyocytes 3.842247e-05 2.212169e-05    12