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' assayobject_filt <-SketchData(object = object_filt,assay = assaytouse,ncells =10000,method ="LeverageScore",sketched.assay ="sketch")
Code
# switch analysis to sketched cellsDefaultAssay(object_filt) <-"sketch"# perform clustering workflowobject_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 usedobject_filt <-FindNeighbors(object_filt, assay ="sketch", reduction ="pca.sketch", dims =1:50)# you may want to tweak resolution parameter in your own dataobject_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 UMAPobject_filt <-RunUMAP(object_filt, reduction ="pca.sketch", reduction.name ="umap.sketch", return.model = T, dims =1:50)
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 dataobject_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 resolutionsobject_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
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 resultsROI <-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.
ref_subset <-inputRead(scrna_refF)# Check the label column you want to use from the scRNA-seq obs dataIdents(ref_subset) <-"subclass_label"counts <- ref_subset[["RNA"]]$countscluster <-as.factor(ref_subset$subclass_label)nUMI <- ref_subset$nCount_RNAlevels(cluster) <-gsub("/", "-", levels(cluster))cluster <-droplevels(cluster)# create the RCTD reference objectreference <-Reference(counts, cluster, nUMI)
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 objectROI <-AddMetaData(ROI, metadata = RCTD@results$results_df)
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 ROIexcitatory_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")'.
---title: "Visium post-QC clustering and cell type annotation"author: "Harvard Chan Bioinformatics Core"date: "`r Sys.Date()`"format: html: code-fold: true code-tools: true code-overflow: wrap df-print: paged highlight-style: pygments number-sections: true self-contained: true theme: default toc: true toc-location: left toc-expand: false lightbox: true page-layout: fullparams: project_file: ../information.R results_dir: ./results visium_postQCF: "../visium.qs" scrna_refF: "../allen_ref_subset.qs"---This code is in this  revision.```{r versioncheck, cache = FALSE, message = FALSE, warning=FALSE}# library(rstudioapi)# setwd(fs::path_dir(getSourceEditorContext()$path))stopifnot(R.version$major>= 4) # requires R4if (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)``````{r load_libraries, cache = FALSE, message = FALSE, warning=FALSE, echo=FALSE,}library(import)library(knitr)# analysis-specific packagelibrary(Seurat)library(SeuratWrappers)library(Banksy)library(quadprog)library(spacexr)# General data-wrangling library(glue)library(qs2)library(tidyverse)# Plottinglibrary(patchwork)library(ggprism)import::from(magrittr,set_colnames,set_rownames,"%<>%")invisible(list2env(params,environment()))source(project_file)ggplot2::theme_set(theme_prism(base_size = 12))# https://grafify-vignettes.netlify.app/colour_palettes.html# NOTE change colors here if you wishscale_colour_discrete <- function(...) scale_colour_manual(..., values = as.vector(grafify:::graf_palettes[["kelly"]]))scale_fill_discrete <- function(...) scale_fill_manual(..., values = as.vector(grafify:::graf_palettes[["kelly"]]))opts_chunk[["set"]]( cache = F, cache.lazy = FALSE, dev = c("png", "pdf"), error = TRUE, highlight = TRUE, message = FALSE, prompt = FALSE, tidy = FALSE, warning = FALSE, echo = T, fig.height = 4)# set seed for reproducibilityset.seed(1234567890L)options(future.globals.maxSize= 2000000000)inputRead <- function(f){ if(R.utils::isUrl(f)){f <- url(f)} if(sum(endsWith(f,c("rds","RDS")))>0){ return(readRDS(f)) }else if(sum(endsWith(f,c("qs","QS")))>0){ return(qs_read(f)) }else{ print("Check file extension and choose appropriate functions!") }}```## Project details - Project: `r project`- PI: `r PI`- Analyst: `r analyst`- Experiment: `r experiment`- Aim: `r aim`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 thresholdsNow, we will move to the analysis of Visium data, we will perform:- Normalization- Unsupervised sketch clustering- scRNA-seq data project# 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.```{r normalize}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)```# Unsupervised Clustering The authors of the `Seurat` package recommend the `Seurat` v5 **[sketch clustering](https://satijalab.org/seurat/articles/seurat5_sketch_analysis)** 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.```{r create sketch assay}object_filt <- FindVariableFeatures(object_filt)# we select 10,000 cells and create a new 'sketch' assayobject_filt <- SketchData( object = object_filt, assay = assaytouse, ncells = 10000, method = "LeverageScore", sketched.assay = "sketch")``````{r perform sketched clustering}# switch analysis to sketched cellsDefaultAssay(object_filt) <- "sketch"# perform clustering workflowobject_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 usedobject_filt <- FindNeighbors(object_filt, assay = "sketch", reduction = "pca.sketch", dims = 1:50)# you may want to tweak resolution parameter in your own dataobject_filt <- FindClusters(object_filt, cluster.name = "seurat_cluster.sketched", resolution = .65)# Use the same dimension of PCs for both FindNeighbors and run UMAPobject_filt <- RunUMAP(object_filt, reduction = "pca.sketch", reduction.name = "umap.sketch", return.model = T, dims = 1:50)``````{r project clusters}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"))``````{r visualize clusters,fig.height = 5,fig.width = 10}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 datasetDefaultAssay(object_filt) <- assaytouseIdents(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``````{r visualize clusters on image}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```# 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.```{r run banksy}# 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 dataobject_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 resolutionsobject_filt <- FindClusters(object_filt, cluster.name = "banksy_cluster", resolution = 0.5)``````{r}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```# 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.```{r}# change the list of clusters to your interest regions based on previous clustering resultsROI <-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.```{r sketch cortex}DefaultAssay(ROI) <- assaytouseROI <- 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"]]$countsROI_cells_hd <- colnames(ROI[["sketch"]])coords <- GetTissueCoordinates(ROI)[ROI_cells_hd, 1:2]# create the RCTD query objectquery <- SpatialRNA(coords, counts_hd, colSums(counts_hd))```# Reference projection ```{r load ref, prep for RCTD}ref_subset <- inputRead(scrna_refF)# Check the label column you want to use from the scRNA-seq obs dataIdents(ref_subset) <- "subclass_label"counts <- ref_subset[["RNA"]]$countscluster <- as.factor(ref_subset$subclass_label)nUMI <- ref_subset$nCount_RNAlevels(cluster) <- gsub("/", "-", levels(cluster))cluster <- droplevels(cluster)# create the RCTD reference objectreference <- Reference(counts, cluster, nUMI)``````{r run RCTD}RCTD <- create.RCTD(query, reference, max_cores = 6)RCTD <- run.RCTD(RCTD, doublet_mode = "doublet") # this command takes ~15 mins to run# add results back to Seurat objectROI <- AddMetaData(ROI, metadata = RCTD@results$results_df)``````{r project to all cortical cells}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```{r visualize labels}Idents(ROI) <- "full_first_type"cells <- CellsByIdentities(ROI)# Layered (starts with L), excitatory neurons in the ROIexcitatory_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)```## Methods### Citation```{r citations}citation("Seurat")citation("Banksy")citation("quadprog")citation("spacexr")```### Session Information```{r}sessionInfo()```