--- title: "GSVA on spatial omics data" author: - name: Robert Castelo affiliation: - &idupf Dept. of Medicine and Life Sciences, Universitat Pompeu Fabra, Barcelona, Spain email: robert.castelo@upf.edu - name: Pere Moles affiliation: *idupf email: pere.moles@upf.edu - name: Justin Guinney affiliation: - Tempus Labs, Inc. email: justin.guinney@tempus.com abstract: > Here we illustrate how to use GSVA with spatial omics data. date: "`r BiocStyle::doc_date()`" package: "`r pkg_ver('GSVA')`" vignette: > %\VignetteIndexEntry{GSVA on spatial omics data} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} %\VignetteKeywords{GeneExpression, Microarray, RNAseq, GeneSetEnrichment, Pathway} output: BiocStyle::html_document: toc: true toc_float: true number_sections: true fig_captions: yes bibliography: GSVA.bib --- **License**: `r packageDescription("GSVA")[["License"]]` ```{r setup, include=FALSE} options(width=80) knitr::opts_chunk$set(collapse=TRUE, message=FALSE, warning=FALSE, comment="", fig.align="center", fig.wide=TRUE) ``` # Introduction GSVA provides now specific support for spatial omics data in the algorithm that runs through the `gsvaParam()` parameter constructor, and originally described in the publication by @haenzelmann2013gsva. At the moment, this specific support consists of the following features: * The currently available container for spatial omics data that allows one to input additional row and column metadata, and spatial coordinates, is a `SpatialExperiment` object. * While the input spatial omics data is always sparse, the output of enrichment scores will be always dense, and therefore, the container storing those scores will be different from the input data, typically a `matrix` or a dense `DelayedMatrix` object using an `HDF5Matrix` backend. The latter will be particularly used when the total number of values exceeds 2^31, which is the largest 32-bit standard integer value in R. * By default, when the input expression data is stored in a sparse data container, as it typically happens with spatial omics data, then a sparse regime of the GSVA algorithm will run, if GSVA is the chosen method, by which nonzero values are treated differently from zero values, leading to slightly different results than those obtained by applying the classical GSVA algorithm. If we set the parameter `sparse=FALSE` in the call to `gsvaParam()`, the classical GSVA algorithm will be used, which for a typical spatial omics data set will result in longer running times and larger memory consumption than running it in the default sparse regime for this type of data. * The GSVA algorithm can be run either at once through a call to `gsva()` with a parameter object or in two steps: (1) row normalization and column rank transformation with `gsvaRanks()`; and (2) column enrichment scores calculation with `gsvaScores()`. Splitting the GSVA algorithm into these two steps allows one to reuse the output of the first step, which is independent of the gene sets, to calculate enrichment scores for different collections of gene sets, without having to repeat the first step. In what follows, we will illustrate the use of GSVA on a publicly available spatial transcriptomics transcriptomics data set published by @maynard2021transcriptome, which is available through the `r Biocpkg("spatialLIBD")` package, and which contains 12 samples of human dorsolateral prefrontal cortex (DLPFC) tissue. # Import data We import the DLPFC data using the `r Biocpkg("spatialLIBD")` package, as a `r Biocpkg("SpatialExperiment")` object, following the instructions at http://research.libd.org/spatialLIBD/articles/spatialLIBD.html. ```{r, message=FALSE, warning=FALSE} library(SpatialExperiment) library(ExperimentHub) library(spatialLIBD) ehub <- ExperimentHub() spe <- fetch_data(type="spe", eh=ehub) ``` Note that this data set contains 12 samples. ```{r} table(spe$sample_id) ``` For the purpose of speeding up processing this vignette, we will analyze here only the sample with identifier `151673`. ```{r} spe <- spe[, spe$sample_id == "151673"] dim(spe) ``` There is an indicator column called `in_tissue` for identifying spots overlapping the tissue, but apparently all spots have such a condition. ```{r} stopifnot(all(colData(spe)$in_tissue)) ``` The column data includes a column called `spatialLIBD` with anatomically annotated groups of spots (six cortical layers from the grey matter and the white matter). This annotation is apparently missing for a small number of spots, we are going to discard those spots. ```{r} table(colData(spe)$spatialLIBD, useNA="always") spe <- spe[, !is.na(colData(spe)$spatialLIBD)] dim(spe) table(colData(spe)$spatialLIBD, useNA="always") ``` # Quality control and pre-processing Here, we perform a quality control (QC) and pre-processing steps using the package `r Biocpkg("scrapper")` [@lun2022powering]. We start identifying mitochondrial genes. ```{r, message=FALSE, warning=FALSE} library(scrapper) is_mito <- grepl("^MT-", rowData(spe)$gene_name) table(is_mito) ``` Calculate QC metrics and tally the number of low-quality spots. ```{r} spe <- quickRnaQc.se(spe, subsets=list(mito=is_mito)) spe$discard <- !spe$keep table(spe$discard) ``` Plot detected low-quality spots to ensure they are not located in a specific region of the tissue. ```{r dlpfcqcplot, echo=TRUE, message=FALSE, fig.height=5, fig.width=5, out.width="600px", fig.cap="Detected low-quality spots in the DLPFC sample 151673."} library(ggspavis) plotObsQC(spe, plot_type="spot", annotate="discard") ``` Filter out low-quality spots. ```{r} spe <- spe[, spe$keep] dim(spe) ``` Filter out genes that are expressed in less than 1% of the spots. ```{r} spotsxgene <- rowSums(counts(spe) > 0) spe <- spe[spotsxgene > floor(ncol(spe)*0.01), ] dim(spe) ``` Calculate library size factors and normalized units of expression in logarithmic scale. ```{r} spe <- normalizeRnaCounts.se(spe, size.factors=spe$sum) assayNames(spe) ``` # Spot cell-type annotation using GSVA Here we use the a collection of brain cell-type marker genes derived from the single-nuclei RNA-seq (snRNA-seq) data set published by @tran2021single, to estimate the most abundant cell type in each spot of the DLPFC sample 151673. This collection of marker gene sets is available in the `r Biocpkg("GSVAdata")` package, which contains the script that generated the marker gene sets, and can be imported as follows. ```{r, message=FALSE, warning=FALSE} library(GSEABase) library(GSVA) fname <- file.path(system.file("extdata", package="GSVAdata"), "human_brain_snRNAseq_cellType_markers.gmt.gz") gsets <- readGMT(fname) gsets ``` We first build a parameter object using the function `gsvaParam()`. By default, the expression values in the `logcounts` assay will be selected for downstream analysis. ```{r} gsvapar <- gsvaParam(spe, gsets) gsvapar ``` Second, we call the `gsva()` function to calculate the GSVA scores for each spot and each gene set. The output will be a `SpatialExperiment` object with the GSVA scores stored in a new assay called `es`. When the input spatial omics data is very large, you can use the two-step approach described in the single-cell RNA-seq vignette. ```{r} es <- gsva(gsvapar) es ``` Finally, we assign to each spot the gene set with highest GSVA score, and store that assignment as a new metadata column into the `spe` object. ```{r} whmax <- apply(assay(es), 2, which.max) es$cellType <- factor(rownames(es)[whmax]) ``` Figure \@ref(fig:dlpfcgsva) shows the annotated layers of the DLPFC data set and the spatial distribution of the most likely cell type in each spot, based on the GSVA scores. ```{r dlpfcgsva, echo=TRUE, message=FALSE, fig.height=10, fig.width=7, out.width="800px", dpi=100, fig.cap="(a) Annotated layers of the DLPFC sample 151673 and (b) spatial distribution of the most likely cell type in each spot, based on GSVA scores."} library(ggspavis) library(patchwork) library(RColorBrewer) ctpal <- colorRampPalette(brewer.pal(9, "Set1"))(nlevels(es$cellType)) plts <- list(plotVisium(spe, point_size=0.60, facets=NULL, annotate="spatialLIBD", pal=brewer.pal(nlevels(spe$spatialLIBD), "Set1")) + labs(tag="a"), plotVisium(es, assay="es", point_size=0.60, facets=NULL, annotate="cellType", pal=ctpal) + labs(tag="b") ) wrap_plots(plts, nrow=1) ``` We can see that the white matter layer is enriched in oligodendrocytes, while the gray matter layers are enriched in excitatory and inhibitory neurons, and astrocytes. While in this case it is clear that oligodendrocytes are the most abundant cell type in the white matter layer, we can also verify this observation by calculating the Moran's I spatial autocorrelation statistic for the GSVA scores of every gene set and rank them, using the function `spatCor()` implemented also in the GSVA package, which can take as input the `SpatialExperiment` object with the GSVA scores obtained in the previous step. ```{r} I <- spatCor(es) I[order(I$observed, decreasing=TRUE), ] ``` We can observe that the `Oligo` gene set has the highest Moran's I statistic, which follows from the spatial distribution of the GSVA scores for that gene set, which is shown in Figure \@ref(fig:dlpfcgsvaoligo). ```{r dlpfcgsvaoligo, echo=TRUE, message=FALSE, fig.height=10, fig.width=7, out.width="800px", dpi=100, fig.cap="Spatial distribution of the GSVA scores for the Oligodendrocytes gene set in the DLPFC sample 151673."} plotVisium(es, assay="es", point_size=1, facets=NULL, annotate="Oligo") ``` # Benchmarking We are still benchmarking and testing this version of GSVA for spatial omics data. If you encounter problems or have suggestions, do not hesitate to contact us by opening an [issue](https://github.com/rcastelo/GSVA/issues) in the GSVA GitHub repo. # Session information {.unnumbered} Here is the output of `sessionInfo()` on the system on which this document was compiled running pandoc `r rmarkdown::pandoc_version()`: ```{r session_info, cache=FALSE} sessionInfo() ``` # References