Question: Visium Data and pulling out cells with high GO enrichment for a given pathway

[DirtyHarry80]

Hi, I am using your package on Visium data as described in your recent manual, it works very fine! Related to that, I am wondering if there is a way to pull out identities of spots that have a particularly high GO enrichment (for a particular pathway or process)? In other words, where are the results of "fgsea" stored in Seurat object?

[vdsukhov]

@DirtyHarry80 Hi

I will consider the code from our vignette for analysis of spatial datasets (https://bioconductor.org/packages/release/bioc/vignettes/fgsea/inst/doc/geseca-tutorial.html#analysis-of-spatial-rna-seq).
In our vignette we first of all perform enrichment analysis and after that to plot the expressions per pathway we use the following code:

topPathways <- gesecaRes[, pathway] |> head(4)

for (ppn in topPathways) {
    pp <- pathways[[ppn]]
    pp <- intersect(pp, rownames(E))
    
    score <- colSums(brain@assays$SCT@scale.data[pp, ])/sqrt(length(pp))
    brain@meta.data[[ppn]] <- score
}

SpatialFeaturePlot(brain, features = topPathways, )

it measures the pathway expression per each spot. After that, the values can be found in the metadata of object (brain@meta.data). For example:

brain@meta.data[, "KEGG_RIBOSOME", drop = FALSE]

Coordinates for each spot could be found in the following attribute - brain@images$anterior1@coordinates.

[Pedramto89]

Hi @vdsukhov
I could generate some results for my data (Visium) but I do not know why it shows only two pathways and for each has two images? It's because of showing only top pathways? How can I visualize the pathways of interest?

[assaron]

@Pedramto89 I've pushed some updates to the plotting code and the vignette. Please, check it out (you need to install the package from github with devtools::install_github('ctlab/fgsea')), may be it will resolve your question. If not, please provide more details: what code you're running and what it produces.

[Pedramto89]

Thank you @assaron
This is the code I use:
`devtools::install_github('ctlab/fgsea')

library(fgsea)
library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
getwd()
data <- readRDS(file = "savedobject.rds")

data <- RunPCA(data, assay = "SCT", verbose = FALSE,
rev.pca = TRUE, reduction.name = "pca.rev", reduction.key="PCR_")
E <- data@reductions$pca.rev@feature.loadings
library(msigdbr)
pathwaysDF <- msigdbr(species = "human", category = "H")
pathways <- split(pathwaysDF$gene_symbol, pathwaysDF$gs_name)

set.seed(1)

gesecaRes <- geseca(pathways, E, minSize = 15, maxSize = 500, center = FALSE)
topPathways <- gesecaRes[, pathway] |> head(10)

for (ppn in topPathways) {
pp <- pathways[[ppn]]
pp <- intersect(pp, rownames(E))

score <- colSums(data@assays$SCT@scale.data[pp, ])/sqrt(length(pp))
data@meta.data[[ppn]] <- score
}
SpatialFeaturePlot(data, features = topPathways, )`
and what I get is the image I uploaded. I changed the head from 4 to 10 and as you see, the legend is so overlapped and I am not sure about the legend scale. It shows what?Also, I still have this question that how can I see the enrichment of a pathway of interest? Thank you!

[assaron]

@Pedramto89 try to use the new plotCoregulationProfileSpatial as in here: https://github.com/ctlab/fgsea/blob/master/vignettes/geseca-tutorial.Rmd The values there are z-scores of the average gene expression of the gene set.

The titles can be shortened for convenience, and you can use cowplot::plot_grid (or other function) to arrange the plots.

topPathways <- gesecaRes[, pathway] |> head(10)
titles <- sub("HALLMARK_", "", topPathways)
ps <- plotCoregulationProfileSpatial(pathways[topPathways], obj,
                                       title=titles)
cowplot::plot_grid(plotlist=ps[1:4], ncol=2)

An individual pathway can be plotted as this:

plotCoregulationProfileSpatial(pathways$HALLMARK_OXIDATIVE_PHOSPHORYLATION,
                               obj,
                               title=sprintf("HALLMARK_OXIDATIVE_PHOSPHORYLATION (pval=%.2g)",
                                             gesecaRes[
                                                 match("HALLMARK_OXIDATIVE_PHOSPHORYLATION", pathway),
                                                 pval]))

[Pedramto89]

Thanks so much for the prompt response @assaron It worked!
Just wanted to let you know that the legend has been removed. Still scaling and interpreting is the question.

[Pedramto89]

Hi
I tried it after a few months and get this weird image: I do not know why some parts of the image is not showing.
Do you think is there any problem with my code?
`library(msigdbr)
library(fgsea)
library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
set.seed(1)
E <- data@reductions$pca.rev@feature.loadings
E <- obj@reductions$pca@feature.loadings
pathwaysDF <- msigdbr(species = "human", category = "H")
pathways <- split(pathwaysDF$gene_symbol, pathwaysDF$gs_name)

gesecaRes <- geseca(pathways, E, minSize = 15, maxSize = 500, center = FALSE)
topPathways <- gesecaRes[, pathway] |> head(10)
titles <- sub("HALLMARK_", "", topPathways)
ps <- plotCoregulationProfileSpatial(pathways[topPathways], obj,
title=titles)
cowplot::plot_grid(plotlist=ps[1:4], ncol=2)

topPathways <- gesecaRes[, pathway] |> head(10)
titles <- sub("HALLMARK_", "", topPathways)
ps <- plotCoregulationProfileSpatial(pathways[topPathways], obj,
title=titles)
cowplot::plot_grid(plotlist=ps[1:4], ncol=2)

plotCoregulationProfileSpatial(pathways$HALLMARK_OXIDATIVE_PHOSPHORYLATION,
obj,
title=sprintf("HALLMARK_OXIDATIVE_PHOSPHORYLATION (pval=%.2g)",
gesecaRes[
match("HALLMARK_OXIDATIVE_PHOSPHORYLATION", pathway),
pval]))
`

[assaron]

@Pedramto89 the images look OK on the first glance. What exactly are you missing?