GSEA软件的二次作图

我们知道,用GSEA的java软件分析的结果之后,默认的图是png,而且只有72 dpi,比如下面这样的

image.png

其实颜值也还好,但是很模糊,这是很不方便用于发文章的,那么有没有一些二次作图的方法呢?
答案是肯定的,不仅有,而且还有很多方法,不知道你们有没有发现本地文件夹里其实还有一个edb文件夹,而这个文件夹只有三个文件,一个rnk,一个gmt,另外就是results.edb,实际上所有的图片信息都在results.edb里,在这里展示一些我收集的二次绘图办法。

Rtoolbox二次绘图

Rtoolbox这个包目前只能在Github上安装,里面也只有两个函数replotGSEA()OverviewPlot(),由于Github经常访问不了,所以可以导入到Gitee上了,安装起来也很简单。

remotes::install_git('https://gitee.com/swcyo/Rtoolbox')

而关于这个包的replotGSEA()的函数使用也很简单

replotGSEA(path, gene.set, class.name, metric.range, enrichment.score.range)

比如我的本地结果都在~/my_analysis.GseaPreranked.1664378586466 这个文件夹里,里面有很多结果,我们只需要提取相应的通路名称,设置一些简单的函数就可以出一个pdf的图了,比如我要提取PPAR通路的结果,只需要一个函数即可

library(Rtoolbox) ##加载R包
replotGSEA(path = '~/my_analysis.GseaPreranked.1664378586466', ##设置本地文件夹路径
           gene.set = "KEGG_PPAR_SIGNALING_PATHWAY%HSA03320", ## 提取PPAR通路
           class.name = "PPAR_SIGNALING_PATHWAY", ##定义图中间的名称
           ## enrichment.score.range= c(-1, 1) ###设置富集分数范围,一般默认即可
           )

这时候会再弹出一个R的窗口(Mac系统可能提示要安装Quartz),这时候会显示一个图,显示了一个比自带更好看的图,还能显示p、FDR和NES的值,我们适当的拉伸图片的长宽,然后点Save可以保存为pdf,之后再自己编辑结果,见下图所示。与原始图比较简直就是天壤之别吧。

replotGSEA

然而,这个方案有两个缺陷

  1. 不能在一张图片上设置多条通路
  1. 不能使用代码自由保存图片格式和大小

gseaplot_modified函数二次绘图

使用这个函数其实纯属于不讲武德的方法,完完全全就是调用Rtoolbox这个包的replotGSEA()绘图,唯一的区别就是这个函数不需要安装Rtoolbox这个包,而是直接定义函数,要说区别吧,我仔细对比了一些源代码R/ReplotGSEA.R,也就是在图片的设置上有非常非常非常细微的差距而已。。。

因为两个函数没有本质差异,所以我也就不放结果了,需要的还不如直接复制R/ReplotGSEA.R这个链接里的函数,这没有什么好说的了。。。

基于ggplot2的绘图

这个教程来自于GSEA自定义做图 - 简书 (jianshu.com),当然最好看的肯定是使用clusterProfiler计算好的结果,然后使用enrichplot包的gseaplot2()函数来绘图,当然我们也是可以借鉴Y叔的画图思路要成图,但这个要求太高。这个教程其实还是在Rtoolbox的基础上进行二次修改,将replotGSEA()函数的作图取消,改成单独提取rank和ES的值,然后使用ggplot2拼图。原理无非就是在结果文件夹中有个edb文件夹,里面又有一个.edb 和 .rank文件,这个文件就是做图的原始文件,如果你动手能力强,可以封装成一个函数,也可以自己开发一个包。

使用修改的函数直接提取数据作图。然而对于单个GSEA而已,GSEA的文件夹里还有png图和tsv的表格(很久以前是xls),网上当然也有一些教程,我们可以先看看tsv的结果,比如我们继续使用PPAR通路的表格,可以看到表格里有SYMBOLRANK.IN.GENE.LISTRANK.METRIC.SCORERUNNING.ES等信息。
我们先把表格读进来

data<-read.delim("~/my_analysis.GseaPreranked.1664378586466/KEGG_PPAR_SIGNALING_PATHWAY%HSA03320.tsv")
## 看看数据分布
head(data)
##    NAME SYMBOL RANK.IN.GENE.LIST RANK.METRIC.SCORE RUNNING.ES CORE.ENRICHMENT
## 1 row_0   MMP1               265         0.4223908 0.03071519              No
## 2 row_1    ME1               368         0.3963699 0.06376065              No
## 3 row_2   OLR1               567         0.3652790 0.09122010              No
## 4 row_3   SCD5              1024         0.3145396 0.10664627              No
## 5 row_4    UBC              1589         0.2741413 0.11529607              No
## 6 row_5  FABP5              2359         0.2328998 0.11430056              No
##    X
## 1 NA
## 2 NA
## 3 NA
## 4 NA
## 5 NA
## 6 NA

GSEA二次绘图,主要是三部分拼图,第一部分是曲线,第二部分是网格线,第三部分是底下的曲线,可以使用的办法很多,重点是知道图的x轴和y轴是什么,推荐使用ggplot2画图,当然如果你想省事,用ggpubr更简单

我们先画最上面的图,可以使用geom_line画出,见下面所示。

library(ggplot2)
p1<-ggplot(data) +
  aes(x = RANK.IN.GENE.LIST, y = RUNNING.ES) + #x轴是rank值,y轴是ES值
  geom_line(size = 1, colour = "#f87669") +
  labs( y = "Enrichment score (ES)",  title = "PPAR SIGNALING PATHWAY",x=NULL) +
  theme_bw(base_size = 12)+ #设置主题和字体大小
  theme(axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(),## 将x轴文字清空
        plot.title=element_text(hjust=0.5))+ #设置标题居中
  scale_x_continuous(expand = c(0, 0)) + #取消x轴左右的空白
  geom_hline(yintercept = 0, linetype = "dashed") #添加ES为0的基准线
p1

GSEA上部分的曲线图

接着我们画中间的部分,见下图所示。

p2<-ggplot(data, aes(x = RANK.IN.GENE.LIST)) +
    geom_linerange(aes(ymin=-min(RANK.IN.GENE.LIST), ymax=max(RANK.IN.GENE.LIST))) +
    xlab(NULL) + ylab(NULL) + theme_bw()+
    theme(legend.position = "none",
          plot.margin = margin(t=-.1, b=0,unit="cm"),
          axis.ticks = element_blank(),
          axis.text = element_blank(),
          axis.line.x = element_blank()) +
    scale_x_continuous(expand=c(0,0)) +
    scale_y_continuous(expand=c(0,0))
p2
GSEA中间部分

最后下面的rank部分,见下图所示。

p3<-ggplot(data) +
  aes(x = RANK.IN.GENE.LIST, y = RANK.METRIC.SCORE) +
  geom_area(size = 1.5,fill='gray30') +
  theme_bw(base_size = 12)+ ylab("Ranked List Metric")+xlab("Rank in Ordered Dataset") +
  scale_x_continuous(expand=c(0,0))

GSEA下面部分

最后,将三张图拼成一张图即可,见最终所示。

library(patchwork)

p1/p2/p3+plot_layout(heights = c(0.5,0.2,0.3))
ggplot2-GSEA

最后,三张图其实都是一样的


最后的全家福
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 217,406评论 6 503
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,732评论 3 393
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 163,711评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,380评论 1 293
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,432评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,301评论 1 301
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,145评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,008评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,443评论 1 314
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,649评论 3 334
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,795评论 1 347
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,501评论 5 345
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,119评论 3 328
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,731评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,865评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,899评论 2 370
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,724评论 2 354

推荐阅读更多精彩内容