理论知识
GSEA（Gene Set Enrichment Analysis），该方法发表于2005年的Gene set enrichment analysis: a knowledge-based approach forinterpreting genome-wide expression profiles，是一种基于基因集的富集分析方法。
应用该分析的背景：常规的转录组分析会获得大量的差异基因，如何将巨量的差异基因与生物学功能结合在一起，就成了一个问题。当然，后面出现了KEGG 和GO分析，这些分析可以说比较笼统，紧紧将所有的差异基因序列提取出来，与预设的通路比较，得到通路富集的结果，至于上调或下调则不得而知。此外，仅仅关注单个基因的差异表达，并不能说其参与的整个通路受到什么影响？相对而言，整个通路的所有基因的整体上调或下调才会显得更有意义。
为了解决这一问题，就出现了GSEA分析。
GSEA分析首先根据FoldChange值排序（一般是降序），然后根据预设的生物学功能的gene subset进行注释，就很容易获得参与某条通路/途径的关键基因，或者直接看出通路调控整体是上升或下降。
在对基因表达数据分析时，首先确定分析的目的，即选择MSigDB中的一个或多个功能基因集进行分析，然后基于基因表达数据与表型的关联度（也可以理解为表达量的变化）的大小进行排序。然后判断每个基因集内的基因是否富集于表型相关度排序后基因列表的上部或下部，从而判断此基因集内基因的协同变化对表型变化的影响。
GSEA的输入是一个基因表达量矩阵，其中的样本分成了A和B两组，找到两组之间差异表达的基因，然后根据foldchange进行排序，用来表示基因在两组间表达量的变化趋势。Gene list Rank 可以理解成foldchange由大到小，由正到负排序。 排序之后的基因列表其左侧可以看做是上调的差异基因，其右侧是下调的差异基因，中间属于差异不显著的基因。GSEA分析的是一个基因集下的所有基因是否在这个排序列表的顶部或者底部富集，如果在顶部富集，我们可以说，从总体上看，该基因集是上调趋势，反之，如果在底部富集，则是下调趋势。(建议读原文献)
关键点
- 首先根据自己的实验目的选择合适的gmt文件，可以参考[基因集数据库MSigdb] (GSEA | MSigDB (gsea-msigdb.org)。image.png
我主要关注的基因集数据库是：
- H: hallmark gene sets
- C7: immuologic signature gene sets
- C8: cell type signature gene sets.
-
基因表达矩阵
主要使用前两列（Gene和log2FoldChange）
image.png
代码实战
#准备gmt文件
library(tidyverse)
library(org.Hs.eg.db)
library(clusterProfiler)
msigdb_GMTs <- "msigdb_v7.0_GMTs"
msigdb <- "h.all.v7.5.1.entrez.gmt" #根据个人实际需求下载或编辑
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分析
#转换成数据框
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)
#结果汇总
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 可以添加标题
结果输出：
-
单个通路
单个通路整体下调
2.多个通路同时展示
hallmark gene set 注释分析
从这个图可以看出处理组与对照组比较，哪些通路被激活哪些通路被抑制，还可进一步去看单个通路中起到关键作用的基因集：Leading Edge Subset。 针对这一块进一步的分析我还不会，这也是深入研究比较关键的部分。后续再补充吧。
image.png
