听说你只有富集表格还想画 GSEA 图?

1引言

首先感谢使用小伙伴的支持和点赞!

我们知道绘制图形是需要对应存储数据的对象的, clusterProfiler 的 gseGO/gseKEGG 返回一个gseaResult 对象,你可以使用 gseaplot/gseaplot2 来绘制图形。有时候你可能输出完富集的表格后,并没有保存这个对象,然后过段时间回头来又想绘图,但是只有表格却画不了,于是你又重新跑了一次富集分析,但发现和上次富集的结果又不一样了! 难搞。

我一般是一起保存富集的 gseaResult 对象和富集的表格,这样下次就可以继续绘图了,保证两次结果一致。要是能把富集表格再次转化为 gseaResult 岂不是很不错。今天群友分享了一个相关的代码,不过不好用,我们可以根据自己的需求,自己手撸一个这样的代码来做个转换。

2安装

# install.packages("devtools")
devtools::install_github("junjunlab/GseaVis")

3测试

我们先做个简单富集分析:

library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(GseaVis)

# load test data
data(geneList, package="DOSE")

# check
head(geneList)
# 4312     8318    10874    55143    55388      991
# 4.572613 4.514594 4.418218 4.144075 3.876258 3.677857

# KEGG enrich
kk <- gseKEGG(geneList     = geneList,
              organism     = 'hsa',
              minGSSize    = 120,
              pvalueCutoff = 0.05,
              verbose      = FALSE)

go <- gseGO(geneList = geneList,
            ont = "BP",
            keyType = "ENTREZID",
            OrgDb = org.Hs.eg.db,
            pvalueCutoff = 0.05)

KK_df <- data.frame(kk)

# check
head(KK_df[1:3,1:5])
#                ID                             Description setSize enrichmentScore      NES
# hsa04110 hsa04110                              Cell cycle     139       0.6637551 2.827249
# hsa05169 hsa05169            Epstein-Barr virus infection     193       0.4335010 1.927343
# hsa05166 hsa05166 Human T-cell leukemia virus 1 infection     202       0.3893613 1.748318

go_df <- data.frame(go)

# check
head(go_df[1:3,1:5])
#                    ID                          Description setSize enrichmentScore      NES
# GO:0000070 GO:0000070 mitotic sister chromatid segregation     151       0.6795449 2.956225
# GO:0000819 GO:0000819         sister chromatid segregation     184       0.6542189 2.891093
# GO:0098813 GO:0098813       nuclear chromosome segregation     236       0.6307606 2.868515

4dfGO2gseaResult

dfGO2gseaResult 将 gseGO 输出的表格转为 gseaResult 对象。

我们先看看 go 对象:

转换 (你需要提供对应的富集表格和排序的 geneList,还有注释文件):

test_go <- dfGO2gseaResult(enrich.df = go_df,
                           geneList = geneList,
                           OrgDb = org.Hs.eg.db)

对比两个对象的绘图结果(可以看到是一样的):

library(patchwork)

pnew <- dotplot(test_go)
praw <- dotplot(go)

plot_grid(praw,pnew)

绘制 gsea 图对比:

gseaplot2(test_go,geneSetID = "GO:0000070")
gseaplot2(go,geneSetID = "GO:0000070")

用 gseaVis 画个看看:

praw <- gseaNb(object = go,geneSetID = "GO:0000070")
pnew <- gseaNb(object = test_go,geneSetID = "GO:0000070")

plot_grid(praw,pnew)

5dfKEGG2gseaResult

dfKEGG2gseaResult 将 gseGO 输出的表格转为 gseaResult 对象。

看一下 kk 对象:

转换:

test_KK <- dfKEGG2gseaResult(enrich.df = KK_df,
                             geneList = geneList,
                             organism = "hsa")

画图对比:

pnew <- dotplot(test_KK)
praw <- dotplot(kk)

plot_grid(praw,pnew)

praw <- gseaNb(object = kk,geneSetID = "hsa04151")
pnew <- gseaNb(object = test_KK,geneSetID = "hsa04151")

plot_grid(praw,pnew)

6其它

上面我们看到对应的基因是 ENTREZID,我们也可以在转换的对象上应用相关的函数,比如 ID 转换:

test_KK2 <- setReadable(test_KK,
                        OrgDb = "org.Hs.eg.db",
                        keyType = "ENTREZID")

7结尾

路漫漫其修远兮,吾将上下而求索。

欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 (微信交流群需收取 20 元入群费用,一旦交费,拒不退还!(防止骗子和便于管理)) 。QQ 群可免费加入, 记得进群按格式修改备注哦。