之前我們在講轉(zhuǎn)錄組系列的時候，說過差異基因的功能富集，用的是GO和KEGG分析。但是這遠(yuǎn)遠(yuǎn)不夠，很多研究者更喜歡使用GSEA，全名是Gene Set Enrichment Analysis (基因集富集分析)。GSEA在一定程度上與GO一樣，但是兩者具有巨大的差別。GO使用的是差異基因，因為閾值的設(shè)定是人為的，所以很有可能遺漏一些重要基因，僅僅是因為這些基因的變化較小。而GSEA則不同，它需要的是對所有的基因進(jìn)行分析，因此能夠保留更多的信息。

通俗的說，GSEA的適用場景是：在兩種不同的生物學(xué)狀態(tài)下，可以理解為處理組與對照組，判斷某一組基因集其表達(dá)模式更接近于那種過程或者通路，從而推斷這些基因?qū)@個生物學(xué)過程的重要貢獻(xiàn)。在分析結(jié)果上GSEA也是更加可靠，可信度較高的。

富集分析是轉(zhuǎn)錄組必做的內(nèi)容，這里我們通過對一篇Cell文章GSEA的可視化來講講GSEA的分析。GSEA通常以軟件客戶端分析，這樣對于編程的要求低了，但是缺點是導(dǎo)出的圖片清晰度不夠。所以這里我們介紹一下R語言的版本！

參考文獻(xiàn)：

圖片

圖片

GSEA分析需要差異基因和變化倍數(shù)，加載差異分析數(shù)據(jù)：

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

首先定義基因集，但是我們之前講過msigdbr這個包（有了這個包，豬的GSEA和GSVA分析也不在話下（第一集）），就不需要從GSEA官網(wǎng)下載基因集了，只需要加載包即可。這里以GO為例。

## 預(yù)定義基因集
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)#將_轉(zhuǎn)化為空格

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

對基因排序：

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

使用fgsea函數(shù)進(jìn)行分析：

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

篩選上下調(diào)最顯著的作圖，可以一次性呈現(xiàn)多個通路：

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)

圖片

當(dāng)然也可以自定義通路可視化：

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)

圖片

最后圖片可以導(dǎo)出為pdf，在AI中修飾，實現(xiàn)與Cell文章中一樣的效果。此外，在可視化上，也可以像GSEA軟件導(dǎo)出的圖片一樣，導(dǎo)出單個通路的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"))

圖片

你學(xué)會了嗎？學(xué)會了還不點贊關(guān)注！?。?br> 更多精彩請訪問：KS科研分享與服務(wù)