当clusterProfiler遇见stringdb...

上次搞的ProjectYulab，还是有小伙伴响应的。当前最好的一个案例就是我说支持PPI网络的，直接上截图。

那么在当前github版本的clusterProfiler中，我们把这块功能给加进去了。直接上实例吧。

require(aplot)
require(ggplot2)
require(clusterProfiler)

genes <- c("TP53", "BRCA1", "BRCA2", "S100A6")
g1 <- getPPI(genes, taxID="9606")
plot(g1)

一条plot指令就出图了。当然这只是igraph默认的画图，我们留待以后再加强一下画图的功能。

新加入的getPPI函数，你的输入只需要基因向量，以及相应的物种ID即可。

而这个物种ID，我们也搞了一个函数getTaxID来帮助你获取。

> getTaxID("Homo sapiens")
[1] "9606"

那么你肯定要问，你想把一些相互作用的蛋白也纳入进去呢？这有何难，用参数add_nodes指令数量就可以了。

g2 <- getPPI(genes, taxID='9606', add_nodes=10)
plot(g2)

---

接下来，你肯定要问，那和富集分析有什么关系？

我们来跑个例子吧。

data(geneList, package="DOSE")
de = names(geneList)[1:100]

x = enrichKEGG(de)      
y = setReadable(x, 'org.Hs.eg.db', 'ENTREZID')

然后我们拿第一条富集的结果出来看。

> y[1,]
               ID Description GeneRatio  BgRatio       pvalue    p.adjust
hsa04110 hsa04110  Cell cycle     12/58 157/8398 5.077653e-10 6.14396e-08
               qvalue
hsa04110 5.719041e-08
                                                                        geneID
hsa04110 CDC45/CDC20/CCNB2/NDC80/CCNA2/CDK1/MAD2L1/CDT1/TTK/AURKB/CHEK1/TRIP13
         Count
hsa04110    12

假如我就想看看富集于某个通路的基因其编码的蛋白质之间是个什么样的相互作用，那么现在就很简单了。富集分析的结果直接可以做为getPPI的输入，然后可以传一下富集分析中的第几条通路，是想要看的。

g = getPPI(y, 1)
plot(g)

要加蛋白进入，同样指定一下数目，就加进去了。

g2 = getPPI(y, 1, add_nodes=10)
plot(g2)

那么你还想问，能不能几条通路的基因一起进去做PPI网络，这当然可以。

聪明的你可能已经发现，为什么这些PPI网络都这么高度连通？因为stringdb就是这样，它默认是会把功能相关的（比如同表达）也做给连上边，如果我们就只要相互作用，那么边会少一些的，这个可能通过参数设定。

g4 = getPPI(y, 1:3, network_type = "physical")
plot(g4)

最后的最后，我们再配合一个aplot，那么我们就可以拼出来不同的通路相应的PPI了。

pp <- list(
    ~plot(getPPI(y, 1, network_type = "physical")),
    ~plot(getPPI(y, 2, network_type = "physical"))
)
names(pp) <- y$Description[1:2]

pg <- plot_list(gglist = pp, ncol=1)
pg

当然需要强调一下的是，我们还可以自己搞点基于ggplot2画网络图的函数，然后画图会更好一些。感觉是应该这样搞一下，以后网络图就可以成为了大家富集分析结果的一张子图了。

好了，再回到开头，然后可以有更多的小伙伴加入到ProjectYulab中来，练练手还是不错的。