GSVA with a genelist on a "per-cluster" basis

Question

[bio_info]

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?

Answer 1

[yura.grabovska]

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.

[bio_info]

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.

[yura.grabovska]

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

[bio_info]

This is really helpful! Thanks a lot, I'll checkout the vignette and try to implement it on my data :)