复现Nature图表：GSEA分析及可视化包装函数

❝

详情请联系作者：

❞

这篇帖子主要的目的是写一个转录组GSEA分析和可视化通用的函数。起因是我们想要复现一篇文章的GSEA可视化图片，这个Nature文章GSEA可视化挺好的：

（reference：B-cell-specific checkpoint molecules that regulate anti-tumour immunity）

后来干脆一不做二不休，写一个函数吧，有差异分析结果即可，可视化也是一键出结果。算是懒人福音吧。我们此部分一共有两个函数，一个是KS_GSEA，作用是进行转录组数据GSEA分析，提供clusterProfiler和fgsea两种R包分析方式，KS_GSEA适用于human和mouse两个物种，支持KEGG、GO（BP）的GSEA分析。KS_GSEA_plot作用是进行GSEA结果的可视化，可选择需要的通路进行展示。此函数为付费内容，微信vip群已提前更新，需要购买函数的小伙伴微信联系作者！

接下来，我们看看函数的具体作用和参数解释。本函数使用解说视频下周B站上线！

这里我们以一个单细胞的数据为例，首先做一下差异分析。Bulk RNA的差异分析不再多说。GSEA分析是基于已经完成差异分析结果，且纳入所有基因。目前R做GSEA用的比较多的是clusterProfiler和fgsea包，所以这两种包的分析方式我们都包含进去了。使用的时候选择自己需要的即可。

library(Seurat)library(msigdbr)library(fgsea)library(clusterProfiler)library(ggplot2)#加载函数source('./KS_GSEA.R')source('./KS_GSEA_plot.R')
#差异基因-----------------------------------------------------------------------load("D:/KS项目/公众号文章/GSEA分析及可视化函数/sce_mar.RData")df <- FindMarkers(Macrophage, min.pct = 0, logfc.threshold = 0,                    group.by = "group",ident.1 ="BM",ident.2="GM")df$gene = rownames(df)

得到的差异分析结果是一个数据框，包含gene symbol这一列，还有logFC这一列。我们知道GSEA分析需要对基因按照FC进行排序，我们这里的分析不需要事先进行排序，只要提供gene和logFC即可。两种R包的结果我们都做一下。

#GSEA分析-----------------------------------------------------------------------GSEA_CP <- KS_GSEA(gene = df$gene,                   LogFC = df$avg_log2FC,                   analysis = "KEGG",                   package = 'clusterProfiler',                   OrgDb = 'org.Hs.eg.db')class(GSEA_CP)# [1] "gseaResult"# attr(,"package")# [1] "DOSE"

GSEA_F <- KS_GSEA(gene = df$gene,                  LogFC = df$avg_log2FC,                  analysis = "KEGG",                  package = 'fgsea',                  OrgDb = 'org.Hs.eg.db')class(GSEA_F)# [1] "data.table" "data.frame"

可以看到，clusterProfiler包返回的是一个gseaResult，fgsea包返回的是一个"data.table"，"data.frame"。这些结果就是我们下一步可视化的输入文件。运行结束后，分析结果已txt或者csv的形式直接保存到当前环境路径下！

我们对可视化函数进行了设置，NES>0的结果点用红色显示。NES<0的结果用蓝色点显示。这里我们挑选自己感兴趣的通路进行可视化。

#NES>0p1=KS_GSEA_plot(inputType = "clusterProfiler",             analysis = "KEGG",             data = GSEA_CP,             term = 'Oxidative phosphorylation',             gene = df$gene,             LogFC  = df$avg_log2FC,             OrgDb = 'org.Hs.eg.db')

p2=KS_GSEA_plot(inputType = "fgsea",             analysis = "KEGG",             data = GSEA_F,             term = 'Huntingtons disease',             gene = df$gene,             LogFC = df$avg_log2FC,             OrgDb = 'org.Hs.eg.db')p1+p2

NES<0，这里需要注意一个问题，clusterProfiler和fgsea虽然都是GSEA分析，但是两者得到的结果并不是一摸一样完全相同的，总是有一些出入，这是因为数据库，分析方式不一样。选择需要的包使用一个即可。

#NES<0p3=KS_GSEA_plot(inputType = "clusterProfiler",             analysis = "KEGG",             data = GSEA_CP,             term = 'JAK-STAT signaling pathway',             gene = df$gene,             LogFC  = df$avg_log2FC,             OrgDb = 'org.Hs.eg.db')

p4=KS_GSEA_plot(inputType = "fgsea",             analysis = "KEGG",             data = GSEA_F,             term = 'Jak stat signaling pathway',             gene = df$gene,             LogFC = df$avg_log2FC,             OrgDb = 'org.Hs.eg.db')p3+p4

当然了，我们可以利用循环一次性可视化通路。

#批量循环#NES>0pathway <- c('Oxidative phosphorylation',             "Parkinson disease",             "Biosynthesis of amino acids",             "Cardiac muscle contraction")
pathway_list <- list()for (i in 1:length(pathway)) {  p = KS_GSEA_plot(inputType = "clusterProfiler",                   analysis = "KEGG",                   data = GSEA_CP,                   term = pathway[i],                   gene = df$gene,                   LogFC  = df$avg_log2FC,                   OrgDb = 'org.Hs.eg.db')  pathway_list[[i]] <- p  }#组合图CombinePlots(pathway_list, ncol = 2)

用另一个数据批量做一下下调的。

#NES<0pathway1 <- c('Melanoma',              'Prostate cancer',              'Ecm receptor interaction',              'Tgf beta signaling pathway')
pathway_list1 <- list()for (i in 1:length(pathway)) {  p = KS_GSEA_plot(inputType = "fgsea",                   analysis = "KEGG",                   data = GSEA_F,                   term = pathway1[i],                   gene = df$gene,                   LogFC  = df$avg_log2FC,                   OrgDb = 'org.Hs.eg.db')  pathway_list1[[i]] <- p  }#组合图CombinePlots(pathway_list1, ncol = 2)

这就是所有内容了，希望对你有所启发。这个函数其实并不是很完美，首先是物种只支持小鼠和人，其次是分析参数我没有再设置，用的是我常用的。当然了，这个函数框架我提供了，需要更加个性化的可以自行修改。觉得分享有用的，点个赞呗！