As including a more detailed vignette inside the package makes the package exceed the tarball size, more detailed vignettes are hosted on an external website. This is a simplified vignette.

1 Installation

This package can be installed from Bioconductor:

if (!requireNamespace("BiocManager")) install.packages("BiocManager")
BiocManager::install("Voyager")
# Devel version
# install.packages("remotes")
remotes::install_github("pachterlab/Voyager")

2 Introduction

In non-spatial scRNA-seq, the SingleCellExperiment (SCE) package implements a data structure and other packages such as scater implement methods for quality control (QC), basic exploratory data analysis (EDA), and plotting functions, using SCE to organize the data and results. Voyager to SpatialFeatureExperiment (SFE) aims to be analogous scater to SFE, implementing basic exploratory spatial data analysis (ESDA) and plotting. SFE inherits from SCE and SpatialExperiment (SPE), so all methods written for SCE and SPE can be used for SFE as well.

In this first version, ESDA is based on the classic geospatial package spdep, but future versions will incorporate methods from GWmodel, adespatial, and etc.

These are the main functionalities of the Voyager at present:

  • Univariate global spatial statistics, such as Moran’s I, Geary’s C, permutation testing of I and C, correlograms, and global G.
  • Univariate local spatial statistics, such as local Moran’s I, local Geary’s C, Getis-Ord Gi*, Moran scatter plot, and local spatial heteroscedasticity (LOSH).
  • Plotting gene expression and colData along with annotation geometries, with colorblind friendly default palettes. The actual geometries are plotted, not just centroids as in Seurat.
  • Plotting permutation testing results and correlograms, multiple genes in the same plot, can color by gene, sample, or any other attribute.
  • Clustering correlograms and Moran’s scatter plots
  • Plotting local spatial statistics in space
  • Plotting dimension reduction in space
  • Plotting spatial neighborhood graphs

Future versions will add bivariate and multivariate spatial statistics and user friendly wrappers of some successful spatial transcriptomics data analysis packages for spatially variable genes, cell type deconvolution, and spatial regions on CRAN, Bioconductor, pip, and conda, to provide a uniform syntax and avoid object conversion, as is done in Seurat for some non-spatial scRNA-seq methods.

3 Dataset

Here we use a mouse skeletal muscle Visium dataset from Large-scale integration of single-cell transcriptomic data captures transitional progenitor states in mouse skeletal muscle regeneration. It’s in the SFEData package, as an SFE object, which contains Visium spot polygons, myofiber and nuclei segmentations, and myofiber and nuclei morphological metrics.

library(SFEData)
library(SpatialFeatureExperiment)
library(Voyager)
library(scater)
#> Loading required package: SingleCellExperiment
#> Loading required package: SummarizedExperiment
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#> 
#> Attaching package: 'MatrixGenerics'
#> The following objects are masked from 'package:matrixStats':
#> 
#>     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#>     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#>     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#>     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#>     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#>     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#>     colWeightedMeans, colWeightedMedians, colWeightedSds,
#>     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#>     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#>     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#>     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#>     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#>     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#>     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#>     rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#> 
#>     Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
#>     as.data.frame, basename, cbind, colnames, dirname, do.call,
#>     duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
#>     lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
#>     pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
#>     tapply, union, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> 
#> Attaching package: 'S4Vectors'
#> The following objects are masked from 'package:base':
#> 
#>     I, expand.grid, unname
#> Loading required package: IRanges
#> Loading required package: GenomeInfoDb
#> Loading required package: Biobase
#> Welcome to Bioconductor
#> 
#>     Vignettes contain introductory material; view with
#>     'browseVignettes()'. To cite Bioconductor, see
#>     'citation("Biobase")', and for packages 'citation("pkgname")'.
#> 
#> Attaching package: 'Biobase'
#> The following object is masked from 'package:MatrixGenerics':
#> 
#>     rowMedians
#> The following objects are masked from 'package:matrixStats':
#> 
#>     anyMissing, rowMedians
#> Loading required package: scuttle
#> Loading required package: ggplot2
sfe <- McKellarMuscleData()
#> snapshotDate(): 2022-10-31
#> see ?SFEData and browseVignettes('SFEData') for documentation
#> loading from cache
# Only use spots in tissue here
sfe <- sfe[,colData(sfe)$in_tissue]
sfe <- logNormCounts(sfe)
sfe
#> class: SpatialFeatureExperiment 
#> dim: 15123 932 
#> metadata(0):
#> assays(2): counts logcounts
#> rownames(15123): ENSMUSG00000025902 ENSMUSG00000096126 ...
#>   ENSMUSG00000064368 ENSMUSG00000064370
#> rowData names(6): Ensembl symbol ... vars cv2
#> colnames(932): AAACATTTCCCGGATT AAACCTAAGCAGCCGG ... TTGTGTTTCCCGAAAG
#>   TTGTTGTGTGTCAAGA
#> colData names(13): barcode col ... in_tissue sizeFactor
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> spatialCoords names(2) : imageX imageY
#> imgData names(1): sample_id
#> 
#> Geometries:
#> colGeometries: spotPoly (POLYGON) 
#> annotGeometries: tissueBoundary (POLYGON), myofiber_full (POLYGON), myofiber_simplified (POLYGON), nuclei (POLYGON), nuclei_centroid (POINT) 
#> 
#> Graphs:
#> Vis5A:

4 Univariate spatial statistics

A spatial neighborhood graph is required for all spdep analyses.

colGraph(sfe, "visium") <- findVisiumGraph(sfe)

All of the numerous univariate methods can be used with runUnivariate(), which stores global results in rowData(sfe) and local results in localResults(sfe). Here we compute Moran’s I for one gene. While Ensembl IDs are used internally, the user can specify more human readable gene symbols. A warning will be given if the gene symbol matches multiple Ensembl IDs.

features_use <- c("Myh1", "Myh2")
sfe <- runUnivariate(sfe, type = "moran", features = features_use, 
                     colGraphName = "visium")
# Look at the results
rowData(sfe)[rowData(sfe)$symbol %in% features_use,]
#> DataFrame with 2 rows and 8 columns
#>                               Ensembl      symbol            type     means
#>                           <character> <character>     <character> <numeric>
#> ENSMUSG00000033196 ENSMUSG00000033196        Myh2 Gene Expression   0.97476
#> ENSMUSG00000056328 ENSMUSG00000056328        Myh1 Gene Expression   4.82572
#>                         vars       cv2 moran_Vis5A   K_Vis5A
#>                    <numeric> <numeric>   <numeric> <numeric>
#> ENSMUSG00000033196   24.0374   25.2984    0.625500   2.21641
#> ENSMUSG00000056328  302.2385   12.9785    0.635718   2.68736

Since Moran’s I is very commonly used, one can call runMoransI rather than runUnivariate.

Compute a local spatial statistic, Getis-Ord Gi*, which is commonly used to detect hotspots and coldspots. The include_self argument is only for Getis-Ord Gi*; when set to TRUE Gi* is computed as the spatial graph includes self-directing edges, and otherwise Gi is computed.

sfe <- runUnivariate(sfe, type = "localG", features = features_use,
                     colGraphName = "visium", include_self = TRUE)
# Look at the results
DataFrame(localResults(sfe, name = "localG"))
#> DataFrame with 932 rows and 2 columns
#>     ENSMUSG00000056328 ENSMUSG00000033196
#>              <numeric>          <numeric>
#> 1             0.953965           3.526243
#> 2             2.219643           4.718390
#> 3            -1.450029          -1.086702
#> 4             3.545740           5.322070
#> 5             1.194738           0.282678
#> ...                ...                ...
#> 928           1.635331           0.204077
#> 929           1.821201           5.413792
#> 930          -2.799465          -2.537571
#> 931          -2.708677          -1.905788
#> 932           0.982869          -1.231641

Spatial statistics can also be computed for numeric columns of colData(sfe), with colDataUnivariate(), and for numeric attributes of the geometries with colGeometryUnivariate() and annotGeometryUnivariate(), all with very similar arguments.

5 Plotting

Plot gene expression and colData(sfe) together with annotation geometry. Here nCounts is the total UMI counts per spot, which is in colData.

plotSpatialFeature(sfe, c("nCounts", "Myh1"), colGeometryName = "spotPoly", 
                   annotGeometryName = "myofiber_simplified", 
                   aes_use = "color", size = 0.4, fill = NA,
                   annot_aes = list(fill = "area"))