GSEA（Gene Set EnrichmentAnalysis），即基因集富集分析，无需设定阈值来区分上调下调基因，使用所有的基因进行分析。
GSEA分析只需要两列信息，基因列和logFC（不同软件的差异分析这一列的名字会有差别）。
首先构建一个genelist，可以是来自自己测序数据差异分析的结果或者是GEO数据集，
genelist由两列构成，第一列表示差异表达的基因ID（基因ID不能重复的,形式同样为entrezID），
第二列为基因对应的表达量或者是FC等数值型向量，注意按照数值从高到低排列
基因名是symbol，要将之转换为entrezid格式
导入数据
library(clusterProfiler)
library(ggplot2)
setwd( "D:/GEO数据挖掘与meta分析/练习/24.GSEA R分析及画图(代码)/24.GSEA R分析及画图(代码)")
deg <- read.csv("diff_ENTREZID.csv", as.is = T)
deg <-na.omit(deg)
geneList <- deg$logFC
names(geneList) <- deg$ENTREZID
head(geneList)#结果为排序的logGC，names为ENTREZID
3040 3043 2354 2219 3553 3576
4.166594 4.045908 2.531847 2.507169 2.138769 1.913179
GSEA分析
需要gmt文件，http://www.gsea-msigdb.org/gsea/downloads.jsp路径下载，选择合适的
gmtfile <- system.file("extdata", "c5.cc.v5.0.entrez.gmt", package="clusterProfiler")
S<-read.gmt(gmtfile) ##读取gmt文件得到基因集S
egmt <-GSEA(geneList,TERM2GENE=S,pvalueCutoff=0.5
#verbose=FALSE #是否打印信息
) #GSEA
head(egmt)
可视化-点图
dotplot(egmt)
GSEA-GO分析
library(org.Hs.eg.db)
go <- gseGO(geneList,OrgDb=org.Hs.eg.db,pvalueCutoff=0.5,
nPerm=1000#置换检验的置换次数
#minGSSize =100, #用于测试的功能集最小容量
#maxGSSize =500, #用于测试的功能集最大容量
)
write.csv(go,"gseGO.csv",quote = F)
GSEA-KEGG分析
kk <- gseKEGG(geneList,organism = "hsa",nPerm=1000)
write.csv(kk,"gseKEGG.csv",quote = F)
gse.KEGG <- gseKEGG(geneList,
organism = "hsa", # 人 hsa
pvalueCutoff = 1,
pAdjustMethod = "BH",) #具体参数在下面
head(gse.KEGG)
head(gse.KEGG)[1:10]
可视化
#The ridgeplot will visualize expression distributions of core enriched genes for GSEA enriched categories. It helps users to interpret up/down-regulated pathways.
ridgeplot(go)
ggsave(file="ridgeplot.pdf",width= 18,height = 10)
gseaplot2(go,geneSetID=1,title=egmt$Description[1],pvalue_table=T)
展示多个GSEA结果
使用数字的方式
gseaplot2(go,
1:3, #绘制前3个
pvalue_table = T) # 显示p值
使用向量指定通路
gseaplot2(go,
c("PLASMA_MEMBRANE","MEMBRANE"), #指定通路向量
pvalue_table = T) # 显示p值
点 形式
gseaplot2(go,
1:5, #按照第一个作图
ES_geom = "dot",
base_size = 20,
pvalue_table = T)
gseaplot(go, geneSetID = 1, by = "runningScore", title = go$Description[1])
gseaplot(go, geneSetID = 1, by = "preranked", title = go$Description[1])
gseaplot(go, geneSetID = 1, title = go$Description[1])
参考： 生信补给站
