GseaVis 为你的通路基因添加热图 - junjunlab/GseaVis GitHub Wiki

引言

生信技能树 推文 GSEA 确实搭配热图后更直观易懂,展示 GSEA 富集图形下面添加了基因的热图,推文也提供了代码。有粉丝直接把推文发我了,啥也没说,意思我懂,不就是让我操作一波吗?幸好我有 GseaVis 的架子了,往上面再添几个小部件吧。

文献图形展示:

解释:

下面的热图应该是富集到该通路里的 core enrichment 基因的表达量热图。

安装

重新安装获取新功能:

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

欢迎大家在 github 上留下你的小心心。话说真的没人去点吗?

数据准备

肯定是需要表达量文件和差异基因的,这里我用自己的数据,你们用你们自己的,我就不提供了。表达量文件可以是 tpm/fpkm/rpkm:

library(clusterProfiler)
library(org.Mm.eg.db)
library(GseaVis)

# load diff
diff <- read.csv('diffdata.csv') %>%
  arrange(desc(log2FoldChange))

genelist <- diff$log2FoldChange
names(genelist) <- diff$gene_name

# check
head(genelist,3)
#    Aplnr    Foxf1     Bmp5
# 13.45176 13.35322 12.02845

# load tpm
expr <- read.csv('tpm.csv')

# check
head(expr,3)
#   gene_name  day0_rep1  day0_rep2 day4_rep1 day4_rep2
# 1     Gnai3 155.689704 147.752705 361.38421 339.49937
# 2     Cdc45  68.347917  65.870974  61.79624  55.27016
# 3       H19   5.850431   7.558738 794.35331 637.19978

注意: 表达量文件第一列必须是基因名哦!

开始富集分析:

# GO enrich
ego <- gseGO(geneList     = genelist,
             OrgDb        = org.Mm.eg.db,
             ont          = "BP",
             keyType      = 'SYMBOL',
             minGSSize    = 100,
             maxGSSize    = 500,
             pvalueCutoff = 0.05,
             verbose      = FALSE)

单通路绘图

开始画图,默认画一个,geneExpHt参数添加热图,exp参数添加表达量文件:

# default plot
gseaNb(object = ego,
       geneSetID = 'GO:0003205',
       add.geneExpHt = T,
       exp = expr)

保留上面部分:

# subplot
gseaNb(object = ego,
       geneSetID = 'GO:0003205',
       add.geneExpHt = T,
       exp = expr,
       subPlot = 2)

修改热图颜色:

# change heatmap colors
gseaNb(object = ego,
       geneSetID = 'GO:0003205',
       add.geneExpHt = T,
       exp = expr,
       subPlot = 2,
       exp.col = c('green','black','red'))

修改样本顺序:

# change heatmap sample orders
gseaNb(object = ego,
       geneSetID = 'GO:0003205',
       add.geneExpHt = T,
       exp = expr,
       subPlot = 2,
       exp.col = c('green','black','red'),
       sample.order = c(paste('day4_rep',1:2,sep = ''),
                        paste('day0_rep',1:2,sep = '')))

调整基因标签大小:

# change gene text size
gseaNb(object = ego,
       geneSetID = 'GO:0007368',
       add.geneExpHt = T,
       exp = expr,
       subPlot = 2,
       ght.geneText.size = 12)

默认会对表达量文件进行 scale 操作,你也可以设置为false:

# no scale expression
gseaNb(object = ego,
       geneSetID = 'GO:0007368',
       add.geneExpHt = T,
       exp = expr,
       subPlot = 2,
       ght.geneText.size = 12,
       scale.exp = F)

修改热图相对高度:

# relative height
gseaNb(object = ego,
       geneSetID = 'GO:0007368',
       add.geneExpHt = T,
       exp = expr,
       subPlot = 2,
       ght.geneText.size = 12,
       ght.relHight = 0.2)

新样式 GSEA 可视化:

# new style gsea
gseaNb(object = ego,
       geneSetID = 'GO:0007368',
       add.geneExpHt = T,
       exp = expr,
       subPlot = 2,
       ght.geneText.size = 12,
       newGsea = T)

多通路绘图

当然也可以支持多个通路可视化:

# multiple terms
gseaNb(object = ego,
       geneSetID = c('GO:0007368','GO:0006304'),
       add.geneExpHt = T,
       exp = expr,
       subPlot = 2,
       ght.geneText.size = 12,
       addPval = T,
       pvalX = 1,pvalY = 1)

此时热图是 融合了多个通路的基因, 你还可以对热图 按通路进行分面展示:

# facet heatmap
gseaNb(object = ego,
       geneSetID = c('GO:0007368','GO:0006304'),
       add.geneExpHt = T,
       exp = expr,
       subPlot = 2,
       ght.geneText.size = 12,
       addPval = T,
       pvalX = 1,pvalY = 1,
       ght.facet = T)

调整分面,这里选了两个上下调的通路:

# adjust facet scales
gseaNb(object = ego,
       geneSetID = c('GO:0007368','GO:0006304'),
       add.geneExpHt = T,
       exp = expr,
       subPlot = 2,
       ght.geneText.size = 12,
       addPval = T,
       pvalX = 1,pvalY = 1,
       ght.facet = T,
       ght.facet.scale = "fixed")

有时候通路间会共享一些基因:

# some terms share same genes
gseaNb(object = ego,
       geneSetID = c('GO:0007368','GO:0003205'),
       add.geneExpHt = T,
       exp = expr,
       subPlot = 2,
       ght.geneText.size = 6,
       addPval = T,
       pvalX = 1,pvalY = 1,
       ght.facet = T,
       ght.facet.scale = "fixed")

当然你也可以 调整分面的顺序:

# change facet order
gseaNb(object = ego,
       geneSetID = c('GO:0007368','GO:0003205'),
       add.geneExpHt = T,
       exp = expr,
       subPlot = 2,
       ght.geneText.size = 10,
       addPval = T,
       pvalX = 1,pvalY = 1,
       ght.facet = T,
       termID.order = c('GO:0007368','GO:0003205'))

注意: 这里热图基因的顺序是按 rank 来排序的。也就是 foldChange。

结尾

怎么说,拿来吧你!