Introduction
Local Geary’s C (Anselin 1995) is defined as:
\[ c_i = \sum_jw_{ij}(x_i - x_j)^2, \]
where \(w_{ij}\)s are spatial weights from location \(i\) to location \(j\) and \(x\) is a variable at spatial location. This is generalized to multiple variables in (Anselin 2019):
\[ c_{k,i} = \sum_{v=1}^k c_{v,i}, \]
where there are \(k\) variables. This is essentially a spatially weighted sum of squared distances between locations in feature space. This vignette demonstrates usage of multivariate local Geary’s C.
Here we load the packages used:
library(Voyager)
library(SFEData)
library(SpatialFeatureExperiment)
library(scater)
library(scran)
library(ggplot2)
library(spdep)
theme_set(theme_bw())QC was performed in another vignette, so this vignette will not plot QC metrics.
(sfe <- McKellarMuscleData("full"))
#> see ?SFEData and browseVignettes('SFEData') for documentation
#> loading from cache
#> class: SpatialFeatureExperiment
#> dim: 15123 4992
#> metadata(0):
#> assays(1): counts
#> rownames(15123): ENSMUSG00000025902 ENSMUSG00000096126 ...
#> ENSMUSG00000064368 ENSMUSG00000064370
#> rowData names(6): Ensembl symbol ... vars cv2
#> colnames(4992): AAACAACGAATAGTTC AAACAAGTATCTCCCA ... TTGTTTGTATTACACG
#> TTGTTTGTGTAAATTC
#> colData names(12): barcode col ... prop_mito in_tissue
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> spatialCoords names(2) : imageX imageY
#> imgData names(1): sample_id
#>
#> unit: full_res_image_pixels
#> Geometries:
#> colGeometries: spotPoly (POLYGON)
#> annotGeometries: tissueBoundary (POLYGON), myofiber_full (POLYGON), myofiber_simplified (POLYGON), nuclei (POLYGON), nuclei_centroid (POINT)
#>
#> Graphs:
#> Vis5A:The image can be added to the SFE object and plotted behind the geometries, and needs to be flipped to align to the spots because the origin is at the top left for the image but bottom left for geometries.
if (!file.exists("tissue_lowres_5a.jpeg")) {
download.file("https://raw.githubusercontent.com/pachterlab/voyager/main/vignettes/tissue_lowres_5a.jpeg",
destfile = "tissue_lowres_5a.jpeg")
}
sfe <- addImg(sfe, file = "tissue_lowres_5a.jpeg", sample_id = "Vis5A",
image_id = "lowres",
scale_fct = 1024/22208)
#> Warning: The `file` argument of `addImg()` is deprecated as of SpatialFeatureExperiment
#> 1.6.0.
#> ℹ Please use the `imageSource` argument instead.
#> ℹ The deprecated feature was likely used in the SpatialFeatureExperiment
#> package.
#> Please report the issue at
#> <https://github.com/pachterlab/SpatialFeatureExperiment/issues>.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
sfe <- mirrorImg(sfe, sample_id = "Vis5A", image_id = "lowres")
sfe_tissue <- sfe[,colData(sfe)$in_tissue]
sfe_tissue <- sfe_tissue[rowSums(counts(sfe_tissue)) > 0,]
sfe_tissue <- logNormCounts(sfe_tissue)
colGraph(sfe_tissue, "visium") <- findVisiumGraph(sfe_tissue)Gene expression
Here we compute multivariate local C for top highly variagle genes (HVGs) in this dataset:
hvgs <- getTopHVGs(sfe_tissue, fdr.threshold = 0.01)
sfe_tissue <- runMultivariate(sfe_tissue, "localC_perm_multi", subset_row = hvgs)The results are stored in reducedDim although it’s not
really a dimension reduction. It can also go into colData
if dest = "colData". The test is two sided, but the
alternative argument can be set to “greater” to only test
for positive spatial autocorrelation and “less” for negative spatial
autocorrelation.
names(reducedDim(sfe_tissue, "localC_perm_multi"))
#> [1] "localC_perm_multi" "E.Ci" "Var.Ci"
#> [4] "Z.Ci" "Pr(z != E(Ci))" "Pr(z != E(Ci)) Sim"
#> [7] "Pr(folded) Sim" "Skewness" "Kurtosis"
#> [10] "-log10p Sim" "-log10p_adj Sim" "cluster"
spatialReducedDim(sfe_tissue, "localC_perm_multi", c(1, 12),
image_id = "lowres", maxcell = 5e4)
In Geary’s C, a value below 1 indicates positive spatial autocorrelation and a value above 1 indicates negative spatial autocorrelation. Local Geary’s C is not scaled, but from the square difference expression, a low value means a more homogeneous neighborhood and a high value means a more heterogeneous neighborhood. Here considering all 341 top HVGs, the muscle tendon junction and the unjury site are more heterogeneous, which is detected as negative cluster.
Permutation testing was performed, although Anselin noted that the pseudo-p-values should only be taken as indicative of interesting regions and should not be interpreted in a strict sense.
spatialReducedDim(sfe_tissue, "localC_perm_multi", c(11, 12),
image_id = "lowres", maxcell = 5e4,
divergent = TRUE, diverge_center = -log10(0.05))
Warm colors indicate adjusted p < 0.05. This should be interpreted along with the clusters. In this dataset, there are interestingly homogeneous regions in the myofibers, and an interestingly heterogeneous region in the injury site. Most of the significant regions are positive cluster, but the center of the injury site is significant and is negative cluster.
Top principal components
Because multivariate local Geary’s C is a spatially weighted sum of squared distances between locations in feature space, it’s affected by the curse of dimensionality when used on a large number of features, when uniformly distributed data points in higher dimensions become more equidistant to each other with increasing number of dimensions. However, real data is not uniformly distributed and can have a much smaller effective dimension than the number of features, as many genes are co-regulated. Anselin suggested using the main principal components, but the issue of curse of dimensionality remains to be further investigated. Furthermore, as the cosine and Manhattan distances have been suggested to mitigate curse of dimensionality, I wonder what if I use these instead of the Euclidean distance in feature space for multivariate local Geary’s C.
So here we perform multivariate local Geary’s C on the top PCs:
sfe_tissue <- runPCA(sfe_tissue, ncomponents = 20, scale = TRUE)
ElbowPlot(sfe_tissue)
What percentage of variance is explained by the top 20 PCs?
sum(attr(reducedDim(sfe_tissue, "PCA"), "percentVar"))
#> [1] 38.8627
out <- localC_perm(reducedDim(sfe_tissue, "PCA"),
listw = colGraph(sfe_tissue, "visium"))
out <- Voyager:::.localCpermmulti2df(out,
nb = colGraph(sfe_tissue, "visium")$neighbours,
p.adjust.method = "BH")
reducedDim(sfe_tissue, "localC_PCs", withDimnames = FALSE) <- out
spatialReducedDim(sfe_tissue, "localC_PCs", c(1, 12),
image_id = "lowres", maxcell = 5e4)
spatialReducedDim(sfe_tissue, "localC_PCs", c(11, 12),
image_id = "lowres", maxcell = 5e4,
divergent = TRUE, diverge_center = -log10(0.05))
The area that seem significant from the permutation test is larger than that from the HVGs, and the area considered negative clusters is smaller. The significant regions are pretty much all positive cluster. Do the differences in results have anything to do with curse of dimensionality? Twenty dimensions can still exhibit curse of dimensionality, but over 300 HVGs here would be worse. Or is it that we lose a lot of information, including negative spatial autocorrelation, by only using 20 PCs?
Session info
sessionInfo()
#> R version 4.4.0 (2024-04-24)
#> Platform: x86_64-apple-darwin20
#> Running under: macOS Ventura 13.6.6
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: UTC
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] spdep_1.3-3 sf_1.0-16
#> [3] spData_2.3.0 scran_1.32.0
#> [5] scater_1.32.0 ggplot2_3.5.1
#> [7] scuttle_1.14.0 SingleCellExperiment_1.26.0
#> [9] SummarizedExperiment_1.34.0 Biobase_2.64.0
#> [11] GenomicRanges_1.56.0 GenomeInfoDb_1.40.0
#> [13] IRanges_2.38.0 S4Vectors_0.42.0
#> [15] BiocGenerics_0.50.0 MatrixGenerics_1.16.0
#> [17] matrixStats_1.3.0 SFEData_1.6.0
#> [19] Voyager_1.6.0 SpatialFeatureExperiment_1.6.1
#>
#> loaded via a namespace (and not attached):
#> [1] bitops_1.0-7 filelock_1.0.3
#> [3] tibble_3.2.1 R.oo_1.26.0
#> [5] lifecycle_1.0.4 edgeR_4.2.0
#> [7] lattice_0.22-6 magrittr_2.0.3
#> [9] limma_3.60.0 sass_0.4.9
#> [11] rmarkdown_2.27 jquerylib_0.1.4
#> [13] yaml_2.3.8 metapod_1.12.0
#> [15] sp_2.1-4 DBI_1.2.2
#> [17] RColorBrewer_1.1-3 abind_1.4-5
#> [19] zlibbioc_1.50.0 purrr_1.0.2
#> [21] R.utils_2.12.3 RCurl_1.98-1.14
#> [23] rappdirs_0.3.3 GenomeInfoDbData_1.2.12
#> [25] ggrepel_0.9.5 irlba_2.3.5.1
#> [27] terra_1.7-71 units_0.8-5
#> [29] RSpectra_0.16-1 dqrng_0.4.0
#> [31] pkgdown_2.0.9 DelayedMatrixStats_1.26.0
#> [33] codetools_0.2-20 DropletUtils_1.24.0
#> [35] DelayedArray_0.30.1 tidyselect_1.2.1
#> [37] farver_2.1.2 UCSC.utils_1.0.0
#> [39] memuse_4.2-3 ScaledMatrix_1.12.0
#> [41] viridis_0.6.5 BiocFileCache_2.12.0
#> [43] jsonlite_1.8.8 BiocNeighbors_1.22.0
#> [45] e1071_1.7-14 systemfonts_1.1.0
#> [47] dbscan_1.1-12 tools_4.4.0
#> [49] ggnewscale_0.4.10 ragg_1.3.2
#> [51] Rcpp_1.0.12 glue_1.7.0
#> [53] gridExtra_2.3 SparseArray_1.4.3
#> [55] xfun_0.44 EBImage_4.46.0
#> [57] dplyr_1.1.4 HDF5Array_1.32.0
#> [59] withr_3.0.0 BiocManager_1.30.23
#> [61] fastmap_1.2.0 boot_1.3-30
#> [63] rhdf5filters_1.16.0 bluster_1.14.0
#> [65] fansi_1.0.6 digest_0.6.35
#> [67] rsvd_1.0.5 R6_2.5.1
#> [69] mime_0.12 textshaping_0.3.7
#> [71] colorspace_2.1-0 wk_0.9.1
#> [73] jpeg_0.1-10 RSQLite_2.3.6
#> [75] R.methodsS3_1.8.2 utf8_1.2.4
#> [77] generics_0.1.3 data.table_1.15.4
#> [79] class_7.3-22 httr_1.4.7
#> [81] htmlwidgets_1.6.4 S4Arrays_1.4.0
#> [83] pkgconfig_2.0.3 scico_1.5.0
#> [85] gtable_0.3.5 blob_1.2.4
#> [87] XVector_0.44.0 htmltools_0.5.8.1
#> [89] fftwtools_0.9-11 scales_1.3.0
#> [91] png_0.1-8 SpatialExperiment_1.14.0
#> [93] knitr_1.46 rjson_0.2.21
#> [95] curl_5.2.1 proxy_0.4-27
#> [97] cachem_1.1.0 rhdf5_2.48.0
#> [99] BiocVersion_3.19.1 KernSmooth_2.23-24
#> [101] parallel_4.4.0 vipor_0.4.7
#> [103] AnnotationDbi_1.66.0 desc_1.4.3
#> [105] s2_1.1.6 pillar_1.9.0
#> [107] grid_4.4.0 vctrs_0.6.5
#> [109] BiocSingular_1.20.0 dbplyr_2.5.0
#> [111] beachmat_2.20.0 sfheaders_0.4.4
#> [113] cluster_2.1.6 beeswarm_0.4.0
#> [115] evaluate_0.23 zeallot_0.1.0
#> [117] magick_2.8.3 cli_3.6.2
#> [119] locfit_1.5-9.9 compiler_4.4.0
#> [121] rlang_1.1.3 crayon_1.5.2
#> [123] labeling_0.4.3 classInt_0.4-10
#> [125] fs_1.6.4 ggbeeswarm_0.7.2
#> [127] viridisLite_0.4.2 deldir_2.0-4
#> [129] BiocParallel_1.38.0 munsell_0.5.1
#> [131] Biostrings_2.72.0 tiff_0.1-12
#> [133] Matrix_1.7-0 ExperimentHub_2.12.0
#> [135] patchwork_1.2.0 sparseMatrixStats_1.16.0
#> [137] bit64_4.0.5 Rhdf5lib_1.26.0
#> [139] KEGGREST_1.44.0 statmod_1.5.0
#> [141] highr_0.10 AnnotationHub_3.12.0
#> [143] igraph_2.0.3 memoise_2.0.1
#> [145] bslib_0.7.0 bit_4.0.5