Visium post-QC clustering and cell type annotation

Author

Harvard Chan Bioinformatics Core

Published

July 1, 2025

This code is in this revision.

Code
# 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)

0.1 Project details

  • Project: name_hbcXXXXX
  • PI: person name
  • Analyst: person in the core
  • Experiment: short description
  • Aim: short description

In this template, we assume you have already run QC on your visium data:

  • Remove genes that are contamination/poor quality-related (e.g. mitochondria or hemoglobin);
  • Keep only cells that have all QC metrics passing the chosen thresholds

Now, we will move to the analysis of Visium data, we will perform:

  • Normalization
  • Unsupervised sketch clustering
  • scRNA-seq data project

1 Normalize Data

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.

Code
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)

2 Unsupervised Clustering

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.

Code
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"
)
Code
# 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
Code
# 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)
Code
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")
)
Code
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

Code
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

3 Spatially-informed Clustering

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.

Code
# 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 seconds
Code
color_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

4 Cell Type Annotation

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.

Code
# 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.

Code
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))

5 Reference projection

Code
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)
Code
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 
Code
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)
Code
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

Code
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)

5.1 Methods

5.1.1 Citation

Code
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)'.
Code
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},
  }
Code
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")'.
Code
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")'.

5.1.2 Session Information

Code
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