参考：

GSEA和GSVA，再也不用去下载gmt文件咯 - 简书 (jianshu.com)

(104条消息) 一个R包完成单细胞基因集富集分析 （全代码）_悟道西方-CSDN博客

#singleseqgset | 单细胞RNA-Seq基因集富集分析

#叫msigdbr的包，就是把MsigDB的那些基因集合都R语言化了。那么以后想做GSEA和GSVA，就可以不用去下载gmt文件啦！

#清空环境变量

rm(list = ls())

#########################################1.端详msigdbr这个包

library(msigdbr)

msigdbr_species()#Mus musculus,Homo sapiens,囊括了11个物种

#看小鼠有哪些基因集

#table(mouse[,1])

#H: hallmark gene sets;C1: positional gene sets; C2: curated gene sets.KEGG ;C3: motif gene sets;C4: computational gene sets

#C5: GO gene sets;C6: oncogenic signatures;C7: immunologic signatures

#############################################2.获取所需基因集

mouse<-msigdbr(species = "Mus musculus",category = "C7")

h.names <- unique(mouse$gs_name)

head(h.names)

#h.sets变成list

h.sets <- vector("list",length=length(h.names))

names(h.sets) <- h.names

for (i in names(h.sets)) {

  h.sets[[i]] <- pull(mouse[mouse$gs_name==i,"gene_symbol"])

}

head(h.sets)

###########################################3.获取表达矩阵

library(singleseqgset)

library(Seurat)

#载入单细胞数据

Myeloid <- readRDS("~/project/LJ.22.02.sc/rdata/Myeloid.rds")

expr.mat=Myeloid@assays$RNA@data

dim(expr.mat)

group=Myeloid$seurat_clusters

table(group)

####计算每组相对于其他组的显著差异基因

logfc.data <- logFC(cluster.ids=Myeloid$seurat_clusters,expr.mat=expr.mat)

head(names(logfc.data))

#提供分组情况，logFC差异倍数

gse.res <- wmw_gsea(expr.mat=expr.mat,cluster.cells=logfc.data[[1]],log.fc.cluster=logfc.data[[2]],gene.sets=h.sets)

#此函数返回两个列表，第一个包含测试统计信息“GSEA_statistics”，第二个包含p值“GSEA_p_values”。

names(gse.res)

res.stats <- gse.res[["GSEA_statistics"]]

res.pvals <- gse.res[["GSEA_p_values"]]

res.stats["ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_1DY_DN",]

results<-merge(res.stats,res.pvals,by="row.names",all.x=TRUE)

colnames(results)[1]<- "C7"

colnames(results)[2:10]<- paste(0:8,"res.stats",by="")

colnames(results)[11:19]<- paste(0:8,"res.pvals",by="")

write_csv(results,"/home/project/LJ.22.02.sc/rdata/C7.results.csv")