This code is in this 
revision.
# library(rstudioapi)
# setwd(fs::path_dir(getSourceEditorContext()$path))
stopifnot(R.version$major>= 4) # requires R4
if (compareVersion(R.version$minor,"3.1")<0) warning("We recommend >= R4.3.1")
stopifnot(compareVersion(as.character(BiocManager::version()), "3.18")>=0)
stopifnot(compareVersion(as.character(packageVersion("Seurat")), "5.1")>=0)In this template, we assume you have already run QC on your visium data:
Now, we will move to the analysis of Visium data, we will perform:
Normalization is important in order to make expression counts comparable across genes and/or sample. We note that the best normalization methods for spatial data are still being developed and evaluated. Here we use a standard log-normalization.
visiumHD_postQCobj <- inputRead(visium_postQCF)
assaytouse <- DefaultAssay(visiumHD_postQCobj)
message("Default assay: [",assaytouse,
"] used, please change it if another assay is of interest.")
object_filt <- NormalizeData(visiumHD_postQCobj, assay = assaytouse)The authors of the Seurat package recommend the Seurat v5 sketch clustering workflow because it exhibits improved performance, especially for identifying rare and spatially restricted groups.
Sketch-based analyses aim to “subsample” large datasets in a way that preserves rare populations. Here, we sketch the Visium HD dataset, perform clustering on the subsampled cells, and then project the cluster labels back to the full dataset.
object_filt <- FindVariableFeatures(object_filt)
# we select 10,000 cells and create a new 'sketch' assay
object_filt <- SketchData(
object = object_filt,
assay = assaytouse,
ncells = 10000,
method = "LeverageScore",
sketched.assay = "sketch"
)# switch analysis to sketched cells
DefaultAssay(object_filt) <- "sketch"
# perform clustering workflow
object_filt <- FindVariableFeatures(object_filt)
object_filt <- ScaleData(object_filt)
object_filt <- RunPCA(object_filt, assay = "sketch", reduction.name = "pca.sketch")
# default first 50 PCs are used
object_filt <- FindNeighbors(object_filt, assay = "sketch",
reduction = "pca.sketch",
dims = 1:50)
# you may want to tweak resolution parameter in your own data
object_filt <- FindClusters(object_filt, cluster.name = "seurat_cluster.sketched",
resolution = .65)Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 2696
Number of edges: 90359
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8909
Number of communities: 15
Elapsed time: 0 seconds# Use the same dimension of PCs for both FindNeighbors and run UMAP
object_filt <- RunUMAP(object_filt, reduction = "pca.sketch",
reduction.name = "umap.sketch", return.model = T,
dims = 1:50)object_filt <- ProjectData(
object = object_filt,
assay = assaytouse,
full.reduction = "full.pca.sketch",
sketched.assay = "sketch",
sketched.reduction = "pca.sketch",
umap.model = "umap.sketch",
dims = 1:50,
refdata = list(seurat_cluster.projected = "seurat_cluster.sketched")
)object_filt$seurat_cluster.projected <- object_filt$seurat_cluster.projected %>%
as.numeric %>% as.factor()
DefaultAssay(object_filt) <- "sketch"
Idents(object_filt) <- "seurat_cluster.sketched"
p1 <- DimPlot(object_filt, reduction = "umap.sketch", label = F, cols = 'polychrome') +
ggtitle("Sketched clustering") +
theme(legend.position = "bottom")+
guides(color = guide_legend(override.aes = list(size=4), ncol=5) )
# switch to full dataset
DefaultAssay(object_filt) <- assaytouse
Idents(object_filt) <- "seurat_cluster.projected"
p2 <- DimPlot(object_filt, reduction = "full.umap.sketch", label = F, raster = F,
cols = 'polychrome') +
ggtitle("Projected clustering") +
theme(legend.position = "bottom")+
guides(color = guide_legend(override.aes = list(size=4), ncol=5) )
p1 | p2
color_pal = Seurat::DiscretePalette(n = length(unique(object_filt$seurat_cluster.projected)),
palette = "polychrome")
names(color_pal) <- sort(unique(object_filt$seurat_cluster.projected))
image_seurat_clusters <- SpatialDimPlot(object_filt,
group.by = 'seurat_cluster.projected',
pt.size.factor = 8, cols = color_pal) +
theme(legend.position = "bottom",legend.title=element_blank())+
guides(fill = guide_legend(override.aes = list(size=4,name=""), ncol=5) )
image_seurat_clusters
BANKSY is another method for performing clustering.
Unlike Seurat, BANKSY takes into account not only an individual spot’s expression pattern but also the mean and the gradient of gene expression levels in a spot’s broader neighborhood. This makes it valuable for identifying and segmenting spatial tissue domains.
# lambda: (numeric between 0-1) Spatial weight parameter
# k_geom: (integer) kNN parameter - number of neighbors to use, default is 15
# Please consider tweaking those two parameters based on your understanding of your data
object_filt <- RunBanksy(object_filt, lambda = 0.8, verbose = T,
assay = assaytouse, slot = 'data', k_geom = 50)
object_filt <- RunPCA(object_filt, assay = "BANKSY",
reduction.name = "pca.banksy",
features = rownames(object_filt),
npcs = 30)
object_filt <- FindNeighbors(object_filt,
reduction = "pca.banksy",
dims = 1:30)
# again, do not forget to try different resolutions
object_filt <- FindClusters(object_filt, cluster.name = "banksy_cluster",
resolution = 0.5)Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 2696
Number of edges: 60852
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8852
Number of communities: 12
Elapsed time: 0 secondscolor_pal = Seurat::DiscretePalette(n = length(unique(object_filt$banksy_cluster)),
palette = "polychrome")
names(color_pal) <- sort(unique(object_filt$banksy_cluster))
image_banksy_clusters <- SpatialDimPlot(object_filt,
group.by = "banksy_cluster",
pt.size.factor = 7,
cols = color_pal)+
theme(legend.position = "bottom",legend.title=element_blank())+
guides(fill = guide_legend(override.aes = list(size=4,name=""), ncol=5) )
image_seurat_clusters | image_banksy_clusters
Perhaps we are particularly interested in understanding the organization of cell types in the cortical region of the brain.
We first subset our Seurat object to this region of interest.
# change the list of clusters to your interest regions based on previous clustering results
ROI <- subset(object_filt, seurat_cluster.projected %in% c(18, 19, 7, 2, 4))
color_pal = Seurat::DiscretePalette(n = length(unique(object_filt$seurat_cluster.projected)),
palette = "polychrome")
names(color_pal) <- sort(unique(object_filt$seurat_cluster.projected))
SpatialDimPlot(ROI, group.by = 'seurat_cluster.projected',
pt.size.factor = 8, cols = color_pal)+
theme(legend.position = "bottom",legend.title=element_blank())+
guides(fill = guide_legend(override.aes = list(size=4,name=""), ncol=5) )
To perform accurate annotation of cell types, we must also take into consideration that our 16 um spots may contain one or more cells each. The method Robust Cell Type Deconvolution (RCTD) has been shown to accurately annotate spatial data from a variety of technologies while taking into consideration that a single spot may exhibit multiple cell type profiles.
RCTD takes an scRNA-seq dataset as a reference and a spatial dataset as a query. For a reference, we use a subsampled version of the mouse scRNA-seq dataset from the Allen Brain Atlas. We use our cortex Seurat object as the spatial query. For computational efficiency, we sketch the spatial query dataset, apply RCTD to deconvolute the ‘sketched’ cortical cells and annotate them, and then project these annotations to the full cortical dataset.
DefaultAssay(ROI) <- assaytouse
ROI <- FindVariableFeatures(ROI)
ROI <- SketchData(
object = ROI,
ncells = 3000,
method = "LeverageScore",
sketched.assay = "sketch"
)
DefaultAssay(ROI) <- "sketch"
ROI <- ScaleData(ROI)
ROI <- RunPCA(ROI, assay = "sketch", reduction.name = "pca.ROI.sketch", verbose = T)
ROI <- FindNeighbors(ROI, reduction = "pca.ROI.sketch", dims = 1:50)
ROI <- RunUMAP(ROI, reduction = "pca.ROI.sketch", reduction.name = "umap.ROI.sketch", return.model = T, dims = 1:50, verbose = T)
counts_hd <- ROI[["sketch"]]$counts
ROI_cells_hd <- colnames(ROI[["sketch"]])
coords <- GetTissueCoordinates(ROI)[ROI_cells_hd, 1:2]
# create the RCTD query object
query <- SpatialRNA(coords, counts_hd, colSums(counts_hd))ref_subset <- inputRead(scrna_refF)
# Check the label column you want to use from the scRNA-seq obs data
Idents(ref_subset) <- "subclass_label"
counts <- ref_subset[["RNA"]]$counts
cluster <- as.factor(ref_subset$subclass_label)
nUMI <- ref_subset$nCount_RNA
levels(cluster) <- gsub("/", "-", levels(cluster))
cluster <- droplevels(cluster)
# create the RCTD reference object
reference <- Reference(counts, cluster, nUMI)RCTD <- create.RCTD(query, reference, max_cores = 6)
Astro CA1-ProS CA2-IG-FC CA3 Car3
500 500 101 442 500
CR CT SUB DG Endo L2 IT ENTl
208 500 500 500 500
L2 IT ENTm L2-3 IT CTX L2-3 IT ENTl L2-3 IT PPP L2-3 IT RHP
433 500 500 500 500
L3 IT ENT L4 RSP-ACA L4-5 IT CTX L5 IT CTX L5 PPP
500 500 500 500 254
L5 PT CTX L5-6 IT TPE-ENT L5-6 NP CTX L6 CT CTX L6 IT CTX
500 500 500 500 500
L6 IT ENTl L6b CTX L6b-CT ENT Lamp5 Micro-PVM
273 500 500 500 500
NP PPP NP SUB Oligo Pvalb SMC-Peri
500 398 500 500 285
Sncg Sst Sst Chodl SUB-ProS Vip
500 500 495 500 500
VLMC
152 RCTD <- run.RCTD(RCTD, doublet_mode = "doublet") # this command takes ~15 mins to run
# add results back to Seurat object
ROI <- AddMetaData(ROI, metadata = RCTD@results$results_df)ROI$first_type <- as.character(ROI$first_type)
ROI$first_type[is.na(ROI$first_type)] <- "Unknown"
ROI <- ProjectData(
object = ROI,
assay = assaytouse,
full.reduction = "pca.ROI",
sketched.assay = "sketch",
sketched.reduction = "pca.ROI.sketch",
umap.model = "umap.ROI.sketch",
dims = 1:50,
refdata = list(full_first_type = "first_type")
)We can see that the excitatory neurons are located in layers at varying cortical depths, as expected
Idents(ROI) <- "full_first_type"
cells <- CellsByIdentities(ROI)
# Layered (starts with L), excitatory neurons in the ROI
excitatory_names <- sort(grep("^L.* CTX", names(cells), value = TRUE))
SpatialDimPlot(ROI, cells.highlight = cells[excitatory_names],
cols.highlight = c("#FFFF00", "grey50"), facet.highlight = T,
combine = T, ncol = 4, pt.size.factor = 8)
citation("Seurat")To cite Seurat in publications, please use:
Hao et al. Dictionary learning for integrative, multimodal and
scalable single-cell analysis. Nature Biotechnology (2023) [Seurat
V5]
Hao and Hao et al. Integrated analysis of multimodal single-cell
data. Cell (2021) [Seurat V4]
Stuart and Butler et al. Comprehensive Integration of Single-Cell
Data. Cell (2019) [Seurat V3]
Butler et al. Integrating single-cell transcriptomic data across
different conditions, technologies, and species. Nat Biotechnol
(2018) [Seurat V2]
Satija and Farrell et al. Spatial reconstruction of single-cell gene
expression data. Nat Biotechnol (2015) [Seurat V1]
To see these entries in BibTeX format, use 'print(<citation>,
bibtex=TRUE)', 'toBibtex(.)', or set
'options(citation.bibtex.max=999)'.citation("Banksy")To cite BANKSY in publications use:
Singhal, V., Chou, N., et al. BANKSY: A Spatial Omics Algorithm that
Unifies Cell Typing and Tissue Domain Segmentation Preprint at
bioRxiv https://doi.org/10.1101/2022.04.14.488259 (2022)
A BibTeX entry for LaTeX users is
@Article{,
title = {BANKSY: A Spatial Omics Algorithm that Unifies Cell Typing and Tissue Domain Segmentation},
author = {Vipul Singhal and Nigel Chou and Joseph Lee and Jinyue Liu and Wan Kee Chock and Li Lin and Yun-Ching Chang and Erica Teo and Hwee Kuan Lee and Kok Hao Chen and Shyam Prabhakar},
journal = {bioRxiv},
year = {2022},
url = {https://www.biorxiv.org/content/10.1101/2022.04.14.488259},
}citation("quadprog")To cite package 'quadprog' in publications use:
dpodi/LINPACK) SobBATRpbAW<FcfCM (2019). _quadprog: Functions to
Solve Quadratic Programming Problems_. R package version 1.5-8,
<https://CRAN.R-project.org/package=quadprog>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {quadprog: Functions to Solve Quadratic Programming Problems},
author = {S original by Berwin A. Turlach R port by Andreas Weingessel <Andreas.Weingessel@ci.tuwien.ac.at> Fortran contributions from Cleve Moler dpodi/LINPACK)},
year = {2019},
note = {R package version 1.5-8},
url = {https://CRAN.R-project.org/package=quadprog},
}
ATTENTION: This citation information has been auto-generated from the
package DESCRIPTION file and may need manual editing, see
'help("citation")'.citation("spacexr")To cite package 'spacexr' in publications use:
Cable D (2025). _spacexr: SpatialeXpressionR: Cell type
identification and cell type-specific differential expression in
spatial transcriptomics_. R package version 2.2.1, commit
744153c633da7ff70d0988308a160633fcf44af7,
<https://github.com/dmcable/spacexr>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {spacexr: SpatialeXpressionR: Cell type identification and cell
type-specific differential expression in spatial
transcriptomics},
author = {Dylan Cable},
year = {2025},
note = {R package version 2.2.1, commit 744153c633da7ff70d0988308a160633fcf44af7},
url = {https://github.com/dmcable/spacexr},
}
ATTENTION: This citation information has been auto-generated from the
package DESCRIPTION file and may need manual editing, see
'help("citation")'.sessionInfo()R version 4.4.2 (2024-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.5
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/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: America/New_York
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] future_1.40.0 ggprism_1.0.5 patchwork_1.3.0
[4] lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1
[7] dplyr_1.1.4 purrr_1.0.4 readr_2.1.5
[10] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.2
[13] tidyverse_2.0.0 qs2_0.1.5 glue_1.8.0
[16] spacexr_2.2.1 quadprog_1.5-8 Banksy_1.5.6
[19] SeuratWrappers_0.4.0 Seurat_5.3.0 SeuratObject_5.1.0
[22] sp_2.2-0 knitr_1.50 import_1.3.2
loaded via a namespace (and not attached):
[1] RcppHungarian_0.3 RcppAnnoy_0.0.22
[3] splines_4.4.2 later_1.4.2
[5] R.oo_1.27.1 polyclip_1.10-7
[7] fastDummies_1.7.5 lifecycle_1.0.4
[9] aricode_1.0.3 doParallel_1.0.17
[11] globals_0.17.0 lattice_0.22-7
[13] MASS_7.3-65 magrittr_2.0.3
[15] plotly_4.10.4 rmarkdown_2.29
[17] yaml_2.3.10 remotes_2.5.0
[19] httpuv_1.6.16 sctransform_0.4.1
[21] spam_2.11-1 spatstat.sparse_3.1-0
[23] reticulate_1.42.0 cowplot_1.1.3
[25] pbapply_1.7-2 RColorBrewer_1.1-3
[27] abind_1.4-8 zlibbioc_1.52.0
[29] Rtsne_0.17 GenomicRanges_1.58.0
[31] R.utils_2.13.0 BiocGenerics_0.52.0
[33] GenomeInfoDbData_1.2.13 IRanges_2.40.1
[35] S4Vectors_0.44.0 ggrepel_0.9.6
[37] irlba_2.3.5.1 listenv_0.9.1
[39] spatstat.utils_3.1-3 goftest_1.2-3
[41] RSpectra_0.16-2 spatstat.random_3.3-3
[43] fitdistrplus_1.2-2 parallelly_1.43.0
[45] codetools_0.2-20 DelayedArray_0.32.0
[47] tidyselect_1.2.1 UCSC.utils_1.2.0
[49] farver_2.1.2 matrixStats_1.5.0
[51] stats4_4.4.2 spatstat.explore_3.4-2
[53] jsonlite_2.0.0 progressr_0.15.1
[55] ggridges_0.5.6 survival_3.8-3
[57] iterators_1.0.14 foreach_1.5.2
[59] dbscan_1.2.2 tools_4.4.2
[61] ica_1.0-3 Rcpp_1.0.14
[63] gridExtra_2.3 SparseArray_1.6.2
[65] xfun_0.52 MatrixGenerics_1.18.1
[67] GenomeInfoDb_1.42.3 withr_3.0.2
[69] BiocManager_1.30.25 fastmap_1.2.0
[71] digest_0.6.37 rsvd_1.0.5
[73] timechange_0.3.0 R6_2.6.1
[75] mime_0.13 colorspace_2.1-1
[77] scattermore_1.2 sccore_1.0.6
[79] tensor_1.5 dichromat_2.0-0.1
[81] spatstat.data_3.1-6 R.methodsS3_1.8.2
[83] generics_0.1.3 data.table_1.17.0
[85] httr_1.4.7 htmlwidgets_1.6.4
[87] S4Arrays_1.6.0 uwot_0.2.3
[89] pkgconfig_2.0.3 gtable_0.3.6
[91] lmtest_0.9-40 SingleCellExperiment_1.28.1
[93] XVector_0.46.0 htmltools_0.5.8.1
[95] dotCall64_1.2 scales_1.4.0
[97] Biobase_2.66.0 png_0.1-8
[99] SpatialExperiment_1.16.0 spatstat.univar_3.1-2
[101] rstudioapi_0.17.1 tzdb_0.5.0
[103] reshape2_1.4.4 rjson_0.2.23
[105] nlme_3.1-168 zoo_1.8-14
[107] KernSmooth_2.23-26 parallel_4.4.2
[109] miniUI_0.1.2 pillar_1.10.2
[111] grid_4.4.2 vctrs_0.6.5
[113] RANN_2.6.2 promises_1.3.2
[115] stringfish_0.16.0 xtable_1.8-4
[117] cluster_2.1.8.1 evaluate_1.0.3
[119] magick_2.8.6 cli_3.6.5
[121] compiler_4.4.2 rlang_1.1.6
[123] crayon_1.5.3 future.apply_1.11.3
[125] labeling_0.4.3 mclust_6.1.1
[127] plyr_1.8.9 stringi_1.8.7
[129] viridisLite_0.4.2 deldir_2.0-4
[131] lazyeval_0.2.2 spatstat.geom_3.3-6
[133] Matrix_1.7-3 RcppHNSW_0.6.0
[135] hms_1.1.3 shiny_1.10.0
[137] SummarizedExperiment_1.36.0 ROCR_1.0-11
[139] leidenAlg_1.1.5 igraph_2.1.4
[141] RcppParallel_5.1.10