chuiqin / irGSEA

The integration of single cell rank-based gene set enrichment analysis
Other
110 stars 17 forks source link

How to calculate a diy geneset? #5

Closed htqqdd closed 2 years ago

htqqdd commented 2 years ago

Thanks for building such an amazing package. I would like to know how to input my diy gene set for evaluation? Can you provide a data type requirement or an example for a custom gene set? Thank you!

chuiqin commented 2 years ago

If you have your own geneset:

# load PBMC dataset by R package SeuratData
library(Seurat)
library(SeuratData)
# download 3k PBMCs from 10X Genomics
InstallData("pbmc3k")
data("pbmc3k.final")
pbmc3k.final <- SeuratObject::UpdateSeuratObject(pbmc3k.final)

# custom geneset
markers <- list()
markers$Tcell_gd <- c("TRDC+", "TRGC1+", "TRGC2+", "TRDV1+","TRAC-","TRBC1-","TRBC2-")
markers$Tcell_NK <- c("FGFBP2+", "SPON2+", "KLRF1+", "FCGR3A+", "CD3E-","CD3G-")
markers$Tcell_CD4 <- c("CD4","CD40LG")
markers$Tcell_CD8 <- c("CD8A","CD8B")
markers$Tcell_Treg <- c("FOXP3","IL2RA")
pbmc3k.final3 <- irGSEA.score(object = pbmc3k.final, assay = "RNA", slot = "data",
custom = T, geneset = markers, method = c("AUCell", "UCell", "singscore", "ssgsea"),
kcdf = 'Gaussian')

If you download your geneset from MiSigdb (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp):

# read geneset and convert it to a list
library(clusterProfiler)
h.geneset <-  read.gmt("~/h.all.v7.2.symbols.gmt")
h.geneset <- read.gmt("~/h.all.v7.2.symbols.gmt") %>% 
group_split(term,.keep = F) %>% 
set_names(levels(factor(h.geneset$term))) %>% 
map(~.x %>% pull())
pbmc3k.final3 <- irGSEA.score(object = pbmc3k.final, assay = "RNA", slot = "data",
custom = T, geneset = h.geneset, method = c("AUCell", "UCell", "singscore", "ssgsea"),
kcdf = 'Gaussian')