知识分享| 转录组个性化分析(2)——GSEA
上一期中给大家介绍了拟时序分析的意义及具体的分析过程，本期继续给大家带来转录组个性化分析——GSEA。废话不多说，干货直接奉上！
1 GSEA基本概念
基因集富集分析（Gene Set Enrichment Analysis，GSEA）：用一个预先定义的基因集中的基因来评估在与表型相关度排序的基因表中的分布趋势，从而判断其对表型的贡献。
2 GSEA原理
(1) 背景基因排序：将全部基因按照某种指标（差异分析p值，表型相关性，表达量等）进行排序，比如log2FC排序。
(2) 目标基因富集：将某个特定类型的基因在排序表中标出，目标基因可以是某个通路或GO terms的基因等。
(3) 计算富集分数：使用加权法，计算ES值变化。对位于中部（与性状相关性低）的部分采用较小的权值，所以越集中在两端，与表型的相关性越高。ES曲线最大值为富集分数（Enrichment Score）。
(4) Permutation test：对基因集的ES值进行显著性检验及多重假设检验，从而计算出显著富集的基因集。
3 GSEA分析的作用
GSEA和常规的GO、KEGG的差异在于，GSEA使用的是基因集，传统的富集分析不需要考虑基因表达量的变化趋势，其算法的核心只关注这些差异基因的分布是否和随机抽样得到的分布一致，即使后期在可视化时，我们在通路图上用不同颜色标记了上下调的基因，但是由于没有采用有效的统计学手段去分析这条通路下所有差异基因的总体变化趋势，这使得传统的富集分析结果无法回答如下问题：一个富集到的通路下，既有上调差异基因，也有下调差异基因，那么这条通路总体的表现形式究竟是怎样呢，是被抑制还是激活？或者更直观点说，这条通路下的基因表达水平在实验处理后是上升了呢，还是下降了呢？
想要知道进行了差异分析的两组有什么功能和通路的差别，手上有大部分的功能分子以及对应的值，这个值可以是 logFC。可以用这个 logFC作为分子的排序，从而来评估在预先定义的基因集中是否显著富集。
4 GSEA分析代码
library(DO.db)
require(DOSE)
library(clusterProfiler)
library(AnnotationHub)
library(readr)
library("genefilter")
library(pheatmap)
library(tidyverse)
library(DESeq2)
library(ggplot2)
library(export)
library(enrichplot)
library(Rgraphviz)
#构造图片输出函数，need input filename width height
#函数依赖export包
out_img <- function(filename,pic_width=5,pic_height=7){
graph2png(file=filename,width=pic_width,height=pic_height)
graph2ppt(file=filename,width=pic_width,height=pic_height)
graph2tif(file=filename,width=pic_width,height=pic_height)
}
#all_entrez.csv的第一列是entrezid,第二列是FoldChange的值。
#gse需要单独做数据格式
d <- read.csv("all_entrez.csv")
geneList <- d[,2]
names(geneList)=as.factor(d[,1])
geneList <- sort(geneList,decreasing = TRUE)
#gseGO进行GSEA分析
#参考连接https://yulab-smu.github.io/clusterProfiler-book/chapter12.html
###gseBP <- gseGO(geneList=geneList,ont="BP",OrgDb=maize,keyType = 'ENTREZID',nPerm = 50000,minGSSize = 100,maxGSSize = 6000,pvalueCutoff = 0.05,verbose = FALSE)
############# GSEA CC 模式 start
ego3 <- gseGO(geneList = geneList,OrgDb = maize,ont = "CC",nPerm = 1000,minGSSize = 100,maxGSSize = 1000,pvalueCutoff = 0.05,verbose = FALSE)
write.csv(ego3,file = "GESA-GO_CC.csv")
#ridgeline plot for expression distribution of GSEA result
ridgeplot(ego3)
out_img(filename = "ridgeplot_CC",pic_width = 12,pic_height = 12)
#只显示值最高的一组的信息
#gseaplot(ego3,geneSetID = 1,by="runningScore",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,by="preranked",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,title=ego3$Description[1])
#显示前4组信息
gseaplot2(ego3,geneSetID = 1:4, ES_geom = "dot",pvalue_table = TRUE)
out_img(filename = "gseaplot_CC",pic_width=12,pic_height = 10)
#gsearank(ego3,1,title=ego3[1,"Description"])
############GSEA CC 模式end
############# GSEA BP 模式 start
ego2 <- gseGO(geneList = geneList,OrgDb = maize,ont = "BP",pvalueCutoff = 0.05,verbose = FALSE)
write.csv(ego2,file = "GESA-GO_BP.csv")
#ridgeline plot for expression distribution of GSEA result
ridgeplot(ego2)
out_img(filename = "ridgeplot_BP",pic_width = 12,pic_height = 12)
#只显示值最高的一组的信息
#gseaplot(ego3,geneSetID = 1,by="runningScore",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,by="preranked",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,title=ego3$Description[1])
#显示前4组信息
gseaplot2(ego2,geneSetID = 1:4, ES_geom = "dot",pvalue_table = TRUE)
out_img(filename = "gseaplot_BP",pic_width=12,pic_height = 10)
#gsearank(ego3,1,title=ego3[1,"Description"])
############GSEA BP 模式end
############# GSEA MF 模式 start
ego4 <- gseGO(geneList = geneList,OrgDb = maize,ont = "MF",pvalueCutoff = 0.05,verbose = FALSE)
write.csv(ego4,file = "GESA-GO_MF.csv")
#ridgeline plot for expression distribution of GSEA result
ridgeplot(ego4)
out_img(filename = "ridgeplot_MF",pic_width = 12,pic_height = 12)
#只显示值最高的一组的信息
#gseaplot(ego3,geneSetID = 1,by="runningScore",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,by="preranked",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,title=ego3$Description[1])
#显示前4组信息
gseaplot2(ego4,geneSetID = 1:4, ES_geom = "dot",pvalue_table = TRUE)
out_img(filename = "gseaplot_MF",pic_width=12,pic_height = 10)
#gsearank(ego3,1,title=ego3[1,"Description"])
############GSEA MF 模式end
#gsaKEGG基因富集分析
kk2 <- gseKEGG(geneList = geneList,organism = 'zma',pvalueCutoff = 0.05,verbose = FALSE)
write.csv(kk2,file="GSEA_KEGG.csv")
gseaplot2(kk2,geneSetID = 1:4,ES_geom = "dot",pvalue_table = TRUE)
out_img(filename="GSEA_KEGG",pic_width = 12,pic_height = 10)
ridgeplot(kk2)out_img(filename="ridgeplot_GSEA_KEGG",pic_width = 12,pic_height = 12)
5 结果解读
典型结果图由上、中、下三个部分组成：
上：为富集评分的情况，如果 NES 为正，则峰出现在左侧（头部富集）（高表达组富集）基因集中核心分子主要集中在左侧高表达组中；如果 NES 为负，则尾部会出现谷（尾部富集）（低表达组富集），基因集中核心分子主要集中在右侧低表达组中。
中：每一根竖线代表基因集中一个分子，上传数据的分子根据给定的值进行排序，排序后单独提取当前基因集中的定义的分子，分子的位置情况即为中间部分的所示。
下：把上传数据分子给定的值进行归一化后的值进行可视化。下部分的结果可以不用怎么关注。
从该图中可以看出，这个基因集是在MUT这一组高表达的，在折线图中有个峰值，该峰值就是这个基因集的Enrichemnt score，峰值之前的基因就是该基因集下的核心基因。基因集下的所有基因在这个排序列表的顶部富集，我们可以说，从总体上看，该基因集是上调趋势，反之，如果在底部富集，则是下调趋势。
