理論知識
GSEA（Gene Set Enrichment Analysis），該方法發(fā)表于2005年的Gene set enrichment analysis: a knowledge-based approach forinterpreting genome-wide expression profiles，是一種基于基因集的富集分析方法。
應(yīng)用該分析的背景：常規(guī)的轉(zhuǎn)錄組分析會獲得大量的差異基因，如何將巨量的差異基因與生物學功能結(jié)合在一起，就成了一個問題。當然，后面出現(xiàn)了KEGG 和GO分析，這些分析可以說比較籠統(tǒng)，緊緊將所有的差異基因序列提取出來，與預(yù)設(shè)的通路比較，得到通路富集的結(jié)果，至于上調(diào)或下調(diào)則不得而知。此外，僅僅關(guān)注單個基因的差異表達，并不能說其參與的整個通路受到什么影響？相對而言，整個通路的所有基因的整體上調(diào)或下調(diào)才會顯得更有意義。
為了解決這一問題，就出現(xiàn)了GSEA分析。
GSEA分析首先根據(jù)FoldChange值排序（一般是降序），然后根據(jù)預(yù)設(shè)的生物學功能的gene subset進行注釋，就很容易獲得參與某條通路/途徑的關(guān)鍵基因，或者直接看出通路調(diào)控整體是上升或下降。
在對基因表達數(shù)據(jù)分析時，首先確定分析的目的，即選擇MSigDB中的一個或多個功能基因集進行分析，然后基于基因表達數(shù)據(jù)與表型的關(guān)聯(lián)度（也可以理解為表達量的變化）的大小進行排序。然后判斷每個基因集內(nèi)的基因是否富集于表型相關(guān)度排序后基因列表的上部或下部，從而判斷此基因集內(nèi)基因的協(xié)同變化對表型變化的影響。
GSEA的輸入是一個基因表達量矩陣，其中的樣本分成了A和B兩組，找到兩組之間差異表達的基因，然后根據(jù)foldchange進行排序，用來表示基因在兩組間表達量的變化趨勢。Gene list Rank 可以理解成foldchange由大到小，由正到負排序。 排序之后的基因列表其左側(cè)可以看做是上調(diào)的差異基因，其右側(cè)是下調(diào)的差異基因，中間屬于差異不顯著的基因。GSEA分析的是一個基因集下的所有基因是否在這個排序列表的頂部或者底部富集，如果在頂部富集，我們可以說，從總體上看，該基因集是上調(diào)趨勢，反之，如果在底部富集，則是下調(diào)趨勢。(建議讀原文獻)
關(guān)鍵點
- 首先根據(jù)自己的實驗?zāi)康倪x擇合適的gmt文件，可以參考[基因集數(shù)據(jù)庫MSigdb] (GSEA | MSigDB (gsea-msigdb.org)。image.png
我主要關(guān)注的基因集數(shù)據(jù)庫是：
- H: hallmark gene sets
- C7: immuologic signature gene sets
- C8: cell type signature gene sets.
-
基因表達矩陣
主要使用前兩列（Gene和log2FoldChange）
image.png
代碼實戰(zhàn)
#準備gmt文件
library(tidyverse)
library(org.Hs.eg.db)
library(clusterProfiler)
msigdb_GMTs <- "msigdb_v7.0_GMTs"
msigdb <- "h.all.v7.5.1.entrez.gmt" #根據(jù)個人實際需求下載或編輯
kegmt <- read.gmt(file.path(msigdb_GMTs,msigdb))
#準備差異基因列表
geneList = DEG$log2FoldChange
names(geneList) = DEG$ENTREZID
head(geneList)
geneList = sort(geneList, decreasing = TRUE)
#GSEA分析
set.seed(1)
egmt<-GSEA(geneList,TERM2GENE = kegmt) #GSEA分析
#轉(zhuǎn)換成數(shù)據(jù)框
egmt_result_df <- as.data.frame(egmt)
write.table(egmt_result_df,file="GSEA_MSigDb_h.all_result.txt",sep = "\t",row.names = T,col.names = NA,quote = F)
save(egmt,egmt_result_df,file = "GSEA_deg_h.all.rda")
可視化
gseaplot2(egmt,3,color="red",pvalue_table = T)
#結(jié)果匯總
gseaplot2(egmt, geneSetID = c(1,3), subplots = 1:3,pvalue_table = T)
#氣泡圖 展示geneset被激活還是抑制
egmt2<- setReadable(egmt,OrgDb=org.Hs.eg.db, keyType = "ENTREZID")
dotplot(egmt2,split=".sign",showCategory = 10,font.size = 10, title = "",
label_format = 30,)+facet_grid(~.sign)
#edit legends
# +guides(
#reverse color order (higher value on top)
#color = guide_colorbar(reverse = TRUE))
#reverse size order (higher diameter on top)
#size = guide_legend(reverse = TRUE))
# Title 可以添加標題
結(jié)果輸出：
-
單個通路
單個通路整體下調(diào)
2.多個通路同時展示
hallmark gene set 注釋分析
從這個圖可以看出處理組與對照組比較，哪些通路被激活哪些通路被抑制，還可進一步去看單個通路中起到關(guān)鍵作用的基因集：Leading Edge Subset。 針對這一塊進一步的分析我還不會，這也是深入研究比較關(guān)鍵的部分。后續(xù)再補充吧。
image.png
