2021-07-02 clusterProfiler分析GSEA

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])

参考： 生信补给站