Exploratory spatial data analysis with Voyager

Column geometries

Column geometries or colGeometries are the geometries that correspond to columns of the gene count matrix, such as Visium spots and cells in datasets from a single cell resolution technology. Each SFE object can have multiple column geometries. For example, in a dataset with single cell resolution, whole cell segmentation and nuclei segmentation are two different colGeometries. However, for Visium, the spot polygons are the only colGeometry obviously relevant, though users can add other geometries such as results of geometric operations on the spot polygons. The different geometries can be get or set with their names, and “spotPoly” is the standard name for Visium spot polygons.

# Get Visium spot polygons
(spots <- colGeometry(sfe, "spotPoly"))
#> Simple feature collection with 77 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 5000 ymin: 13000 xmax: 7000 ymax: 15000
#> CRS:           NA
#> First 10 features:
#>                                        geometry sample_id
#> AAATTACCTATCGATG POLYGON ((6472.186 13875.23...     Vis5A
#> AACATATCAACTGGTG POLYGON ((5778.291 13635.43...     Vis5A
#> AAGATTGGCGGAACGT POLYGON ((7000 13809.84, 69...     Vis5A
#> AAGGGACAGATTCTGT POLYGON ((6749.535 13874.64...     Vis5A
#> AATATCGAGGGTTCTC POLYGON ((5500.941 13636.03...     Vis5A
#> AATGATGATACGCTAT POLYGON ((6612.42 14598.82,...     Vis5A
#> AATGATGCGACTCCTG POLYGON ((5501.981 14118.62...     Vis5A
#> AATTCATAAGGGATCT POLYGON ((6889.769 14598.22...     Vis5A
#> ACGAGTACGGATGCCC POLYGON ((5084.397 13395.63...     Vis5A
#> ACGCTAGTGATACACT POLYGON ((5639.096 13394.44...     Vis5A

plot(st_geometry(spots))

# Set colGeometry
colGeometry(sfe, "spotPoly") <- spots

To see which colGeometries are present in the SFE object:

colGeometryNames(sfe)
#> [1] "spotPoly"

There are shorthands for some specific column or row geometries. For example, spotPoly(sfe) is equivalent to colGeometry(sfe, "spotPoly") shown above.

# Getter
(spots <- spotPoly(sfe))
#> Simple feature collection with 77 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 5000 ymin: 13000 xmax: 7000 ymax: 15000
#> CRS:           NA
#> First 10 features:
#>                                        geometry sample_id
#> AAATTACCTATCGATG POLYGON ((6472.186 13875.23...     Vis5A
#> AACATATCAACTGGTG POLYGON ((5778.291 13635.43...     Vis5A
#> AAGATTGGCGGAACGT POLYGON ((7000 13809.84, 69...     Vis5A
#> AAGGGACAGATTCTGT POLYGON ((6749.535 13874.64...     Vis5A
#> AATATCGAGGGTTCTC POLYGON ((5500.941 13636.03...     Vis5A
#> AATGATGATACGCTAT POLYGON ((6612.42 14598.82,...     Vis5A
#> AATGATGCGACTCCTG POLYGON ((5501.981 14118.62...     Vis5A
#> AATTCATAAGGGATCT POLYGON ((6889.769 14598.22...     Vis5A
#> ACGAGTACGGATGCCC POLYGON ((5084.397 13395.63...     Vis5A
#> ACGCTAGTGATACACT POLYGON ((5639.096 13394.44...     Vis5A

# Setter
spotPoly(sfe) <- spots

Annotation

Annotation geometries can be get or set with annotGeometry(). In column or row geometries, the number of rows of the sf data frame (i.e. the number of geometries in the data frame) is constrained by the number of rows or columns of the gene count matrix respectively, because just like rowData and colData, each row of a rowGeometry or colGeometry sf data frame must correspond to a row or column of the gene count matrix respectively. In contrast, an annotGeometry sf data frame can have any dimension, not constrained by the dimension of the gene count matrix.

# Getter, by name or index
(tb <- annotGeometry(sfe, "tissueBoundary"))
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 5094 ymin: 13000 xmax: 7000 ymax: 14969
#> CRS:           NA
#>   ID                       geometry sample_id
#> 7  7 POLYGON ((5094 13000, 5095 ...     Vis5A

plot(st_geometry(tb))

# Setter, by name or index
annotGeometry(sfe, "tissueBoundary") <- tb

See which annotGeometries are present in the SFE object:

annotGeometryNames(sfe)
#> [1] "tissueBoundary"      "myofiber_full"       "myofiber_simplified"
#> [4] "nuclei"              "nuclei_centroid"

There are shorthands for specific annotation geometries. For example, tissueBoundary(sfe) is equivalent to annotGeometry(sfe, "tissueBoundary"). cellSeg() (cell segmentation) and nucSeg() (nuclei segmentation) would first query colGeometries (for single cell, single molecule technologies, equivalent to colGeometry(sfe, "cellSeg") or colGeometry(sfe, "nucSeg")), and if not found, they will query annotGeometries (for array capture and microdissection technologies, equivalent to annotGeometry(sfe, "cellSeg") or annotGeometry(sfe, "nucSeg")).

# Getter
(tb <- tissueBoundary(sfe))
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 5094 ymin: 13000 xmax: 7000 ymax: 14969
#> CRS:           NA
#>   ID                       geometry sample_id
#> 7  7 POLYGON ((5094 13000, 5095 ...     Vis5A

# Setter
tissueBoundary(sfe) <- tb

Spatial graphs

The spatial neighborhood graphs for Visium spots are stored in the colGraphs field, which has similar user interface as colGeometries. SFE also wraps all methods to find the spatial neighborhood graph implemented in the spdep package, and triangulation is used here as demonstration.

(g <- findSpatialNeighbors(sfe, MARGIN = 2, method = "tri2nb"))
#> Characteristics of weights list object:
#> Neighbour list object:
#> Number of regions: 77 
#> Number of nonzero links: 428 
#> Percentage nonzero weights: 7.218755 
#> Average number of links: 5.558442 
#> 
#> Weights style: W 
#> Weights constants summary:
#>    n   nn S0      S1       S2
#> W 77 5929 77 28.0096 309.4083

plot(g, coords = spatialCoords(sfe))

# Set graph by name
colGraph(sfe, "graph1") <- g

# Get graph by name
(g <- colGraph(sfe, "graph1"))
#> Characteristics of weights list object:
#> Neighbour list object:
#> Number of regions: 77 
#> Number of nonzero links: 428 
#> Percentage nonzero weights: 7.218755 
#> Average number of links: 5.558442 
#> 
#> Weights style: W 
#> Weights constants summary:
#>    n   nn S0      S1       S2
#> W 77 5929 77 28.0096 309.4083

For Visium, spatial neighborhood graph of the hexagonal grid can be found with the known locations of the barcodes. One SFE object can have multiple colGraphs.

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

plot(colGraph(sfe, "visium"), coords = spatialCoords(sfe))

Which graphs are present in this SFE object?

colGraphNames(sfe)
#> [1] "graph1" "visium"

While this workshop only works with one sample, i.e. tissue section, operations on multiple samples is discussed in the vignette of the SFE package.

Spots in tissue

While the example dataset has all Visium spots whether on tissue or not, only spots that intersect tissue are used for further analyses.

names(colData(sfe))
#>  [1] "barcode"   "col"       "row"       "x"         "y"         "dia"      
#>  [7] "tissue"    "sample_id" "nCounts"   "nGenes"    "prop_mito" "in_tissue"

Total UMI counts (nCounts), number of genes detected per spot (nGenes), and the proportion of mitochondrially encoded counts (prop_mito) have been precomputed and are in colData(sfe). The plotSpatialFeature() function can be used to visualize various attributes in space: the expression of any gene, colData values, and geometry attributes in colGeometry and annotGeometry. The Visium spots are plotted as polygons reflecting their actual size relative to the tissue, rather than as points, as is the case in other packages that plot Visium data. The plotting of geometries is being performed under the hood with geom_sf.

The tissue boundary was found by thresholding the H&E image and removing small polygons that are most likely debris. The in_tissue column of colData(sfe) indicates which Visium spot polygon intersects the tissue polygon; this can be found with SpatialFeatureExperiment::annotPred().

We demonstrate the use of scran (Lun, Bach, and Marioni 2016) for normalization below, although we note that it is not necessarily the best approach to normalizing spatial transcriptomics data. The problem of when and how to normalize spatial transcriptomics data is non-trivial because, as the nCounts plot in space shows above, spatial autocorrelation is evident. Furthermore, in Visium, reverse transcription occurs in situ on the spots, but PCR amplification occurs after the cDNA is dissociated from the spots. Artifacts may be subsequently introduced from the amplification step, and these would not be associated with spatial origin. Spatial artifacts may arise from the diffusion of transcripts and tissue permeablization. However, given how the total counts seem to correspond to histological regions, the total counts may have a biological component and hence should not be treated as a technical artifact to be normalized away as in scRNA-seq data normalization methods. In other words, the issue of normalization for spatial transcriptomics data, and Visium in particular, is complex and is currently unsolved.

sfe_tissue <- sfe[,colData(sfe)$in_tissue]
sfe_tissue <- sfe_tissue[rowSums(counts(sfe_tissue)) > 0,]

#clusters <- quickCluster(sfe_tissue)
#sfe_tissue <- computeSumFactors(sfe_tissue, clusters=clusters)
#sfe_tissue <- sfe_tissue[, sizeFactors(sfe_tissue) > 0]
sfe_tissue <- logNormCounts(sfe_tissue)

Myofiber and nuclei segmentation polygons are available in this dataset in the annotGeometries field. Myofibers were manually segmented, and nuclei were segmented with StarDist trained with a manually segmented subset.

annotGeometryNames(sfe_tissue)
#> [1] "tissueBoundary"      "myofiber_full"       "myofiber_simplified"
#> [4] "nuclei"              "nuclei_centroid"

plotSpatialFeature(sfe_tissue, features = "nCounts", 
                   colGeometryName = "spotPoly",
                   annotGeometryName = "myofiber_simplified", 
                   aes_use = "color", linewidth = 0.5, fill = NA,
                   annot_aes = list(fill = "area"))

colData(sfe_tissue)$n_myofibers <- 
  annotNPred(sfe_tissue, colGeometryName = "spotPoly",
             annotGeometryName = "myofiber_simplified")

plotSpatialFeature(sfe_tissue, features = "n_myofibers", 
                   colGeometryName = "spotPoly", image = "lowres", color = "black",
                   linewidth = 0.1)

colData(sfe_tissue)$mean_myofiber_area <- 
  annotSummary(sfe_tissue, "spotPoly", "myofiber_simplified", 
               annotColNames = "area")[,1] # it always returns a data frame

# The gray spots don't intersect any myofiber
plotSpatialFeature(sfe_tissue, "mean_myofiber_area", "spotPoly", image = "lowres", 
                   color = "black", linewidth = 0.1)

plotColData(sfe_tissue, x = "nCounts", y = "nGenes", colour_by = "mean_myofiber_area")

plotColData(sfe_tissue, x = "nCounts", y = "prop_mito", colour_by = "mean_myofiber_area")

markers <- c(I = "Myh7", IIa = "Myh2", IIb = "Myh4", IIx = "Myh1")

plotSpatialFeature(sfe_tissue, "Myh2", "spotPoly",
                   annotGeometryName = "myofiber_simplified", 
                   annot_aes = list(fill = "area"), swap_rownames = "symbol", 
                   exprs_values = "logcounts", aes_use = "color", linewidth = 0.5,
                   fill = NA)

Moran’s I

As a reference, Pearson correlation is defined as

There are several ways to quantify spatial autocorrelation, the most common of which is Moran’s I (Moran 1950):

where $n$ is the number of spots or locations, $i$ and $j$ are different locations, or spots in the Visium context, $x$ and $y$ are variables with values at each location, and $w_{ij}$ is a spatial weight, which can be inversely proportional to distance between spots or an indicator of whether two spots are neighbors, subject to various definitions of neighborhood and whether to normalize the number of neighbors. The spdep package uses the neighborhood.

Moran’s I is analogous to the Pearson correlation between the value at each location and the average value at its neighbors (but not identical, see (Lee 2001)). Just like Pearson correlation, Moran’s I is generally bound between -1 and 1, where positive value indicates positive spatial autocorrelation, i.e. nearby values tend to be more similar, and negative value indicates negative spatial autocorrelation, i.e. nearby values tend to be more dissimilar.

Upon visual inspection earlier in the workshop, total UMI counts per spot (nCounts) seem to have spatial autocorrelation. For numeric columns of colData(sfe), all univariate methods can be called with colDataUnivariate(). Here we compute Moran’s I for nCounts and nGenes:

sfe_tissue <- colDataUnivariate(sfe_tissue, type = "moran", 
                                features = c("nCounts", "nGenes"),
                                colGraphName = "visium")
colFeatureData(sfe_tissue)[c("nCounts", "nGenes"),]
#> DataFrame with 2 rows and 2 columns
#>         moran_Vis5A   K_Vis5A
#>           <numeric> <numeric>
#> nCounts    0.528705   3.00082
#> nGenes     0.384028   3.88036

For colData, the results are added to colFeatureData(sfe), and features for which Moran’s I is not calculated have NA. The column names of featureData distinguishes between different samples (there’s only one sample in this dataset), and are parsed by plotting functions. Here the first column is the Moran’s I value, which indicates moderate positive spatial autocorrelation for both nCounts and nGenes. The second column is kurtosis of the data.

Compute Moran’s I for attributes of geometries: Here “area” is the area of the cross section of each myofiber as seen in this tissue section and “eccentricity” is the eccentricity of the ellipse fitted to each myofiber.

# Remember zero.policy = TRUE since there're singletons
sfe_tissue <- annotGeometryUnivariate(sfe_tissue, type = "moran",
                                      features = c("area", "eccentricity"), 
                                      annotGeometryName = "myofiber_simplified",
                                      annotGraphName = "myofiber_poly2nb", 
                                      zero.policy = TRUE)
attr(annotGeometry(sfe_tissue, "myofiber_simplified"), "featureData")[c("area", "eccentricity"),]
#> DataFrame with 2 rows and 2 columns
#>              moran_Vis5A   K_Vis5A
#>                <numeric> <numeric>
#> area            0.327888   4.95675
#> eccentricity    0.110938   3.26913

For a non-geometry column in a colGeometry, colGeometryUnivariate() is like annotGeometryUnivariate() here, but none of the colGeometries in this dataset has extra columns.

For gene expression, the logcounts assay is used by default (use the exprs_values argument to change the assay), though this may or may not be best practice. If the metrics are computed for a large number of features, parallel computing is supported, with BiocParallel, with the BPPARAM argument.

sfe_tissue <- runUnivariate(sfe_tissue, type = "moran", features = hvgs, 
                            colGraphName = "visium", BPPARAM = SerialParam())
rowData(sfe_tissue)[head(hvgs),c("moran_Vis5A", "K_Vis5A", "symbol")]
#> DataFrame with 6 rows and 3 columns
#>                    moran_Vis5A   K_Vis5A      symbol
#>                      <numeric> <numeric> <character>
#> ENSMUSG00000029304    0.734937   1.63516        Spp1
#> ENSMUSG00000050708    0.665563   1.81841        Ftl1
#> ENSMUSG00000050335    0.741474   1.68098      Lgals3
#> ENSMUSG00000021939    0.708362   1.86896        Ctsb
#> ENSMUSG00000021190    0.659916   1.66838        Lgmn
#> ENSMUSG00000018893    0.675840   1.82510          Mb

As Moran’s I is very commonly used, runMoransI(sfe_tissue, features = hvgs) is equivalent to runUnivariate(sfe_tissue, type = "moran", features = hvgs).

Exercises

1. Use listSFEMethods() to find the “name” of Geary’s C (Geary 1954). This name should be used in the type argument in runUnivariate().
2. Compute Geary’s C on the highly variable genes, and show the results. Interpretation of Geary’s C: a value below 1 indicates positive spatial autocorrelation, while a value above 1 indicates negative spatial autocorrelation.

Further reading

1. Spatial transcriptomics data is usually much larger than the typical geospatial dataset back in the 1950s when Moran’s I and Geary’s C were devised. See (Luo, Griffith, and Wu 2019) for asymptotic properties of Moran’s I for large datasets with normal and skewed distributions.
2. The negative binomial distribution is often used to model transcriptomics data due to bursts in transcription, although the Poisson distribution is sometimes used instead to simplify the math. See (Griffith and Haining 2006) for a consideration of the Poisson distribution in spatial analyses.
3. Moran’s I is not exactly the same as Pearson correlation between the values themselves and the spatially smoothed values. The bounds of Moran’s I depend on the spatial neighborhood graph. Usually the upper bound is around 1, while the lower bound is closer to -0.5 than -1. See (Jong, Sprenger, and Veen 1984) for a derivation of extreme values of Moran’s I and Geary’s C.
4. Spatial autocorrelation decays at different length scales for different features, and the correlogram is one way to find the length scales. These are the vignettes that use correlograms. Also see this vignette on Moran’s I flipping signs at different length scales.

Local Moran’s I

To recap, global Moran’s I is defined as

Local Moran’s I (Anselin 1995) is defined as

It’s similar to global Moran’s I, but the values at locations $i$ are not summed and there’s no normalization by the sum of spatial weights. Local Moran’s I has been used in spatial transcriptomics in the MERINGUE package (Miller et al. 2021). Here we compute local Moran’s I for the gene Myh2.

sfe_tissue <- runUnivariate(sfe_tissue, type = "localmoran", features = "Myh2",
                            colGraphName = "visium", swap_rownames = "symbol")

It is useful to plot the log normalized Myh2 gene expression as context to interpret the local results:

plotSpatialFeature(sfe_tissue, features = "Myh2", colGeometryName = "spotPoly",
                   swap_rownames = "symbol", image_id = "lowres", color = "black",
                   linewidth = 0.1)

Any local spatial results can be plotted with plotLocalResults(), which is similar to plotSpatialFeature(). Here a divergent palette is used because Moran’s I has a sensible center at 0 (actually the expected value of Moran’s I is -1/(n-1) when the mean is unknown, but it’s very close to 0 as n is typically large in spatial -omics).

plotLocalResult(sfe_tissue, "localmoran", features = "Myh2", 
                colGeometryName = "spotPoly", divergent = TRUE,
                diverge_center = 0, image_id = "lowres", 
                swap_rownames = "symbol", color = "black",
                linewidth = 0.1)

We see that myofiber regions with higher Myh2 expression also have stronger spatial autocorrelation, while the injury site locally has some negative spatial autocorrelation.

The results are stored in the localResults field in the SFE object, with getters and setters analogous to reducedDims, but the name of the local method and the feature/gene for which the local method was run need to be specified as well.

lr <- localResult(sfe_tissue, type = "localmoran", feature = "Myh2", 
                  swap_rownames = "symbol")
head(lr)
#>                          Ii          E.Ii     Var.Ii       Z.Ii Pr(z != E(Ii))
#> AAACATTTCCCGGATT 2.12545883 -0.0012237995 0.37891181 3.45488511   5.505274e-04
#> AAACCTAAGCAGCCGG 3.33088903 -0.0038553468 0.59334776 4.32920244   1.496503e-05
#> AAACGAGACGGTTGAT 0.30430735 -0.0009817045 0.15152269 0.78428227   4.328745e-01
#> AAACGGGCGTACGGGT 4.69775712 -0.0063342069 0.97242481 4.77032239   1.839313e-06
#> AAACTCGGTTCGCAAT 0.01991573 -0.0002817611 0.04351933 0.09681804   9.228709e-01
#> AAACTGCTGGCTCCAA 0.55285063 -0.0009817045 0.15152269 1.42278566   1.547983e-01
#>                       mean    median     pysal    -log10p -log10p_adj
#> AAACATTTCCCGGATT High-High High-High High-High 3.25922109    2.657161
#> AAACCTAAGCAGCCGG High-High High-High High-High 4.82492235    3.979824
#> AAACGAGACGGTTGAT   Low-Low   Low-Low   Low-Low 0.36363800    0.000000
#> AAACGGGCGTACGGGT High-High High-High High-High 5.73534434    4.890246
#> AAACTCGGTTCGCAAT High-High High-High High-High 0.03485905    0.000000
#> AAACTGCTGGCTCCAA   Low-Low   Low-Low   Low-Low 0.81023381    0.000000

It is interesting to see how spatial autocorrelation relates to gene expression level, much as finding how variance relates to mean in the expression of each gene, which usually indicates overdispersion compared to Poisson in scRNA-seq and Visium data:

df <- data.frame(myh2 = logcounts(sfe_tissue)[rowData(sfe_tissue)$symbol == "Myh2",],
                 Ii = localResult(sfe_tissue, "localmoran", "Myh2", 
                                  swap_rownames = "symbol")[,"Ii"])
ggplot(df, aes(myh2, Ii)) + geom_point(alpha = 0.3) +
    geom_vline(xintercept = mean(df$myh2), linetype = 2, color = "gray") +
    geom_hline(yintercept = 0, linetype = 2, color = "gray") +
    labs(x = "Myh2 (log counts)", y = "localmoran")

For this gene, Visium spots with higher expression also tend to have higher local Moran’s I, but this may or may not apply to other genes. The vertical dashed line marks the mean gene expression; note the subtraction of mean in the expression for global and local Moran’s I, leading spots close to mean to have local Moran’s I close to 0.

Local spatial analyses often return a matrix or data frame. The plotLocalResult() function has a default column for each local spatial method, but other columns can be plotted as well. Use the localResultAttrs() function to see which columns are present, and use the attribute argument to specify which column to plot.

localResultAttrs(sfe_tissue, "localmoran", "Myh2", swap_rownames = "symbol")
#>  [1] "Ii"             "E.Ii"           "Var.Ii"         "Z.Ii"          
#>  [5] "Pr(z != E(Ii))" "mean"           "median"         "pysal"         
#>  [9] "-log10p"        "-log10p_adj"

Some local spatial methods return p-values at each location, in a column with name like Pr(z != E(Ii)), where the test is two sided (default, can be changed with the alternative argument in runUnivariate() which is passed to the relevant underlying function in spdep). Negative log of the p-value is computed to facilitate visualization (smaller or more significant p-values are plotted as higher values), and the p-value is corrected for multiple hypothesis testing with p.adjustSP() in spdep, where the number of tests is the number of neighbors of each location rather than the total number of locations (-log10p_adj).

plotLocalResult(sfe_tissue, "localmoran", features = "Myh2", 
                colGeometryName = "spotPoly", attribute = "-log10p_adj", divergent = TRUE,
                diverge_center = -log10(0.05), swap_rownames = "symbol",
                image_id = "lowres", color = "black",
                linewidth = 0.1)

In this plot and all following plots of p-values, a divergent palette is used to show locations that are significant after adjusting for multiple testing and those that are not significant in different colors. The center of the divergent palette is p = 0.05, so the brown spots are significant while a dark blue means really not significant.

The “pysal” column shows the type of neighborhood, such as whether low value is near other low values, or high value is near other high values.

plotLocalResult(sfe_tissue, "localmoran", features = "Myh2", 
                colGeometryName = "spotPoly", attribute = "pysal", 
                swap_rownames = "symbol", image_id = "lowres", color = "black",
                linewidth = 0.1)

Exercises

1. Compute local spatial heteroscedasticity (LOSH) (J. Keith Ord and Getis 2012) on Myh2 and plot the results. A sequential palette is appropriate.
2. Which other columns are returned by LOSH? Plot one of them in space. See documentation of spdep::LOSH() for the meanings of the other columns.
3. How does the spatial pattern of LOSH compare to that of local Moran for the same genes?

Further reading

1. Getis-Ord Gi* (J. K. Ord and Getis 1995) is another commonly used local spatial statistic which identifies hotspots (high values clustered together in space) and coldspots (low values clostered in space). These vignettes use Getis-Ord Gi*.
2. The Moran scatter plot (Anselin 1996) is another ESDA tool. Voyager has a special function to plot the results moranPlot(). See these vignettes for examples of Moran scatter plot applied to spatial transcriptomics.
3. Analogous to the Moran scatter plot is the Geary scatter plot (not yet implemented in Voyager) proposed in (Griffith and Chun 2022) which is said to better detect local negative spatial autocorrelation. This paper also includes other considerations on Moran’s I and Geary’s C.

Lee’s L

Lee’s L (Lee 2001) was developed from relating Moran’s I to Pearson correlation, and is defined as

where $n$ is the number of spots or locations, $i$ and $j$ are different locations, or spots in the Visium context, $x$ and $y$ are variables with values at each location, and $w_{ij}$ is a spatial weight, which can be inversely proportional to distance between spots or an indicator of whether two spots are neighbors, subject to various definitions of neighborhood. The Giotto package has implemented something like Lee’s L (Dries et al. 2021).

Here we compute Lee’s L for top highly variagle genes (HVGs) in this dataset:

hvgs <- getTopHVGs(sfe_tissue, fdr.threshold = 0.01)

Because bivariate global results can have very different formats (matrix for Lee’s L and lists for many other methods), the results are not stored in the SFE object. The calculateBivariate() function is used to perform all bivariate analyses. Analogous to runUnivariate() there is runBivariate() which stores the results in the SFE object, but it only applies to local bivariate methods whose results have a more uniform format and are stored in the localResults field just like local univariate results.

res <- calculateBivariate(sfe_tissue, type = "lee", feature1 = hvgs)

This gives a spatially informed correlation matrix among the genes, which can be plotted as a heatmap:

pal_rng <- getDivergeRange(res)
pal <- scico(256, begin = pal_rng[1], end = pal_rng[2], palette = "vik")

pheatmap(res, color = pal, show_rownames = FALSE, 
         show_colnames = FALSE, cellwidth = 1, cellheight = 1,
         treeheight_col = 0, treeheight_row = 0)

Some coexpression blocks can be seen. Note that unlike in Pearson correlation, the diagonal is not 1, because

which is approximated the ratio between the variance of spatially lagged $x$ and variance of $x$. Because the spatial lag introduces smoothing, the spatial lag reduced variance, making the diagonal less than 1. This is the spatial smoothing scalar (SSS), and Moran’s I is approximately Pearson correlation between $X$ and spatially lagged $X$ ($\tilde X$) multiplied by SSS:

Similarly for Lee’s L, as shown in (Lee 2001),

With more spatial clustering, the variance is less reduced by the spatial lag, leading to a larger SSS. Hence when both $X$ and $Y$ are spatially distributed like salt and pepper while strongly correlated, Lee’s L will be low because the lack of spatial autocorrelation leads to a small SSS.

Weighted correlation network analysis (WGCNA) (Langfelder and Horvath 2008) is a time honored method to find gene co-expression modules, and it can take any correlation matrix. Then it would be interesting to apply WGCNA to the Lee’s L matrix to identify spatially informed gene co-expression modules.

Exercises

Local Lee’s L is analogous to local Moran’s I – a disaggregated form of Lee’s L showing the contribution of each spot to the global Lee’s L. This is derived in (Lee 2001) after global Lee’s L is shown.

Run local Lee’s L on two genes of your choice. You can use the myofiber type marker genes Myh7, Myh2, Myh4, and Myh1. Then plot the results in space. Hint: Use localResultFeatures() to find what name the results are stored under. How would you interpret the results?

See this vignette on bivariate methods in Voyager applied to this dataset, including local Lee’s L and a bivariate version of local Moran’s I.

Non-spatial PCA

First we run regular principal component analysis (PCA), to compare to a type of spatially informed PCA known as MULTISPATI PCA (Stéphane Dray, Saı̈d, and Débias 2008).

hvgs2 <- getTopHVGs(dec, n = 2000)

sfe_tissue <- runPCA(sfe_tissue, ncomponents = 20, subset_row = hvgs2,
                     exprs_values = "logcounts", scale = TRUE,
                     BSPARAM = IrlbaParam())

Use the elbow plot to see variance explained by each PC:

ElbowPlot(sfe_tissue)

Plot top gene loadings in each PC, which are the contribution of each gene to each PC:

plotDimLoadings(sfe_tissue, swap_rownames = "symbol")

Plot the first 4 PCs in space

spatialReducedDim(sfe_tissue, "PCA", 4, divergent = TRUE, diverge_center = 0,
                  image_id = "lowres")

The first PC separates the leukocyte infiltrated injury site from the myofibers, while PC2 and PC3 tease out the muscle tendon junctions.

MULTISPATI PCA

Spatially informed dimension reduction is actually not new, and dates back to at least 1985, with Wartenberg’s crossover of Moran’s I and PCA (Wartenberg 1985), which was generalized and further developed as MULTISPATI PCA (Stéphane Dray, Saı̈d, and Débias 2008), implemented in the adespatial package on CRAN. In short, while PCA tries to maximize the variance explained by each PC, MULTISPATI maximizes the product of Moran’s I and variance explained. Also, while all the eigenvalues from PCA are non-negative, because the covariance matrix is positive semidefinite, MULTISPATI can give negative eigenvalues, which represent negative spatial autocorrelation, which can be present and interesting but is not as common as positive spatial autocorrelation and is often masked by the latter (Griffith 2019).

In single cell -omics conventions, let $X$ denote a gene count matrix whose columns are cells or Visium spots and whose rows are genes, with $n$ columns. Let $W$ denote the row normalized $n\times n$ adjacency matrix of the spatial neighborhood graph of the cells or Visium spots, which does not have to be symmetric. MULTISPATI diagonalizes a symmetric matrix

However, the implementation in adespatial is more general and can be used for other multivariate analyses in the duality diagram paradigm, such as correspondence analysis; the equation above is simplified just for PCA, without having to introduce the duality diagram here.

Here we compute MULTISPATI PCA, with the 20 most positive and 20 most negative eigenvalues.

sfe_tissue <- runMultivariate(sfe_tissue, "multispati", colGraphName = "visium",
                              nfposi = 20, nfnega = 20, subset_row = hvgs2)

Then plot the most positive and most negative eigenvalues. Note that the eigenvalues here are not variance explained. Instead, they are the product of variance explained and Moran’s I. So the most positive eigenvalues correspond to eigenvectors that simultaneously explain more variance and have large positive Moran’s I. The most negative eigenvalues correspond to eigenvectors that simultaneously explain more variance and have negative Moran’s I.

ElbowPlot(sfe_tissue, nfnega = 20, reduction = "multispati")

Here the positive eigenvalues drop sharply after, PC1, while none of the negative eigenvalues seem noteworthy. However, in spatial transcriptomics datasets with single cell resolution, there can be very negative eigenvalues and the corresponding PC is biologically relevant, see this vignette on mouse liver data.

What do these components mean? Each component is a linear combination of genes to maximize the product of variance explained and Moran’s I. The second component maximizes this product provided that it’s orthogonal to the first component, and so on. As the loss in variance explained is usually not huge, these components can be considered axes along which spatially coherent groups of spots are separated from each other as much as possible according to expression of the highly variable genes, so in theory, clustering with positive MULTISPATI components should give more spatially coherent clusters. Because of the spatial coherence, MULTISPATI might be more robust to outliers.

plotDimLoadings(sfe_tissue, dims = 1:4, reduction = "multispati",
                swap_rownames = "symbol")

Plot the these PCs:

spatialReducedDim(sfe_tissue, "multispati", 4, divergent = TRUE, 
                  diverge_center = 0, image_id = "lowres")

Here unlike in non-spatial PCA, PC4 continues to be spatially structured.

set.seed(29)
sfe_tissue$clusts_nonspatial <- clusterCells(sfe_tissue, use.dimred = "PCA", 
                                             BLUSPARAM = NNGraphParam(
                                                 cluster.fun = "leiden",
                                                 cluster.args = list(
                                                     objective_function = "modularity",
                                                     resolution_parameter = 1
                                                 )
                                             ))

set.seed(29)
sfe_tissue$clusts_multispati <- clusterRows(reducedDim(sfe_tissue, "multispati")[,1:20], 
                                            BLUSPARAM = NNGraphParam(
                                                cluster.fun = "leiden",
                                                cluster.args = list(
                                                    objective_function = "modularity",
                                                    resolution_parameter = 1
                                                )
                                            ))

plotSpatialFeature(sfe_tissue, c("clusts_nonspatial", "clusts_multispati"), 
                   colGeometryName = "spotPoly") &
     guides(colour = guide_legend(override.aes = list(size=2), ncol = 2))

# non-spatial
sfe_tissue <- reducedDimMoransI(sfe_tissue, dimred = "PCA", components = 1:20)
# spatial
sfe_tissue <- reducedDimMoransI(sfe_tissue, dimred = "multispati", components = 1:40)

df_moran <- tibble(PCA = reducedDimFeatureData(sfe_tissue, "PCA")$moran_Vis5A[1:20],
                   MULTISPATI_pos = 
                       reducedDimFeatureData(sfe_tissue, "multispati")$moran_Vis5A[1:20],
                   MULTISPATI_neg = 
                       reducedDimFeatureData(sfe_tissue,"multispati")$moran_Vis5A[21:40] |> 
                       rev(),
                   index = 1:20)

data("ditto_colors")

# Takes a while if not using optimized BLAS
(mb <- moranBounds(colGraph(sfe_tissue, "visium")))
#>       Imin       Imax 
#> -0.5762132  1.0021884

df_moran |> 
    pivot_longer(cols = -index, values_to = "value", names_to = "name") |> 
    ggplot(aes(index, value, color = name)) +
    geom_line() +
    scale_color_manual(values = ditto_colors) +
    geom_hline(yintercept = 0, color = "gray") +
    geom_hline(yintercept = mb, linetype = 2, color = "gray") +
    scale_y_continuous(breaks = scales::breaks_pretty()) +
    scale_x_continuous(breaks = scales::breaks_width(5)) +
    labs(y = "Moran's I", color = "Type", x = "Component")

Exercises

Thought experiment: suppose you perform standard PCA and MULTISPATI PCA on your city, such as on buildings, demographics, or when the city is divided into pixels, on whichever spatial features you find relevant to your life. What would the principal components look like? How would MULTISPATI PC’s differ from the standard PC’s?

Further reading

1. The other multivariate spatial method in Voyager as of Bioc 3.17 is a multivariate generalization of local Geary’s C (Anselin 2019). See this vignette on its application to this dataset.
2. Unlike similar EDA packages for spatial -omics data, Voyager is extensible, so you can make the uniform user interface run other spatial methods, just like in Tidymodels. See this vignette on extending Voyager.
3. There are spatially informed dimension reduction methods designed for spatial -omics data, although they tend to only consider positive spatial autocorrelation. For example, see (Shang and Zhou 2022) and (Velten et al. 2022).
4. This paper discusses many other types of multivariate spatial analyses in ecology, besides MULTISPATI PCA (S. Dray et al. 2012).
5. MULTISPATI PCA can be thought of as in between two extremes. One extreme is standard PCA, which diagonalizes the covariance matrix. The other extreme is the Moran eigen map (MEM) (Griffith 1996), which only uses the spatial weights matrix, without the data matrix. MEM’s are made by diagonalizing a double centered spatial weights matrix. The first eigenvector has values that make the largest possible Moran’s I given the spatial neighborhood graph. The second eigenvector also maximizes Moran’s I given that it’s orthogonal to the first eigenvector, and so on. The eigenvalues are Moran’s I multiplied by a constant. The last eigenvector has the most negative Moran’s I given the spatial neighborhood graph. These eigenvectors, or MEM’s, represent spatial structures of different length scales, which can be selected and used as covariates in regression to account for spatial autocorrelation, in a procedure called spatial filtering (Griffith 2000; Griffith and Peres-Neto 2006). See this vignette of adespatial for examples in ecology.