Introduction
Gene Set Enrichment Analysis (GSEA)
是一種用于分析基因表達(dá)數(shù)據(jù)的計(jì)算生物學(xué)方法，旨在揭示與特定生物學(xué)過(guò)程、通路或功能相關(guān)的基因表達(dá)模式。GSEA
最初由麻省理工學(xué)院的研究人員開(kāi)發(fā) (1)，它不同于傳統(tǒng)的基因差異分析方法，如
t 檢驗(yàn)或ANOVA，這些方法通常關(guān)注單個(gè)基因的表達(dá)差異。 相反，GSEA
專注于整個(gè)基因集的表達(dá)趨勢(shì)，這些基因集通常與特定的生物學(xué)過(guò)程或通路相關(guān)聯(lián)。
在典型的實(shí)驗(yàn)中，從屬于兩類之一的樣本集合中生成數(shù)千個(gè)基因的 mRNA
表達(dá)譜，例如，對(duì)藥物敏感和耐藥的腫瘤。
基因可以根據(jù)它們?cè)陬悇e之間的差異表達(dá)在排序列表L中排序。挑戰(zhàn)在于從這個(gè)列表中提取含義。
一種常見(jiàn)的方法是關(guān)注 L
頂部和底部的少數(shù)基因（即那些表現(xiàn)出最大差異的基因），以辨別明顯的生物學(xué)線索。這種方法有一些主要的局限性:
在校正多個(gè)假設(shè)檢驗(yàn)后，沒(méi)有任何單個(gè)基因可以滿足統(tǒng)計(jì)顯著性的閾值，因?yàn)橄鄬?duì)于微陣列技術(shù)固有的噪聲，相關(guān)的生物學(xué)差異是適度（不是很明顯）的。
或者，人們可能會(huì)留下一長(zhǎng)串具有統(tǒng)計(jì)學(xué)意義的基因（太多了），而沒(méi)有任何統(tǒng)一的生物學(xué)主題。解釋可能是令人畏懼的和臨時(shí)的，取決于生物學(xué)家的專業(yè)領(lǐng)域。
單基因分析可能會(huì)錯(cuò)過(guò)對(duì)通路的重要影響。細(xì)胞過(guò)程常常影響協(xié)同作用的基因組。編碼代謝途徑成員的所有基因增加
20% 可能會(huì)顯著改變通過(guò)該途徑的通量，并且可能比單個(gè)基因增加 20
倍更重要。當(dāng)不同的群體研究相同的生物系統(tǒng)時(shí)，兩項(xiàng)研究中具有統(tǒng)計(jì)學(xué)意義的基因列表可能顯示出令人沮喪的很少的重疊（但是通路可能有重疊）。
以下是 GSEA 的簡(jiǎn)要介紹和工作原理：
基因集定義：首先，GSEA
需要一個(gè)事先定義好的基因集合，這些基因通常按照其在特定生物學(xué)通路、功能類別或疾病過(guò)程中的作用進(jìn)行組織。這些基因集可以來(lái)自公共數(shù)據(jù)庫(kù)，如Gene
Ontology、KEGG
Pathway、Reactome，或者是研究者自己根據(jù)文獻(xiàn)和領(lǐng)域知識(shí)構(gòu)建的。基因表達(dá)數(shù)據(jù)：GSEA
需要分析的基因表達(dá)數(shù)據(jù)，通常是從微陣列或RNA測(cè)序?qū)嶒?yàn)中獲得的。這些數(shù)據(jù)包括不同條件或樣本中基因的表達(dá)水平。排列(permutation)檢驗(yàn)：GSEA
的核心思想是通過(guò)對(duì)基因集的成員在整個(gè)基因表達(dá)數(shù)據(jù)中的排列來(lái)確定它們的富集程度。具體來(lái)說(shuō)，GSEA
將所有基因根據(jù)其在不同條件下的表達(dá)水平進(jìn)行排序。然后，它從基因集的一端開(kāi)始，計(jì)算該基因集中的基因在排序列表中的偏離程度。如果基因集中的基因在排序列表的某個(gè)位置中集中排列，說(shuō)明該基因集富集在特定的生物學(xué)過(guò)程或通路中。這個(gè)過(guò)程通過(guò)構(gòu)建一個(gè)富集分?jǐn)?shù)（Enrichment
Score）來(lái)量化。統(tǒng)計(jì)顯著性：對(duì)于每個(gè)基因集，GSEA
計(jì)算一個(gè)富集分?jǐn)?shù)，并基于隨機(jī)排列檢驗(yàn)來(lái)估計(jì)其統(tǒng)計(jì)顯著性。如果一個(gè)基因集的富集分?jǐn)?shù)在隨機(jī)排列中的分布表現(xiàn)出顯著差異，那么就認(rèn)為這個(gè)基因集在樣本中富集。結(jié)果可視化：最后，GSEA
會(huì)生成結(jié)果報(bào)告，其中包括富集分?jǐn)?shù)、基因集的統(tǒng)計(jì)顯著性以及相關(guān)通路或功能的信息。這些結(jié)果可視化為富集圖譜，用于展示不同基因集的富集情況。
總之，GSEA
是一種用于揭示基因表達(dá)數(shù)據(jù)中生物學(xué)意義的強(qiáng)大工具，它可以幫助研究人員理解不同生物學(xué)過(guò)程、通路或功能在不同條件下的活動(dòng)變化。它已經(jīng)廣泛應(yīng)用于生物醫(yī)學(xué)研究中，特別是在基因表達(dá)分析和疾病機(jī)制研究中發(fā)揮了重要作用。
Alogrithm
? Step 1: 計(jì)算富集分?jǐn)?shù)（Enrichment
Score，ES）。計(jì)算一個(gè)富集分?jǐn)?shù)（ES），反映了一個(gè)基因集合S在整個(gè)排序列表L的極端位置（頂部或底部）上過(guò)度表示的程度。分?jǐn)?shù)的計(jì)算是通過(guò)遍歷列表L來(lái)實(shí)現(xiàn)的，在遇到基因集S中的基因時(shí)，增加一個(gè)運(yùn)行總和統(tǒng)計(jì)量；在遇到不在基因集S中的基因時(shí)，減少這個(gè)統(tǒng)計(jì)量。增量的大小取決于基因與表型的相關(guān)性。富集分?jǐn)?shù)是在隨機(jī)漫步中遇到的與零的最大偏差；它對(duì)應(yīng)于一種加權(quán)的Kolmogorov–Smirnov-like統(tǒng)計(jì)量。
? Step 2:
估計(jì)ES的顯著性水平。通過(guò)使用經(jīng)驗(yàn)性的基于表型的置換測(cè)試過(guò)程來(lái)估計(jì)ES的統(tǒng)計(jì)顯著性（名義P值），該過(guò)程保留了基因表達(dá)數(shù)據(jù)的復(fù)雜相關(guān)結(jié)構(gòu)。具體來(lái)說(shuō)，對(duì)表型標(biāo)簽進(jìn)行排列并重新計(jì)算基因集在排列后數(shù)據(jù)中的ES，這樣生成了ES的零分布。然后，計(jì)算觀察到的ES的經(jīng)驗(yàn)性名義P值相對(duì)于這個(gè)零分布。重要的是，類別標(biāo)簽的排列保留了基因與基因之間的相關(guān)性，因此提供了一個(gè)更符合生物學(xué)的顯著性評(píng)估，而不是通過(guò)排列基因而獲得的評(píng)估。
實(shí)際上GSEA好像提供了兩種排列方法 (2)，另一種是基因排列，直接將觀察到的路徑
ES 與通過(guò)使用匹配大小（例如 1,000
次）的隨機(jī)采樣基因集重復(fù)分析而獲得的分?jǐn)?shù)分布進(jìn)行比較，表型排列應(yīng)與大量重復(fù)一起使用（例如，每個(gè)條件至少十次）。
與基因集排列方法相比，表型排列方法的主要優(yōu)點(diǎn)是，它在排列過(guò)程中保持了具有生物學(xué)重要基因相關(guān)性的基因集結(jié)構(gòu)。
表型排列的計(jì)算成本很高，并且對(duì)于當(dāng)前版本的
GSEA，需要自定義編程來(lái)分別計(jì)算數(shù)千個(gè)表型隨機(jī)化的 ES 和差異表達(dá)統(tǒng)計(jì)數(shù)據(jù)。
? Step 3:
多重假設(shè)檢驗(yàn)的調(diào)整。當(dāng)評(píng)估整個(gè)基因集合數(shù)據(jù)庫(kù)時(shí)，會(huì)調(diào)整估計(jì)的顯著性水平以考慮多重假設(shè)檢驗(yàn)。首先，對(duì)每個(gè)基因集的ES進(jìn)行歸一化，以考慮集合的大小，得到一個(gè)歸一化富集分?jǐn)?shù)（NES）。然后，通過(guò)計(jì)算與每個(gè)NES相對(duì)應(yīng)的虛假發(fā)現(xiàn)率（FDR）來(lái)控制假陽(yáng)性比例。FDR是給定NES的一個(gè)集合代表虛假陽(yáng)性發(fā)現(xiàn)的估計(jì)概率；它是通過(guò)比較觀察到的NES和零分布的尾部來(lái)計(jì)算的。
Leading-Edge Subset
基因集可以通過(guò)使用多種方法來(lái)定義，但并非基因集的所有成員通常都會(huì)參與生物過(guò)程。通常，提取有助于
ES 的高分基因集的核心成員是有用的。
將Leading-Edge Subset (前沿子集)定義為基因集 S 中出現(xiàn)在排序列表 L
中或之前運(yùn)行總和達(dá)到與零的最大偏差的那些基因（圖
1B）。前沿子集可以解釋為解釋富集信號(hào)的基因集的核心。
MSigDB
他們?cè)诎l(fā)表GSEA時(shí)同時(shí)創(chuàng)建了MSigDB
(Molecular Signatures Database)，分子特征數(shù)據(jù)庫(kù) (MSigDB)
是一個(gè)包含數(shù)以萬(wàn)計(jì)的注釋基因集的資源，可與 GSEA
軟件一起使用，分為人類和小鼠集合。
比如說(shuō)人類分子特征數(shù)據(jù)庫(kù) (MSigDB) 中的 33591 個(gè)基因集分為 9
個(gè)主要集合和幾個(gè)子集合，可以自己去下載研究相關(guān)的重要基因集：
| Collection | Description |
|---|---|
| H: hallmark gene sets | Hallmark gene sets summarize and represent specific well-defined biological states or processes and display coherent expression. These gene sets were generated by a computational methodology based on identifying overlaps between gene sets in other MSigDB collections and retaining genes that display coordinate expression. details |
| C1: positional gene sets | Gene sets corresponding to human chromosome cytogenetic bands. details |
| C2: curated gene sets | Gene sets in this collection are curated from various sources, including online pathway databases and the biomedical literature. Many sets are also contributed by individual domain experts. The gene set page for each gene set lists its source. The C2 collection is divided into the following two subcollections: Chemical and genetic perturbations (CGP) and Canonical pathways (CP). details |
| C3: regulatory target gene sets | Gene sets representing potential targets of regulation by transcription factors or microRNAs. The sets consist of genes grouped by elements they share in their non-protein coding regions. The elements represent known or likely cis-regulatory elements in promoters and 3’-UTRs. The C3 collection is divided into two subcollections: microRNA targets (MIR) and transcription factor targets (TFT). details |
| C4: computational gene sets | Computational gene sets defined by mining large collections of cancer-oriented microarray data. The C4 collection is divided into two subcollections: CGN and CM. details |
| C5: ontology gene sets | Gene sets that contain genes annotated by the same ontology term. The C5 collection is divided into two subcollections, the first derived from the Gene Ontology resource (GO) which contains BP, CC, and MF components and a second derived from the Human Phenotype Ontology (HPO). details |
| C6: oncogenic signature gene sets | Gene sets that represent signatures of cellular pathways which are often dis-regulated in cancer. The majority of signatures were generated directly from microarray data from NCBI GEO or from internal unpublished profiling experiments involving perturbation of known cancer genes. details |
| C7: immunologic signature gene sets | Gene sets that represent cell states and perturbations within the immune system. details |
| C8: cell type signature gene sets | Gene sets that contain curated cluster markers for cell types identified in single-cell sequencing studies of human tissue. details |
Table 1: Gene sets in the Human Molecular Signatures Database .
Example
這里介紹在R中完成GSEA的方法
- 對(duì)GSEA分析的geneList排序
明確我們用來(lái)排序的指標(biāo)（metric）。
目前大部分分析是只比較兩組間差異的，所以會(huì)將基因按照在兩類樣本中的差異表達(dá)程度（一般是log2FoldChange）排序，但實(shí)際上這里可以使用多種指標(biāo)比如相關(guān)性，p值，也可以直接進(jìn)行多組的分析比較（當(dāng)然一定要明確排序兩端代表的生物學(xué)意義）。盡量不要進(jìn)行原始基因表格的篩選，比如只取DEGs，因?yàn)镚SEA本身就是一個(gè)無(wú)閾值的富集方法，相比f(wàn)isher.test等閾值方法考慮的會(huì)更多一些。
#直接使用DOSE提供的一個(gè)geneList，已經(jīng)通過(guò)log2FoldChange從大到小排序過(guò)了，而且這個(gè)向量的name是每一個(gè)entrez gene id, value是log2FoldChange值。
library(DOSE)
data(geneList, package = "DOSE")
head(geneList)
## 4312 8318 10874 55143 55388 991
## 4.572613 4.514594 4.418218 4.144075 3.876258 3.677857
- 使用R包進(jìn)行GSEA富集
如果要使用GSEA官網(wǎng)提供的帶GUI軟件，可以參考這篇nature
protocol，里面詳細(xì)介紹了每一步的操作 (2)。
這里我們使用更方便一點(diǎn)的R包來(lái)進(jìn)行。clusterProfiler包內(nèi)的gseGO()和gseKEGG()函數(shù)可以很方便地對(duì)GO與KEGG通路進(jìn)行GSEA。
如果明確我們要富集的基因集的話，也可以自己構(gòu)建并使用clusterProfiler::GSEA()富集。
library(clusterProfiler)
KEGG_kk_entrez <- gseKEGG(geneList = geneList,
organism = "hsa",pvalueCutoff = 0.25) #實(shí)際為padj閾值,可調(diào)整。
但上面這行代碼在我的機(jī)器上結(jié)果是– No gene can be
mapped….，可能是我機(jī)器上的clusterProfiler包內(nèi)部數(shù)據(jù)庫(kù)問(wèn)題。
但是沒(méi)有關(guān)系，我一般也會(huì)自己從kegg上把想分析的通路信息都爬下來(lái)，操作見(jiàn)上一篇關(guān)于kegg
api的介紹，
具體的代碼都集成在ReporterScore里面了，下面直接使用ReporterScore進(jìn)行GSEA分析。
library(ReporterScore)
#加載human kegg數(shù)據(jù)庫(kù)
load_org_pathway(org = "hsa")
#這里包含了human所有g(shù)ene-kegg_pathway的對(duì)應(yīng)關(guān)系，可以使用download_org_pathway()進(jìn)行更新。
head(hsa_kegg_pathway$all_org_gene)
## pathway_id kegg_gene_id gene_symbol gene_desc KO_id
## 1 hsa00010 3101 HK3 hexokinase 3 K00844
## 2 hsa00010 3098 HK1 hexokinase 1 K00844
## 3 hsa00010 3099 HK2 hexokinase 2 K00844
## 4 hsa00010 80201 HKDC1 hexokinase domain containing 1 K00844
## 5 hsa00010 2645 GCK glucokinase K12407
## 6 hsa00010 2821 GPI glucose-6-phosphate isomerase K01810
# 假設(shè)我們的gene table如下，行名是gene symbol，使用t.test比較WT-OE兩組差異
ko.test(genedf,"Group",metadata,method = "t.test")->da_res
# 統(tǒng)計(jì)結(jié)果如下
head(da_res)
## KO_id average_WT sd_WT average_OE sd_OE
## MME MME 0.0026728975 0.0011094132 0.0020694937 0.0014922004
## PEX11A PEX11A 0.0004554658 0.0001951678 0.0024480601 0.0014916566
## DSC2 DSC2 0.0027767713 0.0006253559 0.0025519275 0.0004966801
## PLIN5 PLIN5 0.0005779169 0.0008952163 0.0005197504 0.0001435245
## MIRLET7A2 MIRLET7A2 0.0020807307 0.0007731661 0.0014321838 0.0004273716
## NCBP2L NCBP2L 0.0021064422 0.0005243558 0.0017419317 0.0005382770
## diff_mean Highest p.value p.adjust
## MME -0.00060340387 WT 0.2200529827 0.321714887
## PEX11A 0.00199259427 OE 0.0001376795 0.001449258
## DSC2 -0.00022484380 WT 0.2852910552 0.392962886
## PLIN5 -0.00005816647 WT 0.8072016928 0.865027977
## MIRLET7A2 -0.00064854693 WT 0.0095033556 0.028624565
## NCBP2L -0.00036451047 WT 0.0707452991 0.134260284
# GSEA富集,logFC排序,從gene水平富集到hsa的通路
gsea_res=KO_gsea(da_res,weight = "logFC",type = "hsa",feature = "gene")
# GSEA結(jié)果如下：
head(gsea_res@result)
## ID
## hsa05235 hsa05235
## hsa05162 hsa05162
## hsa05323 hsa05323
## hsa05322 hsa05322
## hsa04660 hsa04660
## hsa04142 hsa04142
## Description
## hsa05235 PD-L1 expression and PD-1 checkpoint pathway in cancer - Homo sapiens (human)
## hsa05162 Measles - Homo sapiens (human)
## hsa05323 Rheumatoid arthritis - Homo sapiens (human)
## hsa05322 Systemic lupus erythematosus - Homo sapiens (human)
## hsa04660 T cell receptor signaling pathway - Homo sapiens (human)
## hsa04142 Lysosome - Homo sapiens (human)
## setSize enrichmentScore NES pvalue p.adjust qvalue rank
## hsa05235 13 0.8015443 1.885221 0.005680812 0.3204232 0.3178287 118
## hsa05162 18 0.7355862 1.857939 0.006161984 0.3204232 0.3178287 14
## hsa05323 14 0.7751069 1.850851 0.005438948 0.3204232 0.3178287 1
## hsa05322 15 0.7364826 1.787932 0.011670053 0.3798506 0.3767749 101
## hsa04660 15 0.7348533 1.783977 0.012174699 0.3798506 0.3767749 43
## hsa04142 20 -0.5462876 -1.628928 0.018884445 0.4909956 0.4870199 40
## leading_edge core_enrichment
## hsa05235 tags=46%, list=12%, signal=41% CD28/TICAM1/MAPK13/ALK/RPS6KB2/TLR9
## hsa05162 tags=11%, list=1%, signal=11% CD28/TAB2
## hsa05323 tags=7%, list=0%, signal=7% CD28
## hsa05322 tags=27%, list=10%, signal=24% CD28/H4C13/H3C15/H2BW2
## hsa04660 tags=13%, list=4%, signal=13% CD28/MAPK13
## hsa04142 tags=15%, list=4%, signal=15% AP1M2/SLC17A5/SLC11A2
ID：某條通路基因集的名字
Description：通路的文字描述
setSize：gene set（S）中的基因數(shù)目（經(jīng)過(guò)條件篩選后的值）
enrichmentScore：富集評(píng)分
NES：校正后的歸一化的ES值
pvalue：對(duì)富集得分ES的統(tǒng)計(jì)學(xué)分析，用來(lái)表征富集結(jié)果的可信度
p.adjust：對(duì)pvalue進(jìn)行BH調(diào)整后的值
qvalue：即多重假設(shè)檢驗(yàn)校正之后的pvalue。對(duì)NES可能存在的假陽(yáng)性結(jié)果的概率估計(jì)，因此FDR越小說(shuō)明富集越顯著；
rank：取到ES值時(shí)，對(duì)應(yīng)基因在排序好的基因列表中所處的位置
leading_edge：該處有3個(gè)統(tǒng)計(jì)值
-
ES>0在左邊，ES<0在右邊
tags=xx%表示peak gene在S中的位置。指示有貢獻(xiàn)的gene數(shù)量。
list=xx%表示peak gene在L中的位置。指示ES在哪里得到。
-
signal計(jì)算：
N：L中的gene數(shù)量
Nh：S中的gene數(shù)量
- 結(jié)果可視化
使用enrichplot包對(duì)富集結(jié)果進(jìn)行可視化。
#對(duì)于單條通路
enrichplot::gseaplot(gsea_res,geneSetID = "hsa05235")
#多條通路展示
enrichplot::gseaplot2(gsea_res,geneSetID = c("hsa05235","hsa05162","hsa04660"))
第1部分 - ES折線圖：在ES折線圖中，離垂直距離x=0軸最遠(yuǎn)的峰值便是基因集的ES值。峰出現(xiàn)在排序基因集的前端（ES值大于0）則說(shuō)明通路上調(diào)，出現(xiàn)在后端（ES值小于0）則說(shuō)明通路下調(diào)。
第2部分 - 基因集成員位置圖：在該圖中，用豎線標(biāo)記了基因集中各成員出現(xiàn)在基因排序列表中的位置。若豎線集中分布在基因排序列表的前端或后端，說(shuō)明該基因集通路上調(diào)或下調(diào)；若豎線較均勻分布在基因排序列表中，則說(shuō)明該基因集通路在比較的兩個(gè)數(shù)據(jù)中無(wú)明顯變化。 紅色部分對(duì)應(yīng)的基因在實(shí)驗(yàn)組中高表達(dá)，藍(lán)色部分對(duì)應(yīng)的基因在對(duì)照組中高表達(dá)， leading edge subset 是（0,0）到曲線峰值ES出現(xiàn)對(duì)應(yīng)的這部分基因成員。
第3部分 - 排序后所有基因rank值分布：該圖展示了排序后的所有基因rank值（由log2FoldChange值計(jì)算得出）的分布，以灰色面積圖顯展示。
Reference
links：
https://zhuanlan.zhihu.com/p/581172803
https://zhuanlan.zhihu.com/p/518144716
- A. Subramanian, P. Tamayo, V. K. Mootha, S. Mukherjee, et al., Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences. 102, 15545–15550 (2005).
- J. Reimand, R. Isserlin, V. Voisin, M. Kucera, et al., Pathway enrichment analysis and visualization of omics data using g:Profiler, GSEA, cytoscape and EnrichmentMap. Nature Protocols. 14, 482–517 (2019).
