Hi everyone,
I have a scRNA-seq immune cell dataset (Seurat object) with 8 clusters and I have been asked to perform a GSVA or ssGSEA of a genelist on a "per-cluster" basis. I do not have much experience with enrichment analysis and I tried going through some tutorials but none of them explain how to do this for each cluster.
Could someone please explain how to run a GSVA or ssGSEA for this genelist on all clusters with some steps? Do I need to compute differentially expressed genes for each cluster? What packages can I use for the analysis?
GSVA / ssGSEA were originally designed for bulk expression analysis and so some of the assumptions/caveats with single-cell data don't carry over well. For a very basic approach you could
log.mat <- GetAssayData(seurat.object, assay = "RNA", layer = "data")
And use that to run GSVA/ssGSEA using the GSVA package. IIRC there are some single-cell wrappers for this, but really doesn't extend functionallity very far.
After that you could calculate mean/median scores per cluster using something like
cluster.means <- t(apply(gsva.res, 1, \(x) tapply(x, seurat.object$your.clusters, mean)))
Originally, the advice was for people to use limma on top of GSVA/ssGSEA to do differential expression on the output, but this has caveats for use with single-cell data and I wouldn't advise it. Either way, for an exploratory analysis what you probably want to know more is what expression programmes are active in the clusters rather than have a specific p-value badge to pin on it.
Hi! Thanks a lot for the explanation! How can I integrate my custom genelist in this analysis here? As an input I have my clustered scRNA-seq dataset and a genelist, but I am confused how to carry out a GSVA with these inputs.
library(GSVA)
log.mat <- GetAssayData(seurat.object, assay = "RNA", layer = "data")
genesets <- list(Genelist1 = some.list.of.genes1,
Genelist2 = some.list.of.genes2,
... = ...)
gsva.par <- gsvaParam(log.mat,
genesets,
minSize = 10,
maxSize = 500, ## tailor this to your genelist size
maxDiff = TRUE)
gsva.res <- gsva(gsva.par)
Something like this - I haven't tested it for obvious reasons. See the vignette for GSVA It's pretty self-explanatory
This is really helpful! Thanks a lot, I'll checkout the vignette and try to implement it on my data :)
Log in to answer this question.
