clusterProfiler 做 GSEA 富集分析

原创 JunJunLab 老俊俊的生信笔记 2023-06-15

收录于合集 #R包分享 168个

浮躁

1引言

分享一下如何使用 Y 叔的 clusterProfiler 包做 GSEA 富集分析。

2分析

加载内置数据:

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

# 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

names是 entriz id, value是排序的基因差异表达值。

KEGG or GO 富集分析:

# GO enrich
ego3 <- gseGO(geneList     = geneList,
              OrgDb        = org.Hs.eg.db,
              ont          = "BP",
              minGSSize    = 100,
              maxGSSize    = 500,
              pvalueCutoff = 0.05,
              verbose      = FALSE)

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

gseaplot 绘图展示:

# plot
p1 <- gseaplot(ego3, geneSetID = 1,
               by = "runningScore",
               title = ego3$Description[1])

p2 <- gseaplot(ego3, geneSetID = 1,
               by = "preranked",
               title = ego3$Description[1])

p3 <- gseaplot(ego3, geneSetID = 1,
               title = ego3$Description[1])

# combine
cowplot::plot_grid(p1, p2, p3, ncol=2, labels=LETTERS[1:3])

gseaplot2 绘图展示:

# gseaplot2
gseaplot2(ego3, geneSetID = 1,
          title = ego3$Description[1])

多个通路一起展示:

gseaplot2(ego3, geneSetID = 1:3)

展示不同部分:

p1 <- gseaplot2(ego3, geneSetID = 1:3, subplots = 1)
p2 <- gseaplot2(ego3, geneSetID = 1:3, subplots = 1:2)

cowplot::plot_grid(p1, p2, ncol=1, labels=LETTERS[1:2])

gsearank 绘图展示:

gsearank(ego3, 1, title = ego3[1, "Description"])

绘制多个通路:

library(ggplot2)
library(cowplot)

pp <- lapply(1:3, function(i) {
  anno <- ego3[i, c("NES", "pvalue", "p.adjust")]
  lab <- paste0(names(anno), "=",  round(anno, 3), collapse="\n")

  gsearank(ego3, i, ego3[i, 2]) + xlab(NULL) +ylab(NULL) +
    annotate("text", 10000, ego3[i, "enrichmentScore"] * .75, label = lab, hjust=0, vjust=0)
})

plot_grid(plotlist=pp, ncol=1)