本節(jié)概覽：
1.GSEA簡單介紹
2.創(chuàng)建GSEA分析所需的geneList，包含log2FoldChange和ENTREZID信息
3.利用clusterProfiler進(jìn)行GSEA富集GO與KEGG通路
4.GSEA富集結(jié)果可視化：GSEA結(jié)果圖咱扣、 gsearank plot 闹伪、ridgeplot

1. GSEA簡單介紹

以下對GSEA涉及的一些重要概念進(jìn)行了簡單介紹偏瓤，詳細(xì)介紹見：
一文掌握GSEA椰憋，超詳細(xì)教程 - 云+社區(qū) - 騰訊云 (tencent.com)
史上最全GSEA可視化教程橙依，今天讓你徹底搞懂GSEA窗骑！ - 知乎 (zhihu.com)

1.1 GSEA定義與基本原理：

• 定義：
  基因集富集分析(Gene Set Enrichment Analysis, GSEA)是一種計算方法创译，用來確定一組先驗定義的基因集是否在兩種生物狀態(tài)之間顯示出統(tǒng)計學(xué)上顯著的软族、一致的差異。
  官網(wǎng)地址：GSEA (gsea-msigdb.org)

• 基本原理：
  使用預(yù)定義的基因集（通常來自功能注釋或先前實驗的結(jié)果）掖疮，將基因按照在兩類樣本中的差異表達(dá)程度排序氮墨，然后檢驗預(yù)先設(shè)定的基因集合是否在這個排序表的頂端或者底端富集规揪。基因集合富集分析檢測基因集合而不是單個基因的表達(dá)變化字支，因此可以包含這些細(xì)微的表達(dá)變化奸忽，預(yù)期得到更為理想的結(jié)果

• 與GO\KEGG差異基因富集分析區(qū)別：
  差異基因富集分析是先篩選差異基因，再判斷差異基因在哪些注釋的通路存在富集欠雌；這涉及到閾值的設(shè)定富俄，存在一定主觀性并且只能用于表達(dá)變化較大的基因霍比，即我們定義的顯著差異基因悠瞬。而GSEA則不局限于差異基因浅妆，從基因集的富集角度出發(fā)障癌，理論上更容易囊括細(xì)微但協(xié)調(diào)性的變化對生物通路的影響混弥。
  

  gsea.png

1.2 MSigDB（Molecular Signatures Database）:

分子特征數(shù)據(jù)庫。一般進(jìn)行GSEA或GSVA使用的就是該數(shù)據(jù)庫中的基因集蒿涎，我們也可以自定義基因集劳秋。MSigDB所包含的基因集如下所示:

• KEGG信息包含在C2中，GO信息包含在C5中嗽冒。

  
  

1.3 GSEA中關(guān)鍵概念

• ES（Enrichment Score）：富集得分
  ES反應(yīng)基因集成員s在排序列表L的兩端富集的程度剿另。計算方式是贬蛙，從基因集L的第一個基因開始阳准，計算一個累計統(tǒng)計值野蝇。當(dāng)遇到一個落在s里面的基因，則增加統(tǒng)計值乱灵。遇到一個不在s里面的基因痛倚，則降低統(tǒng)計值澜躺。 每一步統(tǒng)計值增加或減少的幅度與基因的表達(dá)變化程度（fold-change值）是相關(guān)的掘鄙。富集得分ES最后定義為最大的峰值操漠。正值ES表示基因集在列表的頂部富集，負(fù)值ES表示基因集在列表的底部富集撞秋。
  p-value用來評估富集得分(ES)的顯著性吻贿，通過排列檢驗 (permutation test)計算觀察到的富集得分(ES)出現(xiàn)的可能性舅列。

• NES (Normalized Enrichment Score)：標(biāo)準(zhǔn)化富集得分
  每個基因子集s計算得到的ES根據(jù)基因集的大小進(jìn)行標(biāo)準(zhǔn)化得到標(biāo)準(zhǔn)化富集得分Normalized Enrichment Score (NES)。隨后會針對NES計算假陽性率FDR把敞。

• Leading-edge subset：領(lǐng)頭基因亞集
  對富集貢獻(xiàn)最大的基因成員

• 一般認(rèn)為|NES|>1先巴，p-value<0.05伸蚯，F(xiàn)DR<0.25的通路是顯著富集的剂邮。
  |NES|值越大挥萌，F(xiàn)DR值就越小枉侧，說明分析的結(jié)果可信度越高榨馁。

2. 創(chuàng)建GSEA分析所需的geneList

在了解了GSEA基本概念后就可以正式開始實操了翼虫，首先需要將基因按照在兩類樣本中的差異表達(dá)程度排序珍剑。
下面我們構(gòu)建包含了geneList，里面含有從大到小排序的log2FoldChange和對應(yīng)的ENTREZID信息：

rm(list = ls())  
options(stringsAsFactors = F)
# library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(clusterProfiler)
library(enrichplot)
library(tidyverse)
library(ggstatsplot)

setwd("C:/Users/Lenovo/Desktop/test")
load(list.files(path = "./3.DEG",pattern = 'DEG_results.Rdata',full.names = T))
dir.create("5.GSEA_kegg_go")
setwd("5.GSEA_kegg_go")

## 物種設(shè)置
organism = 'mmu'    #  人類'hsa' 小鼠'mmu'   
OrgDb = 'org.Mm.eg.db'#人類"org.Hs.eg.db" 小鼠"org.Mm.eg.db"

#### 按照需要可選擇不同的DEG方法數(shù)據(jù)集 ####
need_DEG <- DEG_DESeq2
need_DEG <- need_DEG[,c(2,5)] #選擇log2FoldChange和pvalue（湊成數(shù)據(jù)框）

colnames(need_DEG) <- c('log2FoldChange','pvalue')
need_DEG$SYMBOL <- rownames(need_DEG)

##### 創(chuàng)建gsea分析的geneList（包含從大到小排列的log2FoldChange和ENTREZID信息）####
#轉(zhuǎn)化id  
df <- bitr(rownames(need_DEG), 
           fromType = "SYMBOL",
           toType =  "ENTREZID",
           OrgDb = OrgDb) #人數(shù)據(jù)庫org.Hs.eg.db 小鼠org.Mm.eg.db
need_DEG <- merge(need_DEG, df, by='SYMBOL')  #按照SYMBOL合并注釋信息
geneList <- need_DEG$log2FoldChange
names(geneList) <- need_DEG$ENTREZID
geneList <- sort(geneList, decreasing = T)   #從大到小排序

3. 利用clusterProfiler包進(jìn)行GSEA富集

clusterProfiler包內(nèi)的gseGO()和gseKEGG()函數(shù)可以很方便地對GO與KEGG通路進(jìn)行GSEA， 再使用DOSE::setReadable轉(zhuǎn)化id 饰序。

##### gsea富集 ####
KEGG_kk_entrez <- gseKEGG(geneList     = geneList,
                   organism     = organism, #人hsa 鼠mmu
                   pvalueCutoff = 0.25)  #實際為padj閾值可調(diào)整 
KEGG_kk <- DOSE::setReadable(KEGG_kk_entrez, 
                             OrgDb=OrgDb,
                             keyType='ENTREZID')#轉(zhuǎn)化id             
  
GO_kk_entrez <- gseGO(geneList     = geneList,
               ont          = "ALL",  # "BP"、"MF"和"CC"或"ALL"
               OrgDb        = OrgDb,#人類org.Hs.eg.db 鼠org.Mm.eg.db
               keyType      = "ENTREZID",
               pvalueCutoff = 0.25)   #實際為padj閾值可調(diào)整
GO_kk <- DOSE::setReadable(GO_kk_entrez, 
                           OrgDb=OrgDb,
                           keyType='ENTREZID')#轉(zhuǎn)化id 

save(KEGG_kk_entrez, GO_kk_entrez, file = "GSEA_result.RData")

4. GSEA富集結(jié)果可視化

GSEA的可視化主要是GSEA結(jié)果圖由缆、 gsearank plot和ridgeplot山脊圖。
同樣也可以進(jìn)行其他可視化如barplot是晨、dotplot罩缴、cnetplot等等层扶，詳見RNA-seq入門的簡單實戰(zhàn)（六）：GO镜会、KEGG富集分析與超全可視化攻略
或者參閱說明書Chapter 15 Visualization of functional enrichment result | Biomedical Knowledge Mining using GOSemSim and clusterProfiler (yulab-smu.top)戳表，這里就不再進(jìn)行展示啦

4.1 gseaplot GSEA結(jié)果圖

下面選取KEGG通路的富集結(jié)果進(jìn)行g(shù)seaplot繪圖示范

• 首先對富集結(jié)果進(jìn)行條件篩選匾旭，一般認(rèn)為|NES|>1，NOM pvalue<0.05女蜈，F(xiàn)DR（padj）<0.25的通路是顯著富集的鞭光；還可以從結(jié)果中細(xì)分出上下調(diào)通路單獨繪圖惰许，以下代碼僅展示KEGG通路富集結(jié)果的上調(diào)通路汹买。
• gseaplot2(）函數(shù)既可以對單獨的通路繪圖聊倔，也可以合并幾個通路一起繪圖耙蔑；各類詳細(xì)參數(shù)設(shè)置見以下代碼處

##選取富集結(jié)果
kk_gse <- KEGG_kk
kk_gse_entrez <- KEGG_kk_entrez

###條件篩選 
#一般認(rèn)為|NES|>1甸陌，NOM pvalue<0.05盐股，F(xiàn)DR（padj）<0.25的通路是顯著富集的
kk_gse_cut <- kk_gse[kk_gse$pvalue<0.05 & kk_gse$p.adjust<0.25 & abs(kk_gse$NES)>1]
kk_gse_cut_down <- kk_gse_cut[kk_gse_cut$NES < 0,]
kk_gse_cut_up <- kk_gse_cut[kk_gse_cut$NES > 0,]

#選擇展現(xiàn)NES前幾個通路 
down_gsea <- kk_gse_cut_down[tail(order(kk_gse_cut_down$NES,decreasing = T),10),]
up_gsea <- kk_gse_cut_up[head(order(kk_gse_cut_up$NES,decreasing = T),10),]
diff_gsea <- kk_gse_cut[head(order(abs(kk_gse_cut$NES),decreasing = T),10),]


#### 經(jīng)典的GSEA圖 
up_gsea$Description
i=2
gseap1 <- gseaplot2(kk_gse,
                    up_gsea$ID[i],#富集的ID編號
                    title = up_gsea$Description[i],#標(biāo)題
                    color = "red", #GSEA線條顏色
                    base_size = 20,#基礎(chǔ)字體大小
                    rel_heights = c(1.5, 0.5, 1),#副圖的相對高度
                    subplots = 1:3,   #要顯示哪些副圖 如subplots=c(1,3) #只要第一和第三個圖
                    ES_geom = "line", #enrichment score用線還是用點"dot"
                    pvalue_table = T) #顯示pvalue等信息
ggsave(gseap1, filename = 'GSEA_up_1.pdf', width =10, height =8)
 
#### 合并 GSEA通路 
gseap2 <- gseaplot2(kk_gse,
                    up_gsea$ID,#富集的ID編號
                    title = "UP_GSEA_all",#標(biāo)題
                    color = "red",#GSEA線條顏色
                    base_size = 20,#基礎(chǔ)字體大小
                    rel_heights = c(1.5, 0.5, 1),#副圖的相對高度
                    subplots = 1:3, #要顯示哪些副圖 如subplots=c(1,3) #只要第一和第三個圖
                    ES_geom = "line",#enrichment score用線還是用點"dot"
                    pvalue_table = T) #顯示pvalue等信息
ggsave(gseap2, filename = "GSEA_up_all.pdf",width =12,height =12)

下面解釋一下GSEA圖的含義：

• 第1部分是ES折線圖卵酪，離垂直距離x=0軸最遠(yuǎn)的峰值便是基因集的ES值溃卡，峰出現(xiàn)在排序基因集的前端（ES值大于0）則說明通路上調(diào)，出現(xiàn)在后端（ES值小于0）則說明通路下調(diào)沫换。
• 第二部分為基因集成員位置圖讯赏，用豎線標(biāo)記了基因集中各成員出現(xiàn)在基因排序列表中的位置漱挎。若豎線集中分布在基因排序列表的前端或后端雀哨，說明該基因集通路上調(diào)或下調(diào)雾棺；若豎線較均勻分布在基因排序列表中捌浩，則說明該基因集通路在比較的兩個數(shù)據(jù)中無明顯變化。
  紅色部分對應(yīng)的基因在實驗組中高表達(dá)进统，藍(lán)色部分對應(yīng)的基因在對照組中高表達(dá)螟碎，
  leading edge subset 是（0,0）到曲線峰值ES出現(xiàn)對應(yīng)的這部分基因成員掉分。
• 第三部分是排序后所有基因rank值（由log2FoldChang值計算得出）的分布，以灰色面積圖顯展示尔崔。

4.2 gsearank plot 繪制特定基因集的基因排序列表

gsearank()展示特定基因集的排序，橫坐標(biāo)為基因排序消返，縱坐標(biāo)為ES值耘拇，利用cowplot和ggplot2包可以批量出圖惫叛。

## gsearank plot 繪制出屬于特定基因集的基因排序列表
##繪制up_gsea前3個富集通路
library(cowplot)
library(ggplot2)
pp <- lapply(1:3, function(i) {
  anno <- up_gsea[i, c("NES", "pvalue", "p.adjust")]
  lab <- paste0(names(anno), "=",  round(anno, 3), collapse="\n")
  
  gsearank(kk_gse,
           up_gsea$ID[1], 
           title = up_gsea$Description[i]) + 
    xlab(NULL) +
    # ylab(NULL) +
    annotate("text", 10000,
             up_gsea[i, "enrichmentScore"] * .75, 
             label = lab, 
             hjust=0, vjust=0)
})
rankp <- plot_grid(plotlist=pp, ncol=1)
ggsave(rankp, filename = "gsearank_up.pdf",width=8,height=10)

4.3 ridgeplot山脊圖

展示富集通路的核心富集基因的表達(dá)分布妻熊，x軸為富集通路的核心富集基因表達(dá)變化的log2倍仑最，值為正值表示表達(dá)上調(diào)警医，值為負(fù)值表示表達(dá)下調(diào)。

## ridgeplot
ridgep <- ridgeplot(kk_gse_entrez,
                    showCategory = 15,
                    fill = "p.adjust",
                    core_enrichment = TRUE,
                    label_format = 30, #設(shè)置軸標(biāo)簽文字的每行字符數(shù)長度侈玄，過長則會自動換行序仙。
                    orderBy = "NES",
                    decreasing = F) 
ggsave(ridgep,filename = 'ridgeplot.pdf',width =10,height =10)

（之前運行報錯解決方法見ridgeplot報錯：Error in ans[ypos] <- rep(yes, length.out = len)[ypos] : replacement has ... ）

4.4 其他富集結(jié)果可視化圖

dotplot cnetplot emapplot treeplot heatplot upsetplot
詳見RNA-seq入門的簡單實戰(zhàn)（六）：GO诱桂、KEGG富集分析與超全可視化攻略

GSEA分析和可視化到這就結(jié)束啦挥等，下一節(jié)介紹GSVA的使用

參考資料
?? Introduction | Biomedical Knowledge Mining using GOSemSim and clusterProfiler (yulab-smu.top)
一文掌握GSEA肝劲，超詳細(xì)教程 - 云+社區(qū) - 騰訊云 (tencent.com)
史上最全GSEA可視化教程辞槐，今天讓你徹底搞懂GSEA！ - 知乎 (zhihu.com)
GitHub - jmzeng1314/GEO
【生信技能樹】轉(zhuǎn)錄組測序數(shù)據(jù)分析_嗶哩嗶哩_bilibili
【生信技能樹】GEO數(shù)據(jù)庫挖掘_嗶哩嗶哩_bilibili

RNA-seq實戰(zhàn)系列文章：
RNA-seq入門實戰(zhàn)（零）：RNA-seq流程前的準(zhǔn)備——Linux與R的環(huán)境創(chuàng)建
RNA-seq入門實戰(zhàn)（一）：上游數(shù)據(jù)下載卜范、格式轉(zhuǎn)化和質(zhì)控清洗
RNA-seq入門實戰(zhàn)（二）：上游數(shù)據(jù)的比對計數(shù)——Hisat2+ featureCounts 與 Salmon
RNA-seq入門實戰(zhàn)（三）：從featureCounts與Salmon輸出文件獲取counts矩陣
RNA-seq入門實戰(zhàn)（四）：差異分析前的準(zhǔn)備——數(shù)據(jù)檢查
RNA-seq入門實戰(zhàn)（五）：差異分析——DESeq2 edgeR limma的使用與比較
RNA-seq入門實戰(zhàn)（六）：GO海雪、KEGG富集分析與enrichplot超全可視化攻略
RNA-seq入門實戰(zhàn)（七）：GSEA——基因集富集分析
RNA-seq入門實戰(zhàn)（八）：GSVA——基因集變異分析
RNA-seq入門實戰(zhàn)（九）：PPI蛋白互作網(wǎng)絡(luò)構(gòu)建（上）——STRING數(shù)據(jù)庫的使用
RNA-seq入門實戰(zhàn)（十）：PPI蛋白互作網(wǎng)絡(luò)構(gòu)建（下）——Cytoscape軟件的使用
RNA-seq入門實戰(zhàn)（十一）：WGCNA加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析——關(guān)聯(lián)基因模塊與表型