❝详情请联系作者: ❞ 这篇帖子主要的目的是写一个转录组GSEA分析和可视化通用的函数。起因是我们想要复现一篇文章的GSEA可视化图片,这个Nature文章GSEA可视化挺好的:![](http://image109.360doc.com/DownloadImg/2023/07/1909/269344239_2_2023071909122688_wm.png) (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站上线!
![](http://image109.360doc.com/DownloadImg/2023/07/1909/269344239_3_20230719091226307_wm.png)
![](http://image109.360doc.com/DownloadImg/2023/07/1909/269344239_4_20230719091226448_wm.png) 这里我们以一个单细胞的数据为例,首先做一下差异分析。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的形式直接保存到当前环境路径下! ![](http://image109.360doc.com/DownloadImg/2023/07/1909/269344239_5_20230719091226604_wm.png)
我们对可视化函数进行了设置,NES>0的结果点用红色显示。NES<0的结果用蓝色点显示。这里我们挑选自己感兴趣的通路进行可视化。 #NES>0 p1=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
![](http://image109.360doc.com/DownloadImg/2023/07/1909/269344239_6_20230719091226666_wm.png)
NES<0,这里需要注意一个问题,clusterProfiler和fgsea虽然都是GSEA分析,但是两者得到的结果并不是一摸一样完全相同的,总是有一些出入,这是因为数据库,分析方式不一样。选择需要的包使用一个即可。 #NES<0 p3=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
![](http://image109.360doc.com/DownloadImg/2023/07/1909/269344239_7_20230719091226775_wm.png)
当然了,我们可以利用循环一次性可视化通路。 #批量循环 #NES>0 pathway <- 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)
![](http://image109.360doc.com/DownloadImg/2023/07/1909/269344239_8_20230719091226885_wm.png)
用另一个数据批量做一下下调的。 #NES<0 pathway1 <- 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)
![](http://image109.360doc.com/DownloadImg/2023/07/1909/269344239_9_20230719091227275_wm.png)
这就是所有内容了,希望对你有所启发。这个函数其实并不是很完美,首先是物种只支持小鼠和人,其次是分析参数我没有再设置,用的是我常用的。当然了,这个函数框架我提供了,需要更加个性化的可以自行修改。觉得分享有用的,点个赞呗!
|