我們知道凰慈，用GSEA的java軟件分析的結(jié)果之后，默認的圖是png驼鹅，而且只有72 dpi微谓，比如下面這樣的
其實顏值也還好，但是很模糊输钩，這是很不方便用于發(fā)文章的豺型，那么有沒有一些二次作圖的方法呢？
答案是肯定的买乃，不僅有姻氨，而且還有很多方法，不知道你們有沒有發(fā)現(xiàn)本地文件夾里其實還有一個edb文件夾剪验，而這個文件夾只有三個文件肴焊，一個rnk，一個gmt功戚，另外就是results.edb娶眷，實際上所有的圖片信息都在results.edb里，在這里展示一些我收集的二次繪圖辦法疫铜。
Rtoolbox二次繪圖
Rtoolbox這個包目前只能在Github上安裝茂浮，里面也只有兩個函數(shù)replotGSEA()和OverviewPlot()，由于Github經(jīng)常訪問不了壳咕，所以可以導(dǎo)入到Gitee上了席揽，安裝起來也很簡單。
remotes::install_git('https://gitee.com/swcyo/Rtoolbox')
而關(guān)于這個包的replotGSEA()的函數(shù)使用也很簡單
replotGSEA(path, gene.set, class.name, metric.range, enrichment.score.range)
比如我的本地結(jié)果都在~/my_analysis.GseaPreranked.1664378586466 這個文件夾里谓厘，里面有很多結(jié)果幌羞，我們只需要提取相應(yīng)的通路名稱，設(shè)置一些簡單的函數(shù)就可以出一個pdf的圖了竟稳，比如我要提取PPAR通路的結(jié)果属桦，只需要一個函數(shù)即可
library(Rtoolbox) ##加載R包
replotGSEA(path = '~/my_analysis.GseaPreranked.1664378586466', ##設(shè)置本地文件夾路徑
gene.set = "KEGG_PPAR_SIGNALING_PATHWAY%HSA03320", ## 提取PPAR通路
class.name = "PPAR_SIGNALING_PATHWAY", ##定義圖中間的名稱
## enrichment.score.range= c(-1, 1) ###設(shè)置富集分數(shù)范圍熊痴，一般默認即可
)
這時候會再彈出一個R的窗口（Mac系統(tǒng)可能提示要安裝Quartz），這時候會顯示一個圖聂宾，顯示了一個比自帶更好看的圖果善，還能顯示p、FDR和NES的值系谐，我們適當(dāng)?shù)睦靾D片的長寬巾陕，然后點Save可以保存為pdf，之后再自己編輯結(jié)果纪他，見下圖所示鄙煤。與原始圖比較簡直就是天壤之別吧。
然而茶袒，這個方案有兩個缺陷
- 不能在一張圖片上設(shè)置多條通路
- 不能使用代碼自由保存圖片格式和大小
gseaplot_modified函數(shù)二次繪圖
使用這個函數(shù)其實純屬于不講武德的方法梯刚，完完全全就是調(diào)用Rtoolbox這個包的replotGSEA()繪圖，唯一的區(qū)別就是這個函數(shù)不需要安裝Rtoolbox這個包薪寓，而是直接定義函數(shù)亡资，要說區(qū)別吧，我仔細對比了一些源代碼R/ReplotGSEA.R预愤，也就是在圖片的設(shè)置上有非常非常非常細微的差距而已沟于。。植康。
因為兩個函數(shù)沒有本質(zhì)差異，所以我也就不放結(jié)果了展懈，需要的還不如直接復(fù)制R/ReplotGSEA.R這個鏈接里的函數(shù)销睁，這沒有什么好說的了。存崖。冻记。
基于ggplot2的繪圖
這個教程來自于GSEA自定義做圖 - 簡書 (jianshu.com)，當(dāng)然最好看的肯定是使用clusterProfiler計算好的結(jié)果来惧，然后使用enrichplot包的gseaplot2()函數(shù)來繪圖冗栗，當(dāng)然我們也是可以借鑒Y叔的畫圖思路要成圖，但這個要求太高供搀。這個教程其實還是在Rtoolbox的基礎(chǔ)上進行二次修改隅居，將replotGSEA()函數(shù)的作圖取消，改成單獨提取rank和ES的值葛虐，然后使用ggplot2拼圖胎源。原理無非就是在結(jié)果文件夾中有個edb文件夾，里面又有一個.edb 和 .rank文件屿脐，這個文件就是做圖的原始文件涕蚤，如果你動手能力強宪卿，可以封裝成一個函數(shù)，也可以自己開發(fā)一個包万栅。
使用修改的函數(shù)直接提取數(shù)據(jù)作圖佑钾。然而對于單個GSEA而已，GSEA的文件夾里還有png圖和tsv的表格（很久以前是xls）烦粒，網(wǎng)上當(dāng)然也有一些教程休溶，我們可以先看看tsv的結(jié)果，比如我們繼續(xù)使用PPAR通路的表格撒遣，可以看到表格里有SYMBOL邮偎，RANK.IN.GENE.LIST，RANK.METRIC.SCORE义黎，RUNNING.ES等信息禾进。
我們先把表格讀進來
data<-read.delim("~/my_analysis.GseaPreranked.1664378586466/KEGG_PPAR_SIGNALING_PATHWAY%HSA03320.tsv")
## 看看數(shù)據(jù)分布
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二次繪圖，主要是三部分拼圖廉涕，第一部分是曲線泻云，第二部分是網(wǎng)格線，第三部分是底下的曲線狐蜕，可以使用的辦法很多宠纯，重點是知道圖的x軸和y軸是什么，推薦使用ggplot2畫圖层释，當(dāng)然如果你想省事婆瓜，用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)+ #設(shè)置主題和字體大小
theme(axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank(),## 將x軸文字清空
plot.title=element_text(hjust=0.5))+ #設(shè)置標題居中
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
最后下面的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))
最后聋溜，三張圖其實都是一樣的
