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

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

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

参考文献:

R语言 github man什么意思 r语言gsea代码_开发语言

R语言 github man什么意思 r语言gsea代码_r语言_02

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)

R语言 github man什么意思 r语言gsea代码_R语言 github man什么意思_03

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

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)

R语言 github man什么意思 r语言gsea代码_加载_04

最后图片可以导出为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"))

R语言 github man什么意思 r语言gsea代码_加载_05

你学会了吗?学会了还不点赞关注!!!