0 收听	12 听众	360 主题

GSEA分析及结果图

发表于 2025-3-23 09:48:47 | 查看: 52| 回复: 1

背景

用R代码画出PMID: 27789795文章里四个通路画在一起的GSEA结果图

使用场景

在一个图的空间展示富集的多个通路的GSEA结果。

输入数据的预处理

如果你已经做完了GSEA，也建议你用clusterProfiler再做一次GSEA，因为：

富集分析需要不断地维护，更新知识库；不像其它软件，写好了可以用几年。富集分析的软件非常多，但多半不能用，而多半软件是不维护的，这就是问题之所在。如果你使用Broad Institute的GSEA软件，要小心了，它的KEGG注释还停留在2012年。

所以这里我们使用 clusterProfiler，它解决了三个问题：

一是支持大量的非模式生物；

二是支持多种知识库；

三是支持在线数据。

这是大量同类软件的短板，很多只支持GO/KEGG，只支持少量模式生物，基本不更新维护，所以这里我们使用 clusterProfiler，而不用别的软件的结果，直接以实验结果（全基因组表达量）做为输入，衔接 clusterProfiler，最后出图，完事。

做GSEA，需要一个排序文件，此处输入数据是基因组上所有基因的变化倍数：easy_input_rnk.txt，等同于java版本GSEA的rnk文件。

参考《听说你有RNAseq数据却不知道怎么跑GSEA》

library(clusterProfiler)
df = read.table("easy_input_rnk.txt",header = T,as.is = T)
dim(df) # 12444     2
head(df)

df前几行如下：

此处第一列是基因名，需要先转换成ENTREZID。

如果你的基因ID是其他ID，也可以用clusterProfiler转换，TCGA也用它来做转换。

df.id<-bitr(df$SYMBOL, fromType = "SYMBOL", toType = "ENTREZID",OrgDb = "org.Hs.eg.db")
head(df.id)
dim(df.id)# 12145     2
#让基因名、ENTREZID、foldchange对应起来
easy.df<-merge(df,df.id,by="SYMBOL",all=F)
#按照foldchange排序
sortdf<-easy.df[order(easy.df$foldChange, decreasing = T),]
head(sortdf)
gene.expr = sortdf$foldChange
names(gene.expr) <- sortdf$ENTREZID
head(gene.expr) # 有name的数字型向量

有了这个ENTREZID和foldchange信息的 gene.expr之后，就可以使用 clusterProfiler进行GSEA分析了。

进行GSEA分析

require(enrichplot)
require(clusterProfiler)
kk <- gseKEGG(gene.expr, organism = "hsa")
# Reading KEGG annotation online: "https://rest.kegg.jp/link/hsa/pathway"...
# Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/hsa"...
# using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
# 
# preparing geneSet collections...
# GSEA analysis...
# leading edge analysis...
# done...

这一条语句就做完了KEGG的GSEA分析。

另外，还可以用GO、MSigDB、以及自定义的gene set做富集分析。

GSEA结果保存

按照enrichment score从高到低排序，便于查看富集通路

sortkk<-kk[order(kk$enrichmentScore, decreasing = T),]
write.table(sortkk,"gsea_output.txt",sep = "\t",quote = F,col.names = T,row.names = F)

通过查看gsea_output.txt,选择你想画的pathway，记下通路的term ID。

开始画图

不需要另外准备输入文件，clusterProfiler直接可以用 gseaplot出图：

gseaplot(kk, "hsa04510")

当然你可以改线条颜色：

gseaplot(kk, "hsa04510", color.line='steelblue')

另一种出图方式是模拟Broad institute的GSEA软件的图，可以用 gseaplot2出图：

gseaplot2(kk, "hsa04510")

gseaplot2(kk, "hsa04510", color = "firebrick", rel_heights=c(1, .2, .6))

它还有更外一个功能，就是支持同时展示多个pathways的结果

画出enrichment score排名前4的pathway

gseaplot2(kk, row.names(sortkk)[1:4])

也可以通过查看gsea_output.txt，挑出你想show的pathway，给它以hsa开头的KEGG ID

paths <- c("hsa04510", "hsa04512", "hsa04974", "hsa05410")
gseaplot2(kk, paths)

gseaplot2(kk, paths, color = colorspace::rainbow_hcl(4))

如果你不想同时展示3个subplots，而只想选取其中某一些，也是可以的，通过 subplots参数：

gseaplot2(kk, paths, subplots=1)

gseaplot2(kk, paths, subplots=1:2)

gseaplot2(kk, paths, subplots=c(1,3))

当然还可以直接把pvalue table注释在图片上

gseaplot2(kk, paths, color = colorspace::rainbow_hcl(4), pvalue_table = TRUE)

gseaplot2(kk, paths, pvalue_table = TRUE)

enrichplot提供了多种可视化方式，不单单是上面展示的这些，请移步在线文档https://yulab-smu.top/biomedical-knowledge-mining-book/enrichplot.html#running-score-and-preranked-list-of-gsea-result.