Analysis of 10X visium spatial data
Similarly to single cell RNA-seq, GESECA can be used for gene set enrichment analysis of spatial transcriptomics profiling. As an example we will use glioblastoma sample from Ravi et al (https://pubmed.ncbi.nlm.nih.gov/35700707/).
library(Seurat)
obj <- readRDS(url("https://alserglab.wustl.edu/files/fgsea/275_T_seurat.rds"))
As for scRNA-seq, we apply normalization for 10000 genes and run a do a PCA reduction.
obj <- SCTransform(obj, assay = "Spatial", verbose = FALSE, variable.features.n = 10000)
obj <- RunPCA(obj, assay = "SCT", verbose = FALSE,
rev.pca = TRUE, reduction.name = "pca.rev",
reduction.key="PCR_", npcs = 50)
E <- obj@reductions$pca.rev@feature.loadings
We will use HALLMARK pathways as the gene set collection.
library(msigdbr)
pathwaysDF <- msigdbr("human", collection="H")
pathways <- split(pathwaysDF$gene_symbol, pathwaysDF$gs_name)
Now we can run the analysis (remember, we set center=FALSE because we use the reduced matrix):
set.seed(1)
gesecaRes <- geseca(pathways, E, minSize = 15, maxSize = 500, center = FALSE)
#> Warning in geseca(pathways, E, minSize = 15, maxSize = 500, center = FALSE):
#> For some pathways, in reality P-values are less than 1e-50. You can set the
#> `eps` argument to zero for better estimation.
head(gesecaRes, 10)
#> pathway pctVar pval
#> <char> <num> <num>
#> 1: HALLMARK_HYPOXIA 1.2764439 1.000000e-50
#> 2: HALLMARK_TNFA_SIGNALING_VIA_NFKB 0.5327938 5.698507e-28
#> 3: HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION 0.4649753 1.343028e-24
#> 4: HALLMARK_INTERFERON_GAMMA_RESPONSE 0.4574821 1.520615e-24
#> 5: HALLMARK_INTERFERON_ALPHA_RESPONSE 0.4749296 1.412115e-22
#> 6: HALLMARK_OXIDATIVE_PHOSPHORYLATION 0.2667845 5.482915e-14
#> 7: HALLMARK_MTORC1_SIGNALING 0.2383888 5.287042e-12
#> 8: HALLMARK_GLYCOLYSIS 0.1878601 1.439702e-09
#> 9: HALLMARK_ALLOGRAFT_REJECTION 0.1644770 3.630902e-09
#> 10: HALLMARK_INFLAMMATORY_RESPONSE 0.1670550 7.168505e-09
#> padj log2err size
#> <num> <num> <int>
#> 1: 2.400000e-49 NA 150
#> 2: 6.838209e-27 1.3651796 150
#> 3: 9.123689e-24 1.2790340 145
#> 4: 9.123689e-24 1.2790340 148
#> 5: 6.778151e-22 1.2295041 77
#> 6: 2.193166e-13 0.9545416 181
#> 7: 1.812700e-11 0.8870750 179
#> 8: 4.319105e-09 0.7881868 144
#> 9: 9.682405e-09 0.7614608 98
#> 10: 1.720441e-08 0.7477397 112
Finally, let us plot spatial expression of the top four pathways:
topPathways <- gesecaRes[, pathway] |> head(4)
titles <- sub("HALLMARK_", "", topPathways)
pt.size.factor <- 1.6
# Starting from Seurat version 5.1.0, the scaling method for spatial plots was changed.
# As mentioned here: https://github.com/satijalab/seurat/issues/9049
# This code provides a workaround to adapt to the new scaling behavior.
if (packageVersion("Seurat") >= "5.1.0"){
sfactors <- ScaleFactors(obj@images$slice1)
pt.size.factor <- 1.5 * sfactors$fiducial / sfactors$hires
}
ps <- plotCoregulationProfileSpatial(pathways[topPathways], obj,
title=titles,
pt.size.factor=pt.size.factor)
cowplot::plot_grid(plotlist=ps, ncol=2)
Consistent with the Ravi et al, we see a distinct hypoxic region (defined by HALLMARK_HYPOXIA) and a reactive immune region (defined by HALLMARK_INTERFERON_GAMMA_RESPONSE). Further, we can explore behavior of HALLMARK_OXIDATIVE_PHOSPHORYLATION pathway to see that oxidative metabolism is more characteristic to the “normal” tissue region:
plotCoregulationProfileSpatial(pathways$HALLMARK_OXIDATIVE_PHOSPHORYLATION,
obj,
pt.size.factor=pt.size.factor,
title=sprintf("HALLMARK_OXIDATIVE_PHOSPHORYLATION\n(pval=%.2g)",
gesecaRes[
match("HALLMARK_OXIDATIVE_PHOSPHORYLATION", pathway),
pval]))
#> [[1]]
Analysis of 10X xenium spatial transcriptomics data
GESECA can also be applied to spatial transcriptomics data generated by high-plex in situ technologies with subcellular resolution, such as 10X Genomics’ Xenium platform. In this example, we analyze a Human Ovarian Cancer sample profiled using the Xenium 5K panel.
The raw data used in this analysis is publicly available and can be downloaded from the following sources:
To prepare the data for GESECA analysis, we carried out the following custom preprocessing steps using Seurat:
fldr <- "" # Path to downloaded Xenium output
annf_path <- "" # Path to downloaded cell annotation CSV
# Load the Xenium data as a Seurat object
xobj <- Seurat::LoadXenium(fldr, molecule.coordinates = FALSE)
# Remove assays not needed for downstream analysis
for (nm in names(xobj@assays)[-1]) {
xobj@assays[[nm]] <- NULL
}
# Add spatial coordinates to metadata
coords <- data.table::as.data.table(xobj@images$fov@boundaries$centroids@coords)
xobj@meta.data <- cbind(xobj@meta.data, coords)
# Integrate external cell annotations
annot <- data.table::fread(annf_path)
xobj <- xobj[, rownames(xobj@meta.data) %in% annot$cell_id]
xobj@meta.data$annotation <- annot[match(rownames(xobj@meta.data), annot$cell_id), ]$group
xobj@meta.data$annotation <- factor(xobj@meta.data$annotation)
# Filter out low-quality cells
xobj <- subset(xobj, subset = nCount_Xenium > 50)
# Optional: Crop the region and subsample the object
# These steps reduce size and complexity for a more compact vignette
# Helps keep the vignette lightweight and quick to render
xobj <- xobj[, xobj$x < 3000 & xobj$y < 4000]
set.seed(1)
xobj <- xobj[, sample.int(ncol(xobj), 40000)]
# Normalize and scale data
xobj <- NormalizeData(xobj, normalization.method = "LogNormalize", scale.factor = 1000, verbose = FALSE)
xobj <- FindVariableFeatures(xobj, nfeatures = 2000, verbose = FALSE)
xobj <- ScaleData(xobj, verbose = FALSE)
# Perform PCA and UMAP dimensionality reduction
xobj <- RunPCA(xobj, verbose = FALSE)
xobj <- RunUMAP(xobj, dims = 1:20)
Additionally, before saving the object for use in the vignette, we remove the following components from the Seurat object to reduce its size and keep it tidy:
xobj@reductions$pca <- NULL
xobj@assays$Xenium@layers$scale.data <- NULL
xobj@assays$Xenium@layers$data <- NULL
For this tutorial, we directly use the Seurat object that was preprocessed using the steps described above:
xobj <- readRDS(url("https://alserglab.wustl.edu/files/fgsea/xenium-human-ovarian-cancer.rds"))
xobj <- NormalizeData(xobj, verbose = FALSE)
xobj <- ScaleData(xobj, verbose = FALSE)
The plot below compares the spatial localization of annotated cells (left) with their transcriptional similarity via UMAP (right):
p1 <- ImageDimPlot(xobj, group.by = "annotation",
dark.background = FALSE, size = 0.5,
flip_xy = TRUE) +
theme(legend.position = "none")
p1 <- suppressMessages(p1 + coord_flip() + scale_x_reverse() + theme(aspect.ratio = 1.0))
p2 <- DimPlot(xobj, group.by = "annotation", raster = FALSE, pt.size = 0.5, stroke.size = 0.0) +
coord_fixed()
p1 | p2
The command suppressMessages(p1 + coord_flip() + scale_x_reverse() + theme(aspect.ratio = 1.0)) adjusts the plot orientation to match the spatial layout shown in the Xenium Explorer web viewer, ensuring consistency between the analysis and the original data presentation.
We will now perform the GESECA analysis in the same manner as previously demonstrated with both scRNA-seq and Visium datasets. For consistency and interpretability, we will once again use the HALLMARK gene sets from the MSigDB collection.
xobj <- RunPCA(xobj, verbose = F, rev.pca = TRUE, reduction.name = "pca.rev",
reduction.key="PCR_", npcs=30)
E <- xobj@reductions$pca.rev@feature.loadings
gsets <- msigdbr("human", collection = "H")
gsets <- split(gsets$gene_symbol, gsets$gs_name)
set.seed(1)
gesecaRes <- geseca(gsets, E, minSize = 15, maxSize = 500, center = FALSE, eps = 1e-100)
head(gesecaRes, 10)
#> pathway pctVar pval
#> <char> <num> <num>
#> 1: HALLMARK_E2F_TARGETS 4.8099688 1.010310e-47
#> 2: HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION 4.2252864 6.218037e-45
#> 3: HALLMARK_G2M_CHECKPOINT 4.4643380 6.716647e-44
#> 4: HALLMARK_INTERFERON_GAMMA_RESPONSE 2.0646007 2.455835e-22
#> 5: HALLMARK_INTERFERON_ALPHA_RESPONSE 2.2680902 3.409619e-20
#> 6: HALLMARK_MITOTIC_SPINDLE 2.0814609 6.390584e-19
#> 7: HALLMARK_MTORC1_SIGNALING 1.1236865 4.151984e-12
#> 8: HALLMARK_MYC_TARGETS_V1 1.2533092 2.328815e-11
#> 9: HALLMARK_HYPOXIA 1.0042856 1.929976e-10
#> 10: HALLMARK_TNFA_SIGNALING_VIA_NFKB 0.9767425 4.630348e-10
#> padj log2err size
#> <num> <num> <int>
#> 1: 2.727837e-46 1.7973425 62
#> 2: 8.394350e-44 1.7387813 88
#> 3: 6.044982e-43 1.7208244 59
#> 4: 1.657689e-21 1.2210538 70
#> 5: 1.841194e-19 1.1601796 33
#> 6: 2.875763e-18 1.1239150 26
#> 7: 1.601480e-11 0.8870750 41
#> 8: 7.859752e-11 0.8634154 15
#> 9: 5.789928e-10 0.8266573 58
#> 10: 1.250194e-09 0.8012156 74
Next, we visualize the top enriched pathways identified by GESECA using UMAP embeddings. This allows us to assess transcriptional patterns of gene set activity across cells.
topPathways <- gesecaRes[, pathway] |> head(6)
titles <- sub("HALLMARK_", "", topPathways)
pvals <- gesecaRes[match(topPathways, pathway), pval]
pvals <- paste("p-value:", formatC(pvals, digits = 1, format = "e"))
titles <- paste(titles, pvals, sep = "\n")
plots <- fgsea::plotCoregulationProfileReduction(
gsets[topPathways], xobj,
reduction = "umap", title = titles,
raster = TRUE, raster.dpi = c(500, 500),
pt.size = 2.5
)
plots <- lapply(plots, function(p){
p + theme(plot.title = element_text(hjust=0), text = element_text(size=10))
})
cowplot::plot_grid(plotlist = plots, ncol = 2)
In addition to UMAP plots, you can visualize the spatial expression patterns of the most enriched pathways using the plotCoregulationProfileImage function. This enables a detailed view of how pathway activities vary across spatial coordinates:
imagePlots <- plotCoregulationProfileImage(
gsets[topPathways], object = xobj,
dark.background = FALSE,
title = titles,
size=0.8
)
imagePlots <- lapply(imagePlots, function(p){
suppressMessages(
p + coord_flip() + scale_x_reverse() +
theme(plot.title = element_text(hjust=0), text = element_text(size=10), aspect.ratio = 1.0)
)
})
cowplot::plot_grid(plotlist = imagePlots, ncol = 2)
