不好好作图的NCS系列(五)：从这篇Cell学习GSEA的R语言

之前我们在讲转录组系列的时候，说过差异基因的功能富集，用的是GO和KEGG分析。但是这远远不够，很多研究者更喜欢使用GSEA，全名是Gene Set Enrichment Analysis (基因集富集分析)。GSEA在一定程度上与GO一样，但是两者具有巨大的差别。GO使用的是差异基因，因为阈值的设定是人为的，所以很有可能遗漏一些重要基因，仅仅是因为这些基因的变化较小。而GSEA则不同，它需要的是对所有的基因进行分析，因此能够保留更多的信息。

通俗的说，GSEA的适用场景是：在两种不同的生物学状态下，可以理解为处理组与对照组，判断某一组基因集其表达模式更接近于那种过程或者通路，从而推断这些基因对这个生物学过程的重要贡献。在分析结果上GSEA也是更加可靠，可信度较高的。

富集分析是转录组必做的内容，这里我们通过对一篇Cell文章GSEA的可视化来讲讲GSEA的分析。GSEA通常以软件客户端分析，这样对于编程的要求低了，但是缺点是导出的图片清晰度不够。所以这里我们介绍一下R语言的版本！

参考文献：

图片 图片

GSEA分析需要差异基因和变化倍数，加载差异分析数据：

library(msigdbr)
library(fgsea)
library(dplyr)
library(ggplot2)
library(Seurat)
library(tidyverse)
setwd("F:/生物信息学")
A <- read.csv("DEGs_trans.csv",header = T)

首先定义基因集，但是我们之前讲过msigdbr这个包（有了这个包，猪的GSEA和GSVA分析也不在话下（第一集）），就不需要从GSEA官网下载基因集了，只需要加载包即可。这里以GO为例。

## 预定义基因集
geneset_GO = msigdbr(species = "Homo sapiens",
                     category = "C5", 
                     subcategory = "GO:BP") %>% dplyr::select(gs_name,gene_symbol)
#geneset_GO$gs_name <- gsub('GOBP_','',geneset_GO$gs_name)#去除前缀KEGG_
#geneset_GO$gs_name <- tolower(geneset_GO$gs_name)#将大写换为小写
#geneset_GO$gs_name <- gsub('_',' ',geneset_GO$gs_name)#将_转化为空格

gsea_geneset_GO <- geneset_GO %>% split(x = .$gene_symbol, f = .$gs_name)

对基因排序：

## 根据logfc降序排列基因
Agsea <- A[order(A$log2FoldChange,decreasing = T),]
ranks <- Agsea$log2FoldChange
names(ranks) <- Agsea$gene

使用fgsea函数进行分析：

Agsea_res <- fgsea(pathways = gsea_geneset_GO, 
                           stats = ranks,
                           minSize=5,
                           maxSize=500,
                           nperm=1000)

筛选上下调最显著的作图，可以一次性呈现多个通路：

Up <- Agsea_res[ES > 0][head(order(pval), n=10), pathway]
Down <- Agsea_res[ES < 0][head(order(pval), n=10), pathway]
topPathways <- c(Up, rev(Down))
plotGseaTable(gsea_geneset_GO[topPathways], ranks, Agsea_res, 
              gseaParam = 0.5)

图片

当然也可以自定义通路可视化：

GO_secleted <- c("GOBP_PROTEIN_LOCALIZATION_TO_KINETOCHORE",
                 "GOBP_PROTEIN_LOCALIZATION_TO_CHROMOSOME",
                 "GOBP_CHRONIC_INFLAMMATORY_RESPONSE",
                 "GOBP_NECROPTOTIC_PROCESS",
                 "GOBP_ACTIN_FILAMENT_ORGANIZATION",
                 "GOBP_DNA_DEALKYLATION",
                 "GOBP_REGULATION_OF_GERMINAL_CENTER_FORMATION",
                 "GOBP_MEMBRANE_REPOLARIZATION",
                 "GOBP_VIRION_ATTACHMENT_TO_HOST_CELL",
                 "GOBP_NEGATIVE_REGULATION_OF_ACTIN_NUCLEATION",
                 "GOBP_REGULATION_OF_PROTEIN_KINASE_A_SIGNALING",
                 "GOBP_ORGANELLE_FISSION",
                 "GOBP_MUSCLE_SYSTEM_PROCESS",
                 "GOBP_REGULATION_OF_CELL_DEVELOPMENT",
                 "GOBP_POSITIVE_REGULATION_OF_NEUROGENESIS",
                 "GOBP_DENDRITIC_CELL_APOPTOTIC_PROCESS")
plotGseaTable(gsea_geneset_GO[GO_secleted], ranks, Agsea_res, 
              gseaParam = 0.5)

图片

最后图片可以导出为pdf，在AI中修饰，实现与Cell文章中一样的效果。此外，在可视化上，也可以像GSEA软件导出的图片一样，导出单个通路的NES图。

library(ggplot2)
plotEnrichment(gsea_geneset_GO[["GOBP_DNA_REPAIR"]],
               ranks,ticksSize = 0.5) + labs(title="GOBP_DNA_REPAIR")+
  theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"))

图片

你学会了吗？学会了还不点赞关注！！！
更多精彩请访问：KS科研分享与服务