GSEA 富集分析

生信喵

一、GSEA 简介
       Gene Set Enrichment Analysis (基因集富集分析)用来评估一个预先定义的基因集的基因在与表型相关度排序的基因表中的分布趋势，从而判断其对表型的贡献。其输入数据包含两部分，一是已知功能的基因集 (可以是 GO 注释、MsigDB 的注释或其它符合格式的基因集定义)，一是表达矩阵，软件会对基因根据其于表型的关联度(可以理解为表达值的变化)从大到小排序，然后判断基因集内每条注释下的基因是否富集于表型相关度排序后基因表的上部或下部，从而判断此基因集内基因的协同变化对表型变化的影响。
      
       GSEA

1. GSWA 网站：http://www.gsea-msigdb.org/gsea/index.jsp
   
2. JavaGSEA：
   
3. http://www.gsea-msigdb.org/gsea/login.jsp;jsessionid=473B7BFEFD0F6DF30707EF9D51F73148
   
4. 教程：https://enrichmentmap.readthedocs.io/en/docs-2.2/Tutorial_GSEA.html
   
5. https://baderlab.github.io/Cytoscape_workflows/EnrichmentMapPipeline/
   
6. 案例数据：https://www.gsea-msigdb.org/gsea/datasets.jsp
   
7. 注释文件下载：http://www.gsea-msigdb.org/gsea/msigdb/index.jsp
   

复制代码

二、GSEA 原理
       给定一个排序的基因表 L 和一个预先定义的基因集 S (比如编码某个代谢通路的产物的基因,基因组上物理位置相近的基因，或同一 GO 注释下的基因)，GSEA 的目的是判断 S 里面的成员 s 在 L 里面是随机分布还是主要聚集在 L 的顶部或底部。这些基因排序的依据是其在不同表型状态下的表达差异，若研究的基因集 S 的成员显著聚集在 L 的顶部或底部，则说明此基因集成员对表型的差异有贡献，也是我们关注的基因集。
      
       GSEA 原理
GSEA 计算中几个关键概念：
       1、计算富集得分 (ES, enrichment score). ES 反应基因集成员 s 在排序列表 L 的两端富集的程度。计算方式是，从基因集 L 的第一个基因开始，计算一个累计统计值。当遇到一个落在 s里面的基因，则增加统计值。遇到一个不在 s 里面的基因，则降低统计值。
       每一步统计值增加或减少的幅度与基因的表达变化程度（更严格的是与基因和表型的关联度，可能是 fold-change，也可能是 pearson corelation 值，后面有介绍几种不同的计算方式）是相关的，可以是线性相关，也可以是指数相关 （具体见后面参数选择）。富集得分 ES 最后定义为最大的峰值。正值 ES 表示基因集在列表的顶部富集，负值 ES 表示基因集在列表的底部富集。
       2、评估富集得分(ES)的显著性。通过基于表型而不改变基因之间关系的排列检验(permutation test)计算观察到的富集得分(ES)出现的可能性。若样品量少，也可基于基因集做排列检验，计算 p-value。
       3、多重假设检验校正。首先对每个基因子集 s 计算得到的 ES 根据基因集的大小进行标准化得到 Normalized Enrichment Score (NES)。随后针对 NES 计算假阳性率。（计算 NES 也有另外一种方法，是计算出的 ES 除以排列检验得到的所有 ES 的平均值）
       4、Leading-edge subset，对富集得分贡献最大的基因成员。

三、MSigDB
       MSigDB（Molecular Signatures Database）数据库是 GSEA 官方提供的基因列表数据库。
      
       MSigDB 数据库分类
目前 MSigDB 分为九个等级，分别为
       H：癌症特征基因集合，共 50 gene sets，最常用。
       C1：染色体位置基因集合，共 299 gene sets，用的很少。
       C2：包含了已知数据库，文献和专家支持的基因集信息，包含 5529 gene sets。
       C3：调控靶基因集合，包括 miRNA-target 以及转录因子-target 调控模式，3735 gene sets。
       C4：计算机软件预测出来的基因集合，主要是和癌症相关的基因，858 gene sets。
       C5：Gene Ontology 对应的基因集合，10192 gene sets。
       C6：致癌基因集合，大部分来源于 NCBI GEO 发表芯片数据，189 gene sets。
       C7：免疫相关基因集，4872 gene sets。
       C8：single cell identitiy gene sets, 302 gene sets
       MSgIDB 格式非常简单，里面是基因名字的一个列表，为 gmt 格式。


四、输入文件
       GSEA 分析一般需要三种文件，分别是基因表达文件，是一个矩阵格式，扩展名可以是*gct,*res,*pcl 或者*txt。第二个是分组信息，需要手工填写，扩展名为*.cls。最后就是注释文件，为 gmt 或者 gmx 格式，可以直接从 MSigDB（Molecular Signatures Database）数据库下载，也可以自己准备。
详细文件格式见下面链接：

1. https://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats

复制代码

之前的推文也有介绍，详细参考：

1. https://mp.weixin.qq.com/s/l-oFTtSqCjHp94LcZu2N_g

复制代码

五、利用 R 实现 GSEA
       虽然 GSEA 客户端可以非常方便的完成 GSEA 的分析，但 JAVA 版的 GSEA 软件图形输出格式是 png 格式，图片分辨率较低，不方便编辑，所以，也可以使用 R 包进行 GSEA 分析。能做GSEA 分析的 R 包有很多，例如 clusterProfiler，GSEAbase，fgsea 等。

1. library(GSEABase)
   
2. library(clusterProfiler)
   
3. library(org.Hs.eg.db)
   
4. 
   
5. #使用示例基因排序列表，之前的差异分析结果，这里就不分享了，大家跑下之前的示例RNAseq差异表达就行
   
6. x <- read.csv("res.csv",row.names = 1)
   
7. head(x)
   
8. 
   
9. geneList <- x$log2FoldChange
   
10. names(geneList) <- toupper(rownames(x))
    
11. geneList <- sort(geneList, decreasing = T)
    
12. head(geneList)
    
13. names(geneList) <- mapIds(org.Hs.eg.db, keys=names(geneList), column="SYMBOL", keytype="ENSEMBL",multiVals = "first")
    
14. head(geneList)
    
15. #去除NA
    
16. geneList <- geneList[!is.na(names(geneList))]#20505
    
17. 
    
18. #读取gmt文件,从官网下载
    
19. gmtfile <- "msigdb.v2022.1.Hs.symbols.gmt"
    
20. geneset <- read.gmt(gmtfile)
    
21. 
    
22. #GSEA分析
    
23. egmt <- GSEA(geneList, TERM2GENE = geneset, minGSSize = 1, pvalueCutoff = 0.99, verbose = F)
    
24. egmt <- GSEA(geneList, TERM2GENE = geneset, minGSSize = 1, pvalueCutoff = 0.05, verbose = F, eps=0)
    
25. #查看结果
    
26. gsea.out.df <- egmt@result
    
27. View(head(gsea.out.df))
    
28. gsea.out.df$ID
    
29. 
    
30. #绘制GSEA图
    
31. library(enrichplot)
    
32. options(repr.plot.width=6,repr.plot.height=4)
    
33. gseaplot2(egmt, geneSetID = "CHEN_LVAD_SUPPORT_OF_FAILING_HEART_UP", pvalue_table = T)

复制代码

      
       GSEA 可视化

1. library(msigdbr)
   
2. library(fgsea)
   
3. library(tidyverse)
   
4. library(clusterProfiler)
   
5. library(org.Hs.eg.db)
   
6. 
   
7. ## 预定义基因集
   
8. mdb_c2 <- msigdbr(species = "Homo sapiens")
   
9. mdb_c2
   
10. mdb_kegg = mdb_c2 [grep("^KEGG",mdb_c2$gs_name),]
    
11. mdb_kegg
    
12. 
    
13. fgsea_sets <- mdb_kegg %>% split(x = .$gene_symbol, f = .$gs_name)
    
14. fgsea_sets
    
15. 
    
16. #读取基因差异表达文件
    
17. x <- read.csv("res.csv",row.names = 1)
    
18. head(x)
    
19. 
    
20. #基因按logFC排序
    
21. rownames(x)
    
22. x$genes <- mapIds(org.Hs.eg.db, keys=rownames(x), column="SYMBOL", keytype="ENSEMBL",multiVals = "first")
    
23. x <- na.omit(x)
    
24. gene_log2 <- x %>% arrange(desc(log2FoldChange)) %>% dplyr::select(genes,log2FoldChange)
    
25. 
    
26. generanks<- deframe(gene_log2)
    
27. generanks
    
28. 
    
29. #fgsea分析
    
30. fgseaRes <- fgsea(pathways = fgsea_sets,
    
31.                   stats = generanks,
    
32.                   minSize=15,
    
33.                   maxSize=500,
    
34.                   nperm=10000)
    
35. 
    
36. #fgsea绘图
    
37. plotEnrichment(fgsea_sets[["KEGG_P53_SIGNALING_PATHWAY"]],generanks) +
    
38.   labs(title="KEGG_P53_SIGNALING_PATHWAY")

复制代码