不好好作图的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$log2FoldChangenames(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"))

你学会了吗？学会了还不点赞关注！！！