1 什么是GSEA
基因集富集分析(Gene Set Enrichment Analysis, GSEA)是是一种计算方法，用于确定事先定义的一组基因是否在不同的样品中差异表达。
GSEA首先将基因在样品中的差异倍数值（logFC）由大到小排序，然后判断来自功能注释等预定义的基因集或自定义的基因集中的基因是富集在这个排序列表的顶部还是底部，如果在富集顶部，则该基因集是上调趋势，反之，如果富集在底部，则是下调趋势。
GSEA官网提供了详细说明，以及对应软件的下载地址。
2 GSEA特点
传统的KEGG（通路富集分析）和GO（功能富集）分析时，针对总体的差异基因，不区分哪些差异基因是上调还是下调。所以会出现同一通路下富集到的的既有上调差异基因，也有下调差异基因，无法判断这条通路表现形式究竟是怎样。
而GSEA考虑了基因的表达水平，不需要明确指定差异基因阈值，检验的是基因集而非单个基因的表达变化，算法会根据实际数据的整体趋势进行分析，以判断这条通路的表达情况，激活或者抑制。
3 GSEA结果解读
- 第1部分：Enrichment Score的折线图
横轴为排序后的基因，纵轴为对应的ES, 在折线图中出现的峰值就是这个基因集的富集分数（Enrichment Score，ES）。ES是从排序后的表达基因集的第一个基因开始，如果排序后表达基因列表中的基因出现在功能基因数据集中则加分，反之则减分。正值说明在顶部富集，峰值左边的基因为核心基因，负值则相反。
- 第2部分：基因位置图
黑线代表排序后表达基因列表中的基因位于当前分析的功能注释基因集的位置，红蓝相间的热图是表达丰度排列，红色越深的表示该位置的基因logFC越大 ，蓝色越深表示logFC越小。如果研究的功能注释基因集的成员显著聚集在表达数据集的顶部或底部，则说明功能基因数据集中的基因在数据集中高表达或低表达，若随机分配，则说明表达数据集与该通路无关。
- 第3部分：每个基因对应的信噪比（Signal2noise）
以灰色面积图展示。灰色阴影的面积比，可以从整体上反映组间的Signal2noise的大小。
一般认为校正后的富集分数（NES）|NES|>1,p<0.05, qvalue(即FDR)<0.25的结果有意义。
4 GSEA实战
#加载数据 数据链接：https://wwu.lanzouf.com/idmJB07bcefa
load(file = 'GSEA_test.Rdata')
colnames(deg)
head(deg)
#得到差异基因列表后取出 ，p值和logFC
nrDEG=deg[,c(2,1)]
colnames(nrDEG)=c('log2FoldChange','pvalue')
head(nrDEG)
log2FoldChange pvalue
LYZ -1.812078 1.228136e-144
FCGR3A 2.617707 3.977801e-119
S100A9 -2.286734 2.481257e-95
S100A8 -2.610696 3.626489e-92
IFITM2 1.445772 7.942512e-87
LGALS2 -2.049431 1.275856e-75
#注：最后需要为文件如上：一列为基因名，一列为FC值，一列为p值
#加载Y叔的R包，把SYMBOL转换为ENTREZID，后面可以直接做 KEGG 和 GO
library(org.Hs.eg.db)
library(clusterProfiler)
gene <- bitr(rownames(nrDEG), #转换的列是nrDEG的列名
fromType = "SYMBOL", #需要转换ID类型
toType = "ENTREZID", #转换成的ID类型
OrgDb = org.Hs.eg.db) #对应的物种，小鼠的是org.Mm.eg.db
#让基因名、ENTREZID、logFC对应起来
gene$logFC <- nrDEG$log2FoldChange[match(gene$SYMBOL,rownames(nrDEG))]
head(gene)
geneList=gene$logFC
names(geneList)=gene$ENTREZID
#按照logFC的值来排序geneList
geneList=sort(geneList,decreasing = T)
head(geneList)
以上即完成了数据的准备工作，开始进行GSEA分析
GSEA-KEGG
library(clusterProfiler)
#clusterProfiler包的妙用
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 10,
pvalueCutoff = 0.9,
verbose = FALSE)
#提取GSEA-KEGG结果
tmp=kk_gse@result
kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
tmp=kk@result
pro='test_gsea'
write.csv(kk@result,paste0(pro,'_kegg.gsea.csv'))
#按照enrichment score从高到低排序，便于查看富集通路
kk_gse=kk
sortkk<-kk_gse[order(kk_gse$enrichmentScore, decreasing = T),]
head(sortkk)
dim(sortkk)
#write.table(sortkk,"gsea_output2.txt",sep = "\\t",quote = F,col.names = T,row.names = F)
#可以根据自己想要的通路画出需要的图
library(enrichplot)
gseaplot(kk_gse, "hsa04510")
gseaplot2(kk_gse, "hsa04510", color = "firebrick",
rel_heights=c(1, .2, .6))#改变更多参数，为了美观
#同时展示多个pathways的结果。
#画出排名前四的通路
gseaplot2(kk_gse, row.names(sortkk)[1:4])
#上图用的是ES排名前4个画图，也可以用你自己感兴趣的通路画图
# 自己在刚才保存的txt文件里挑选就行。
paths <- c("hsa04510", "hsa04512", "hsa04974")
paths <- row.names(sortkk)[5:8]
paths
gseaplot2(kk_gse, paths)
#这里的GSEA分析其实由三个图构成，GSEA分析的runningES折线图
# 还有下面基因的位置图，还有所谓的蝴蝶图。如果不想同时展示，还可以通过subplots改变。
gseaplot2(kk_gse, paths, subplots=1)#只要第一个图
gseaplot2(kk_gse, paths, subplots=1:2)#只要第一和第二个图
gseaplot2(kk_gse, paths, subplots=c(1,3))#只要第一和第三个图
#如果想把p值标上去，也是可以的。
gseaplot2(kk_gse, paths, color = colorspace::rainbow_hcl(4),
pvalue_table = TRUE)
#最后的总结代码
gseaplot2(kk_gse,#数据
row.names(sortkk)[1:5],#画哪一列的信号通路
title = "Prion disease",#标题
base_size = 15,#字体大小
color = "green",#线条的颜色
pvalue_table = TRUE,#加不加p值
ES_geom="line")#是用线，还是用d点
###当然，这里不知道具体需要什么通路，就全部都画出来
# 这里找不到显著下调的通路，可以选择调整阈值，或者其它。
down_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore < -0.4,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore > 0.4,];up_kegg$group=1
dat=rbind(up_kegg,down_kegg)
colnames(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group
dat=dat[order(dat$pvalue,decreasing = F),]
library(ggplot2)
library(enrichplot)
gesa_res=kk@result
###画出每张kegg通路
lapply(1:nrow(down_kegg), function(i){
gseaplot2(kk,down_kegg$ID[i],
title=down_kegg$Description[i],pvalue_table = T)
ggsave(paste0(pro,'_down_kegg_',
gsub('/','-',down_kegg$Description[i])
,'.pdf'))
})
lapply(1:nrow(up_kegg), function(i){
gseaplot2(kk,up_kegg$ID[i],
title=up_kegg$Description[i],pvalue_table = T)
ggsave(paste0(pro,'_up_kegg_',
gsub('/','-',up_kegg$Description[i]),
'.pdf'))
})
GSEA-GO
GO和KEGG分析流程基本相同，除了函数名和变量名的变化
ego <- gseGO(geneList = geneList,
OrgDb = org.Hs.eg.db,
ont = "ALL",
nPerm = 1000, ## 排列数
minGSSize = 100,
maxGSSize = 500,
pvalueCutoff = 0.9,
verbose = FALSE) ## 不输出结果
go=DOSE::setReadable(ego, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
go=ego@result
pro='test_gsea'
go_gse=go
sortgo<-go_gse[order(go_gse$enrichmentScore, decreasing = T),]
head(sortgo)
dim(sortgo)
#write.table(sortkk,"gsea_output2.txt",sep = "\\t",quote = F,col.names = T,row.names = F)
#可以根据自己想要的通路画出需要的图
library(enrichplot)
gseaplot2(go_gse,#数据
row.names(sortgo)[1:5],#画那一列的信号通路
title = "Prion disease",#标题
base_size = 15,#字体大小
color = "green",#线条的颜色
pvalue_table = TRUE,#加不加p值
ES_geom="line")#是用线，还是用d点
write.csv(go,file = 'gse_go.csv')
其他可视化方法
气泡图
dotplot(kk,split=".sign")+facet_wrap(~.sign,scales="free")
条形图
down_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore < -0.4,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore > 0.4,];up_kegg$group=1
dat=rbind(up_kegg,down_kegg)
colnames(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group
dat=dat[order(dat$pvalue,decreasing = F),]
library(ggplot2)
g_kegg<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) +
geom_bar(stat="identity") +
scale_fill_gradient(low="blue",high="red",guide = FALSE) +
scale_x_discrete(name ="Pathway names") +
scale_y_continuous(name ="log10P-value") +
coord_flip() + theme_bw(base_size = 15)+
theme(plot.title = element_text(hjust = 0.5), axis.text.y = element_text(size = 15))+
ggtitle("Pathway Enrichment")
g_kegg
网络图
library(ggplot2)
library(enrichplot)
cnetplot(kk,showCategory= 5, foldChange= geneList, colorEdge="TRUE")
#colorEdge不同的term展示不同的颜色，如果希望标记节点的子集，可以使用node_label参数，它支持4种可能的选择(即“category”、“gene”、“all”和“none”).
ego2<-pairwise_termsim(ego)
emapplot(ego2, showCategory= 10, color="pvalue", cex_category=1, layout="kk")
#cex_category调节节点大小，showCategory展示条目个数
参考来源
#section 3已更新#「生信技能树」单细胞公开课2021_哔哩哔哩_bilibili
致谢
I thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
THE END
网友评论