package.list <- c("tidyverse", "semla", "Seurat", "patchwork", "RColorBrewer", "harmony")
invisible(lapply(package.list, function(element) {
library(element, character.only = TRUE)
}))
rm(package.list)Step-by-step discovery of Visium data analysis using embryonic heart sections.
2026 V4SDB Student Summer School
Visium workshop
In this workshop we will go step by step through the analysis of Visium data from 10x Genomics. All the code is provided so you are not asked to write it yourself, but to read each step carefully and understand what it does and why. We will load the data, filter it, process it, and cluster it on a single tissue section. If time allows, all the steps are then repeated on three sections of different developmental ages simultaneously, where more advanced analyses are suggested. If we do not get there together, you are welcome to work through it on your own and contact me with any questions.
What Visium measures
A Visium slide carries a grid of ~5,000 spots, each 55 µm across and spaced 100 µm apart (distance from center to center of spots). Every spot carries millions of capture oligos sharing one spatial barcode unique to that position. When a tissue section is placed on the slide and permeabilised, mRNA diffuses down into the spots and is captured by a poly-T sequence that binds the poly-A tail of the transcript. After sequencing, the barcode tells us where on the slide a transcript came from, and the rest of the read tells us which gene it was.
Two consequences shape everything that follows.
- First, a spot is a fixed position, not a cell: at 55 µm it typically covers a handful of cells, so what we measure is a local mixture.
- Second, poly-A capture is indiscriminate: any polyadenylated transcript is sequenced, whether or not we are interested in it. A large part of this workshop is deciding what to keep.
Where this data comes from
The sections we analyse today are a small subset of a published atlas of the developing human heart, produced in our lab:
Lázár, E.*, Mauron, R.*, Andrusivová, Ž.* et al. Spatiotemporal gene expression and cellular dynamics of the developing human heart. Nature Genetics 57, 2756–2771 (2025). doi:10.1038/s41588-025-02352-6 (*equal contribution)
The full study combined Visium spatial transcriptomics on 16 hearts (38 sections, 69,114 spots, PCW 6–12) with single-cell RNA-seq on 15 hearts (76,991 cells, PCW 5.5–14) and in situ sequencing of 150 transcripts. Unsupervised clustering of 25,208 spots from 17 sections yielded 23 spatial clusters corresponding to the major structural components of the developing heart, and deconvolution against the single-cell data mapped 72 fine-grained cell states onto tissue positions.
Today we work with three sections. Everything we do is a simplified version of the first steps of that pipeline.
Resources you can use during and after the workshop:
- Interactive viewer: https://hdcaheart.serve.scilifelab.se
- Full analysis code: https://github.com/rmauron/HDCA_heart_dev
- Processed data (Mendeley): https://doi.org/10.17632/fhtb99mdzd.1 and https://doi.org/10.17632/w65jtfsvpr.1
The tissue was donated after elective medical abortion with written informed consent, under Swedish ethical permit 2018/769-31.
By the end of the session you will have
On a single section:
- Loaded a Visium section into R and plotted a gene on the tissue
- Decided which genes to filter out, and justified why
- Normalized, scaled, and reduced the data
- Clustered the spots and mapped the clusters back onto the anatomy
- Found the genes that define each cluster, and compared them against the published atlas
On three sections of different ages:
- Repeated the whole pipeline on three sections at once
- Seen what a batch effect looks like, and what correcting it costs
- Looked for genes that change with developmental stage within a cluster
- Asked whether a variable gene is actually spatially structured
Practical notes
Grey boxes are code chunks. Run the chunk your cursor is in with Cmd/Ctrl + Shift + Enter, or a single line with Cmd/Ctrl + Enter.
Run them in order, top to bottom. Many steps deliberately overwrite an object with an updated version, so skipping one, or re-running an earlier one out of order, will produce confusing errors later. Three places to be aware of:
- The single-section analysis builds up an object called
object. - The three-section analysis uses a separate object called
se. The two are independent —objectis not affected by anything in the second half. - In section 7 we cluster twice and
FindClusters()writes to the sameseurat_clusterscolumn both times. We copy the result intoclusters_pcaandclusters_harmonyso that nothing is lost, but if you re-run only part of that chunk,seurat_clustersmay not hold what you expect.
Red text in the console is not always an error; warnings are common and usually harmless. If something breaks, say so and we will look at it together.
Set up
Set project library
When working in R, there is multiple ways to set up your working environment. For the sake of that workshop, you are provided with a fully working environment that contains the code, the data, and the R packages already installed. Below we are just specifying R where to look for the installed packages.
Load library
R packages might be installed, but for us to be able to use them later in the analysis, we have to attach them to our session. In R, we do it by using the library() function. This slightly more advanced way below, allow to call the function library() in a loop. You could also do: library(semla) library(Seurat) etc. This works equally well.
Set seed for project
Sometimes, some code and functions rely on some random drawing of numbers. For instance, if I tell R to sample 100 random cells from my single-cell RNA-seq data, everytime I run the function, I will have different randomly picked cells. It does what we ask, but for the sake of reproducibility, we want to have always the same cells picked when we run the same code the next day. To do so, we can specify a seed by using set.seed(). Some functions have that option in the function’s parameter (worth checking because it might not be used by default), but otherwise you can use set.seed() with any number you want in the code chunk where you need it.
set.seed(1)Loading my first Visium section
Our session is now fully ready to start. We have the packages loaded, we are ready to load some data.
To analyse a Visium section, we need to load the data. 10x Genomics provides a software (spaceranger) that takes the raw sequences and the microscopy images, and generates a count matrix (gene counts per barcode), and several files that help us align the microscopy image to the count matrix.
How do we load this count matrix in R?
In R, scientists have developed convenience packages to effortlessly load the spaceranger outputs, store them into an object class, and provide multiple functions to analyse, plot and else.
In the single-cell RNA seq world, such a convenience package is Seurat, and is very well established and maintained by a large team. They are at the origin of the object class that we call Seurat object.
When Spatial Transcriptomics appeared, there were no such toolkit, so other labs put effort into creating such a package. Today we will rely on the semla package developed in our lab (Lundeberg). semla, under the hood, uses the Seurat object class, but added a way to register the spatial information that single-cell data do not have (image, coordinates, etc.).
In the next steps, most functions we will use come from either semla or Seurat.
Load data & create infoTable from metadata + spaceranger outs
Let’s start!
Below we create a table that provides the path so R (and specific packages e.g., semla) knows where the spaceranger outputs are. We have here:
filtered_feature_bc_matrix.h5: the count datatissue_hires_image.png: an H&E image that can be aligned to the countstissue_positions_list.csv: a table that contains the info to align the pixel of the image and the Visium spotsscalefactors_json.json: a ratio value that help us visualize the image and the spots in the right scale.
data_root_directory <- file.path("./data/spaceranger/V10F24-105_C1")
samples <- Sys.glob(paths = file.path(data_root_directory, "filtered_feature_bc_matrix.h5"))
imgs <- Sys.glob(paths = file.path(data_root_directory, "spatial", "tissue_hires_image.png"))
spotfiles <- Sys.glob(paths = file.path(data_root_directory, "spatial", "tissue_positions_list.csv"))
json <- Sys.glob(paths = file.path(data_root_directory, "spatial", "scalefactors_json.json"))We create a table containing all these paths
infoTable <- tibble(samples,
imgs,
spotfiles,
json)We can inspect the table:
infoTable# A tibble: 1 × 4
samples imgs spotfiles json
<chr> <chr> <chr> <chr>
1 ./data/spaceranger/V10F24-105_C1/filtered_feature_bc_ma… ./da… ./data/s… ./da…We can now use ReadVisiumData() function from semla package to create an R-object with all the data we linked with their paths in the infoTable. semla uses the data structure from the well-known Seurat package, originally developed to facilitate the analysis of single-cell RNA-seq data. semla added, on top of this architecture, all the spatial information (image, coordinates, etc.) that are required for us to know where are the spots and to be able to later use this information in analyses and visualization.
object <- ReadVisiumData(infoTable)
object <- LoadImages(object) #we specify semla to load the image in the objectLet’s inspect what we have so far.
objectAn object of class Seurat
36601 features across 2234 samples within 1 assay
Active assay: Spatial (36601 features, 0 variable features)
1 layer present: countsWhat are we actually looking at?
Before we touch the data, it is worth spending two minutes on the tissue. Everything we do from here on produces coloured dots on this image, and those dots only mean something if you know what is underneath them.
This section is a PCW 10 human heart, cryosectioned at 10 µm. At this stage the four chambers are formed and the heart is growing and maturing rather than being built from scratch. Depending on the cutting plane you should be able to make out:
- Compact myocardium — the dense outer wall of the ventricles
- Trabeculated myocardium — the sponge-like inner ventricular wall, with deep recesses opening into the lumen
- Atrial wall — thinner, smoother, above the atrioventricular plane
- Interventricular septum — the wall between the two ventricles
- Valve region / endocardial cushions — pale, loosely packed mesenchyme at the atrioventricular plane
- Chamber lumens — empty space, often full of blood
- Epicardium — the thin layer at the outer surface
ImagePlot(object)
Two features of this image will explain a lot of what we see later. The lumens are holes: there is no tissue there, so spots covering them capture very little RNA. And the wall thickness varies enormously across the section, so the amount of tissue under a 55 µm spot is far from constant. Neither is a technical problem — both are anatomy — but both show up in every quality metric we are about to plot.
We can start orienting ourselves molecularly straight away. MYH6 and MYH7 are the atrial and ventricular myosin heavy chains, so plotting them together separates the upper and lower halves of the heart:
MapFeatures(
object = object,
slot = "counts",
features = c("MYH6", "MYH7"),
pt_size = 1,
image_use = "raw"
)
Answer
- MYH6 marks the atria, MYH7 the ventricles (strongest in the compact myocardium). They are largely mutually exclusive, which we can check numerically rather than by eye:
d <- FetchData(object, vars = c("MYH6", "MYH7"), layer = "counts")
cor(d$MYH6, d$MYH7, method = "spearman") # negative = complementary domains[1] -0.1332847Spots with neither are not myocardium: chamber lumens (little tissue at all), valve mesenchyme, and the epicardial surface.
TNNT2 covers all myocardium (atria and ventricles), COL1A1 picks out the valve region and mesenchyme, PECAM1 traces endocardium and vessels — a thin outline rather than a solid block.
MapFeatures(object, features = c("TNNT2", "COL1A1", "PECAM1"),
slot = "counts", pt_size = 1, ncol = 3)
Keep this image in mind
When we cluster the spots later, come back to this section. Clusters that trace a recognisable structure are believable. Clusters that scatter across unrelated regions usually reflect a technical gradient rather than biology.
Cleaning up a bit…
We see about twice as many unique genes as there are human protein-coding genes. Our data comes from poly-A capture, so every polyadenylated transcript gets sequenced — including a lot we did not set out to measure.
Four families are worth removing. All are very highly expressed, so they eat a large share of each spot’s counts, and all vary for reasons unrelated to the anatomy we want to map:
- **Hemoglobin (HB*)** — an embryonic heart is full of blood, and erythrocytes are packed with hemoglobin transcripts. Keep them and our clusters will report how much blood a spot contains rather than what tissue it is.
- Mitochondrial (MT-) — only 13 genes, but extremely abundant. In heart they track muscle density and metabolic state, and they also rise with RNA degradation across the slide. Technical gradient mixed with global cell state.
- Ribosomal proteins (RPL, RPS) — ~80 housekeeping genes, tightly correlated, often 10–25% of a spot’s counts. They move as a block and can dominate the first principal component with “ribosome content”.
- MALAT1 — a single nuclear lncRNA, often the most-detected feature of all. How much we capture depends mostly on permeabilisation, so it reads out the experiment more than the biology.
This is a project-specific choice, not a rule: if you studied erythropoiesis, hemoglobin would be your signal, not your noise.
In R, we can use simple regular expression or called regex(text patterns that match) to retrieve these genes:
pattern <- paste0(
"^MT-", # 13 mitochondrial genes: MT-ND*, MT-CO*, MT-ATP*, MT-CYB
"|^HB[ABDEGMQZ][0-9]?$", # HBA1/2, HBB, HBD, HBE1, HBG1/2, HBM, HBQ1, HBZ — not HBEGF, HBP1, HBS1L
"|^RP[LS][0-9]+[A-Z]?[0-9]?$|^RPL36AL$|^RPLP[0-9]$|^RPSA$", # ~85 ribosomal proteins — not the RPS6K kinases
"|^MALAT1$" # nuclear lncRNA, often the single most-detected feature
)
undesirable_genes <- grep(pattern, rownames(object), value = TRUE)
undesirable_genes <- setdiff(undesirable_genes, "RPL3L") # muscle-specific paralog, real cardiac marker
length(undesirable_genes) # 110[1] 110# check nothing was missed
grep("^HB|^RP|^MT-", setdiff(rownames(object), undesirable_genes), value = TRUE) [1] "RPS6KA1" "RPA2" "RPE65" "RPF1" "RPAP2"
[6] "RPRD2" "RPTN" "RPS6KC1" "RPIA" "RPRM"
[11] "RPE" "RPUSD3" "RPP14" "RPN1" "RPL34-AS1"
[16] "HBEGF" "RPP40" "RPP21" "RPS10-NUDT3" "RPF2"
[21] "HBS1L" "RPS6KA2" "RPS6KA2-IT1" "RPS6KA2-AS1" "RPA3"
[26] "RP9" "HBP1" "RP1L1" "RP1" "RPP25L"
[31] "RPP38-DT" "RPP38" "RPP30" "RPARP-AS1" "RPEL1"
[36] "RPS6KA4" "RPS6KB2" "RPS6KB2-AS1" "RPUSD4" "RPAP3"
[41] "RPH3A" "RPGRIP1" "RPS6KL1" "RPS6KA5" "RPUSD2"
[46] "RPAP1" "RPP25" "RPUSD1" "RPL3L" "RPGRIP1L"
[51] "RPH3AL" "RPA1" "RPAIN" "RPRML" "RPS6KB1"
[56] "RPTOR" "RPRD1A" "RPN2" "RPRD1B" "RPS19BP1"
[61] "RPS6KA3" "RPGR" "RP2" "RPS6KA6" "RPA4" # check the list of undesirable genes
undesirable_genes [1] "RPL22" "RPL11" "RPS8" "RPL5" "RPS27" "RPS7" "RPS27A"
[8] "RPL31" "RPL37A" "RPL32" "RPL15" "RPSA" "RPL14" "RPL29"
[15] "RPL24" "RPL22L1" "RPL39L" "RPL35A" "RPL9" "RPL34" "RPS3A"
[22] "RPL37" "RPS23" "RPS14" "RPL26L1" "RPS18" "RPS10" "RPL10A"
[29] "RPL7L1" "RPS12" "RPS20" "RPL7" "RPL30" "RPL8" "RPS6"
[36] "RPL35" "RPL12" "RPL7A" "RPS24" "RPLP2" "HBB" "HBD"
[43] "HBG1" "HBG2" "HBE1" "RPL27A" "RPS13" "MALAT1" "RPS3"
[50] "RPS25" "RPS26" "RPL41" "RPL6" "RPLP0" "RPL21" "RPL10L"
[57] "RPS29" "RPL36AL" "RPS27L" "RPL4" "RPLP1" "RPS17" "HBZ"
[64] "HBM" "HBA2" "HBA1" "HBQ1" "RPS2" "RPS15A" "RPL13"
[71] "RPL26" "RPL23A" "RPL23" "RPL19" "RPL27" "RPL38" "RPL17"
[78] "RPS15" "RPL36" "RPS28" "RPL18A" "RPS16" "RPS19" "RPL18"
[85] "RPL13A" "RPS11" "RPS9" "RPL28" "RPS5" "RPS21" "RPL3"
[92] "RPS4X" "RPL36A" "RPL39" "RPL10" "RPS4Y1" "RPS4Y2" "MT-ND1"
[99] "MT-ND2" "MT-CO1" "MT-CO2" "MT-ATP8" "MT-ATP6" "MT-CO3" "MT-ND3"
[106] "MT-ND4L" "MT-ND4" "MT-ND5" "MT-ND6" "MT-CYB" We can get a list of all the other genes (the ones we want to keep), and see the difference compared to the original object.
genes_to_keep <- setdiff(rownames(object), undesirable_genes)
length(rownames(object)) - length(genes_to_keep)[1] 110Now we can proceed with a first filtering. We want to remove the undesirable genes, meaning we want to keep only all the other genes. We use the function SubsetSTData() from the semla package, which ensure that the spatial information is also modified correctly.
object <- SubsetSTData(object, features = genes_to_keep)
objectAn object of class Seurat
36491 features across 2234 samples within 1 assay
Active assay: Spatial (36491 features, 0 variable features)
1 layer present: countsNote that when using SubsetSTData() we assigned the result of the function to the same original name object. That action overwrites the previous state of our object. If you want to keep the unfiltered and filtered versions side by side, assign the result to a different variable name.
More systematic filtering
The extra transcripts we sequence and gather in our object do not limit themselves to the 4 categories we successfully removed above. Databases allow us to annotate the remaining genes we sequenced into some categories. Here, we will use the GRCh38-2020-A genome reference which lets us find out which category each gene belongs to.
The complete GRCh38-2020-A genome reference is quite large so here we provide the single file of that reference that you can load in the next chunk.
You do not need to run the next chunk This chunk commented out is for you to see what has been done from the complete reference to get the file we will use next.
# package that help us load large files
library(rtracklayer)
# import the
gtf <- import("/srv/home/10x_references/refdata-gex-GRCh38-2020-A/genes/genes.gtf")
unique(gtf$type)
# [1] gene transcript exon CDS start_codon
# [6] stop_codon UTR Selenocysteine
annot <- gtf[gtf$type %in% c("gene")] |>
mcols() |>
as.data.frame() |>
dplyr::select(gene_id, gene_name, gene_type) |>
dplyr::distinct()
saveRDS(annot, "./data/GRCh38-2020-A_gene_biotypes.rds")Here we load the lighter version of the reference genome:
annot <- readRDS("./data/GRCh38-2020-A_gene_biotypes.rds")
head(annot) gene_id gene_name gene_type
1 ENSG00000243485 MIR1302-2HG lncRNA
2 ENSG00000237613 FAM138A lncRNA
3 ENSG00000186092 OR4F5 protein_coding
4 ENSG00000238009 AL627309.1 lncRNA
5 ENSG00000239945 AL627309.3 lncRNA
6 ENSG00000239906 AL627309.2 lncRNAunique(annot$gene_type) [1] "lncRNA" "protein_coding" "IG_V_pseudogene" "IG_V_gene"
[5] "IG_C_gene" "IG_J_gene" "TR_C_gene" "TR_J_gene"
[9] "TR_V_gene" "TR_V_pseudogene" "TR_D_gene" "IG_C_pseudogene"
[13] "TR_J_pseudogene" "IG_J_pseudogene" "IG_D_gene" Here we will use the information of the gene types to keep only the protein coding genes. This step is project dependent, in other cases you might want to also keep the long-non-coding RNA or any types of immunoglobulin transcripts.
# subset the reference genes to only the protein coding genes
annot_protein_coding <- subset(annot, gene_type == "protein_coding")
# fetch the gene names from our object that are annotated as protein coding genes in the reference genome
genes_to_keep_2 <- intersect(rownames(object), annot_protein_coding$gene_name)
# use the same subset idea as before
object <- SubsetSTData(object, features = genes_to_keep_2)
# correct metadata after subset
object$nCount_Spatial <- colSums(LayerData(object, "counts"))
object$nFeature_Spatial <- colSums(LayerData(object, "counts") > 0)
objectAn object of class Seurat
19280 features across 2234 samples within 1 assay
Active assay: Spatial (19280 features, 0 variable features)
1 layer present: countsWe can now see that our object size significantly decreased.
Some QC plots
We now cleaned the genes that we are not interested in. We can now analyse the distributions of these genes with the commonly used histograms, but also spatially.
# check within the object more in details, we see the metadata
head(object@meta.data) orig.ident nCount_Spatial nFeature_Spatial
AAACAGAGCGACTCCT-1 SeuratProject 5651 2586
AAACAGTGTTCCTGGG-1 SeuratProject 8249 3613
AAACCGGGTAGGTACC-1 SeuratProject 10590 4006
AAACCGTTCGTCCAGG-1 SeuratProject 7285 3274
AAACCTCATGAAGTTG-1 SeuratProject 6727 2861
AAACGAAGAACATACC-1 SeuratProject 1756 1106nFeature_Spatial= number of unique genesnCount_Spatial= total number of counts (a genes can be present multiple times)
p1 <- ggplot(object@meta.data, aes(nCount_Spatial)) + geom_histogram()
p2 <- MapFeatures(object = object, features = "nCount_Spatial", slot = "counts", pt_size = 1)
p <- p1 | p2
p
We could decide to remove all spots that have less than a certain number of counts since the section we analyse could have holes or empty regions.
We would do it by running:
# Define thresholds for each subset object (adjust the values as needed)
threshold <- 2000
object$to_keep <- factor(ifelse(object$nCount_Spatial > threshold, "keep", "drop"))
sum(object$to_keep == "drop")[1] 164p <- MapLabels(
object = object,
column_name = "to_keep",
pt_size = 1.5
)
p
# spots_to_keep <- names(object$nCount_Spatial[object$nCount_Spatial > threshold])
# object2 <- SubsetSTData(object, spots = spots_to_keep)Before we continue with the processing, we just need to make sure that no spots contains exactly 0 count. That would produces NaN in some function, then breaking others.
Check:
sum(object$nCount_Spatial == 0) # check that no spots have 0 counts[1] 0all(object$nCount_Spatial == colSums(LayerData(object, "counts"))) # check that the metadata is the same data as the layer counts[1] TRUEProcessing
The data collected by the array and sequenced after library preparation are prone to biological and technical variation. This raw data is called “counts” in a Seurat object.
Before we can compare spots and genes to each other, we apply two corrections:
NormalizeData()works per spot. Spots capture different amounts of RNA, so a gene can look highly expressed simply because that spot yielded more molecules overall. Normalization puts every spot on a comparable scale.ScaleData()works per gene. It converts each gene to a z-score, so that a highly expressed gene does not dominate the PCA purely because its values are large.
Note the order: scaling is applied to the normalized values, not to the raw counts.
We checked before by printing the outlook of the object, that we have only the raw counts from after the spaceranger run.
Now, we will need to go through multiple steps:
1. Normalization
Normalize the counts (the amount of RNA and UMIs differ between spots (cells in case of single-cell RNA-seq)). The normalization is done per spot and is a log1p of a per spot scaled count.
Formula:
normalized count = log1p(count / total_counts × 10000)
object_norm <- NormalizeData(object = object)gene <- "MYH6"
p1 <- ggplot(data.frame(x = LayerData(object, "counts")[gene, ]), aes(x)) +
geom_histogram(bins = 50) +
labs(title = paste(gene, "- raw counts"), x = "count per spot")
p2 <- ggplot(data.frame(x = LayerData(object_norm, "data")[gene, ]), aes(x)) +
geom_histogram(bins = 50) +
labs(title = paste(gene, "- normalized"), x = "log-normalized value")
p3 <- MapFeatures(
object = object,
slot = "counts", # Here we use the non-normalized & unscaled layer
features = gene,
pt_size = 1
)
p4 <- MapFeatures(
object = object_norm,
slot = "data", # Here we use the normalized and scaled layer
features = gene,
pt_size = 1
)
p <- (p1 | p3) / (p2 | p4)
p
2. Find variable features
Some of the downstream processing are relatively computing costly, so we gain from restricting the gene space to only the genes that are the most variable across spots. Excluding the invariable genes help for:
- Improved signal-to-noise ratio
- Reduced computational requirements
Using the Seurat suite, we can do that easily by using the function FindVariableFeatures(). That function fits a line to the relationship between the log(variance) and the log(mean) using the local polynomial regression (loess) per gene.
# Find the top 3000 variable features
object_norm <- FindVariableFeatures(
object_norm,
selection.method = "vst",
nfeatures = 3000
)
# Generate the scatter plot (Black = non-variable, Red = top 3000 variable)
p1 <- VariableFeaturePlot(object_norm, cols = c("black", "red"))
# Add gene labels for the top 10 genes
p_final <- LabelPoints(
plot = p1,
points = head(VariableFeatures(object_norm), 10),
repel = TRUE,
xnudge = 0,
ynudge = 0
)
suppressWarnings(print(p_final))
3. Scaling
The scaling correspond to a Z-score transformation which is performed on the normalized count (after normalization step), and is applied per gene (and not per spot).
Formula: Z = (x - mu) / sigma
- x: normalized count
- mu: mean expression of the gene we normalize
- sigma: standard deviation of the gene we normalize
Scaled data have a mean of 0 and a variance of 1. That means that some genes have positive values and some are negative. This scaled data is used mainly when we want to further reduce dimension using PCA.
# let's assign the scaled data again to our variable called object. We will then have all the layers on the same object.
object <- ScaleData(object = object_norm)Now we can inspect the object again:
objectAn object of class Seurat
19280 features across 2234 samples within 1 assay
Active assay: Spatial (19280 features, 3000 variable features)
3 layers present: counts, data, scale.dataThat is all the transformation that we will perform on our data. Let’s visualize how different the same gene can look on the different layers;
p1 <- MapFeatures(object = object,
features = "MYH6",
slot = "counts") &
labs(title = "counts")
p2 <- MapFeatures(object = object,
features = "MYH6",
slot = "scale.data") &
labs(title = "scale.data")
p3 <- MapFeatures(object = object,
features = "MYH6",
slot = "data") &
labs(title = "data")
p <- p1 | p2 | p3
p
These 3 steps allowed us to correct for sequencing depths difference in genes and in spots. We found the genes that have the most amount of variance in the data (the genes that do not change in space in our section are not the most interesting), and scaled the data for PCA.
Clustering
Our data is now ready for downstream analysis. While we could explore individual genes of interest, a state-of-the-art approach borrowed from single-cell RNA-seq is cluster analysis. The primary goal of clustering is to group data points based on overall similarity. In spatial transcriptomics, we group spots that share similar gene expression profiles, which helps differentiate spots corresponding to distinct anatomical structures or cellular microenvironments.
It is again a multi-step processes:
- Principal Component Analysis (PCA)
- Find the neighbors & Build graph
- Find the clusters
- Visualize the clusters on a UMAP embedding
- Visualize the clusters on the tissue section
- Differentially Expressed Genes (DEG)
1. Principal Component Analysis (PCA)
PCA is a data-simplification technique. Instead of tracking 3000 individual genes across every spot, PCA groups together genes that behave similarly, genes telling the same story. This condenses 3000 separate measurements into a few master trends (Principal Components) that still capture the main biological differences in our sample.
The RunPCA() function from Seurat allows us to do this easily. The function will use the scale.data layer under the hood, and use only the top 3000 genes.
object <- RunPCA(object)p <- ElbowPlot(object = object,
ndims = 50,
reduction = "pca",
plot_type = "variance"
)
p
2. Find the neighbors & Build graph
Calculate distances between spots in PCA space to construct a Shared Nearest Neighbor (SNN) graph. In that example, we decide to take the first 20 principal components.
object <- FindNeighbors(object = object, reduction = "pca", dims = 1:20)3. Find the clusters
Partition the SNN graph into distinct groups using community detection algorithms (e.g., Louvain or Leiden). You see that there is a resolution parameter for the FindClusters() function.
object <- FindClusters(object = object,
algorithm = 1, #1 = original Louvain algorithm
resolution = 1) #higher number gets more clustersModularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 2234
Number of edges: 88122
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8059
Number of communities: 10
Elapsed time: 0 seconds4. Visualize the clusters on a UMAP embedding
Project the clusters onto a 2D layout to inspect global relationships. Here we want to visualize each of our spots as clusters that we just computed. So we take the first 20 principal components from our PCA analysis, and project them in a 2D space using the Uniform Manifold Approximation and Projection (UMAP), “a popular data science tool used for machine learning and data visualization. It shrinks complex information with many variables down into a simple 2D or 3D picture.”
object <- RunUMAP(object = object,
reduction = "pca",
dims = 1:20,
seed.use = 42
)Once UMAP is ran, we can visualize the results in a scatter plot. We need to specify the metadata column that contains the label of our clusters. In that case, we use the column called seurat_clusters.
p1 <- DimPlot(object = object,
reduction = "umap",
group.by = "seurat_clusters") &
labs(title = "UMAP - each point corresponds to a spot")
p1
5. Visualize the clusters on the tissue section
Overlay cluster colors directly onto the histology image to verify anatomical alignment.
p2 <- MapLabels(object = object,
column_name = "seurat_clusters",
pt_size = 1.5)
p2
It might be easier to visualize them side-by-side
p1 <- DimPlot(object = object,
reduction = "umap",
group.by = "seurat_clusters")
p2 <- MapLabels(object = object,
column_name = "seurat_clusters",
pt_size = 1)
p <- p1 | p2
p
6. Differentially Expressed Genes (DEG)
We would like now to identify which genes are the most specific to each of the clusters we identified. To do so we use an analysis method called Differentially Expressed Genes or DEG.
We rely mainly on the FindAllMarkers() function from the Seurat package. If you check the function documentation online, you can find out about the different test that one can choose from. We will use the default Wilcoxon test. The function accounts for the multiple testing. After running the test, we will explore quickly its results.
# Set identity to the clustering column corresponding to the current resolution
object <- SetIdent(object, value = "seurat_clusters")
# Differential expression analysis
detable <- FindAllMarkers(object,
only.pos = TRUE,
max.cells.per.ident = 300,
logfc.threshold = 0.1,
assay = "Spatial",
test.use = "wilcox",
min.pct = 0.05)
# Optionally, print a glimpse of the markers table
print(head(detable)) p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
MYL2 2.768468e-49 0.9670315 1 0.949 5.337607e-45 0 MYL2
TPM1 5.060174e-45 0.4885766 1 0.998 9.756016e-41 0 TPM1
CKM 1.978728e-39 0.6731430 1 0.949 3.814988e-35 0 CKM
MYL3 2.991804e-37 0.6973320 1 0.973 5.768198e-33 0 MYL3
TNNC1 1.139420e-36 0.4573792 1 0.992 2.196802e-32 0 TNNC1
ACTC1 1.748147e-36 0.4189044 1 0.995 3.370427e-32 0 ACTC1We can see that the results is a dataframe containing several metrics in the columns such as the average log fold change (avg_log2FC) or the adjusted p-values (p_val_adj), but also the name of the genes and the clusters it belongs to. We expect to recover known anatomical markers since the test is supposed to pull out which genes are the most specific to our clusters.
We hope to find some of these: - MYH6/NPPA atria, - MYH7/IRX4 ventricle, - TNNT2 pan-cardiomyocyte, - COL1A1/POSTN valve mesenchyme, - PECAM1 endothelium
gene_to_look_up <- "NPPA"
detable |> filter(gene == gene_to_look_up) p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
NPPA 5.448989e-66 2.3398905 1 0.997 1.050565e-61 3 NPPA
NPPA1 7.498068e-52 2.2608631 1 0.997 1.445627e-47 7 NPPA
NPPA2 7.471038e-06 0.8491217 1 0.997 1.440416e-01 9 NPPALet’s find out more:
dim(detable)[1] 10686 7detable |> filter(gene == "MYH6") p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
MYH6 3.544970e-63 2.866364 1.000 0.776 6.834702e-59 3 MYH6
MYH61 7.426592e-54 3.031206 0.984 0.784 1.431847e-49 7 MYH6
MYH62 4.528977e-12 1.554290 0.902 0.792 8.731867e-08 9 MYH6We can see that the dimension of the table is quite large, which means that some genes appear in multiple clusters. We also see that when checking all the rows that contain the gene MYH6.
We will see in a second how we can visualize these results easily.
But first, we want to filter these results and add some more metrics
# Filter and add extra columns
detable <- detable[detable$p_val_adj < 0.05, ] #keep only the rows that have an adjusted p-val smaller than 0.05
detable$pct.diff <- detable$pct.1 - detable$pct.2 # add a column called "pct.diff" which is the difference in percentage 1 (%age of spots expressing that gene in cluster of interest) and percentage 2 (%age of spots expressing that gene in all other clusters). We can think that the largest the pct.diff, the more specific the gene is to that cluster.
detable$log.pct.diff <- log2((detable$pct.1 * 99 + 1) / (detable$pct.2 * 99 + 1)) #calculate the log of the pct.diffThen we can select the top5 markers per clusters
# Select top genes by p_val_adj, pct.diff and avg_log2FC
top5 <- detable |>
group_by(cluster) |>
slice_min(p_val_adj, n = 20, with_ties = FALSE) |>
slice_max(pct.diff, n = 10, with_ties = FALSE) |>
slice_max(avg_log2FC, n = 5, with_ties = FALSE) |>
ungroup()
top5# A tibble: 50 × 9
p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene pct.diff log.pct.diff
<dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr> <dbl> <dbl>
1 1.76e-21 2.10 0.505 0.134 3.39e-17 0 FABP4 0.371 1.84
2 2.51e-31 1.32 0.908 0.503 4.83e-27 0 PDK4 0.405 0.839
3 6.96e-21 1.04 0.847 0.443 1.34e-16 0 MASP1 0.404 0.920
4 6.98e-24 0.855 0.994 0.734 1.35e-19 0 LBH 0.26 0.432
5 2.40e-25 0.744 1 0.797 4.63e-21 0 COX6… 0.203 0.324
6 1.43e-47 1.83 0.935 0.416 2.76e-43 1 HEY2 0.519 1.15
7 4.24e-42 1.59 0.845 0.402 8.17e-38 1 MT3 0.443 1.05
8 1.21e-30 1.25 0.906 0.576 2.33e-26 1 RASD1 0.33 0.644
9 1.58e-33 0.967 0.978 0.735 3.05e-29 1 LGAL… 0.243 0.407
10 2.17e-30 0.798 0.993 0.816 4.18e-26 1 GJA1 0.177 0.280
# ℹ 40 more rowsAnd finally visualize it
# Generate a dot plot for the selected genes with a title indicating the grid info
plts <- DotPlot(object, features = unique(top5 |> pull(gene)), assay = "Spatial") +
coord_flip() +
theme(axis.text.y = element_text(angle = 0, hjust = 1, vjust = 0.5)) +
scale_color_gradientn(colours = RColorBrewer::brewer.pal(n = 11, name = "RdBu") |> rev()) +
scale_x_discrete(limits = rev) +
scale_y_discrete() +
labs(x = "", y = "", title = "Dot Plot of spatial clusters") +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5))
plts
Comparing our clusters with the published atlas
Our clustering came from one section, ~2,200 spots, default parameters, and about 20 minutes of work. The published atlas clustered 25,208 spots from 17 sections and resolved 23 spatial clusters covering the major structural components of the developing heart. Comparing the two is a good way to calibrate what a quick analysis can and cannot recover.
Below is a small set of markers for structures the atlas resolved. This is a teaching subset — the complete per-cluster marker tables are in Supplementary Table 1 of the paper.
reference_markers <- list(
"Atrial myocardium" = c("MYH6", "NR2F2"),
"Ventricular myocardium" = c("MYH7", "IRX4", "MYL2"),
"Trabecular myocardium" = c("NPPA", "BMP10"),
"Cardiomyocyte (general)" = c("TNNT2", "ACTC1", "TTN"),
"Endocardium / endothelium"= c("PECAM1", "NPR3", "CDH5"),
"Valve mesenchyme" = c("POSTN", "COL1A1", "PENK"),
"Epicardium" = c("WT1", "UPK3B", "TBX18"),
"Conduction system" = c("SHOX2", "TBX3", "HCN4"),
"Neural / autonomic" = c("PRPH", "STMN2", "SOX10"),
"Proliferating" = c("MKI67", "TOP2A")
)
# keep only what survived our filtering
reference_markers <- lapply(reference_markers, intersect, y = rownames(object))
reference_markers <- reference_markers[lengths(reference_markers) > 0]
# gene → structure lookup, de-duplicated
g2s <- stack(reference_markers) |>
setNames(c("gene", "structure")) |>
dplyr::filter(!duplicated(gene))
p <- DotPlot(object, features = as.character(g2s$gene), assay = "Spatial")
d <- p$data
d$structure <- factor(g2s$structure[match(d$features.plot, g2s$gene)],
levels = names(reference_markers))
ggplot(d, aes(x = id, y = features.plot,
size = pct.exp, colour = avg.exp.scaled)) +
geom_point() +
facet_grid(structure ~ ., scales = "free_y", space = "free_y") +
scale_colour_gradientn(colours = rev(RColorBrewer::brewer.pal(11, "RdBu")),
name = "Avg. expression") +
scale_size(range = c(0, 6), name = "% expressed") +
labs(x = "cluster", y = NULL,
title = "Our clusters against published structural markers") +
theme_bw(base_size = 11) +
theme(
strip.text.y.right = element_text(angle = 0, hjust = 0, size = 9),
strip.background = element_rect(fill = "grey95", colour = NA),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 9),
panel.grid.minor = element_blank()
)
Passing a named list to DotPlot() groups the genes by structure, so each block of columns is one candidate identity. A cluster with a strong, isolated block is straightforward to name. A cluster with moderate signal in several blocks is probably a mixture — which is expected, since a 55 µm spot rarely contains only one cell type.
Now try annotating. Fill in your best guess for each cluster:
# Edit the labels, then run. Use NA where you are unsure.
annotation <- c(
"0" = "?",
"1" = "?",
"2" = "?",
"3" = "?",
"4" = "?",
"5" = "?",
"6" = "?",
"7" = "?",
"8" = "?",
"9" = "?"
)
object$manual_annotation <- unname(annotation[as.character(object$seurat_clusters)])
MapLabels(object, column_name = "manual_annotation", pt_size = 1.5)
Why our result is not the published one, and that is fine
Three reasons, all worth internalising.
Statistical power. Rare structures need spots. The conduction system occupies a tiny fraction of any section, so with ~2,200 spots there may be only a handful of spots covering the sinoatrial node — not enough to form their own cluster at any resolution.
Sampling. Whether a structure appears at all depends on the cutting plane. The atlas selected 17 sections that each contained at least three chambers; ours is whatever we were given.
Resolution and parameters. We used resolution = 1 and 20 principal components without examining alternatives. Try resolution = 0.4 and 1.5 and watch structures merge and split. There is no correct value — there is a value appropriate to the question.
None of this makes our clustering wrong. It makes it coarse, and knowing which of your results are robust to these choices is most of the skill in this kind of analysis.
Now we will do the same but with 3 sections
Researchers usually have specific questions which require the joint analysis of multiple sections simultaneously. Here we will add 2 extra heart sections of 2 different time points. We will repeat all the steps we did above, this time less explanation to read, and you can run quicker the code that allows you to do some temporal exploration.
1. Build infoTable
set.seed(1)
# Join metadata on sample_id rather than trusting that file order matches CSV
# row order — otherwise sections can silently get the wrong age.
metadata <- read.csv("./data/metadata.csv")
samples <- Sys.glob("./data/spaceranger/*/filtered_feature_bc_matrix.h5")
imgs <- Sys.glob("./data/spaceranger/*/spatial/tissue_hires_image.png")
spotfiles <- Sys.glob("./data/spaceranger/*/spatial/tissue_positions_list.csv")
json <- Sys.glob("./data/spaceranger/*/spatial/scalefactors_json.json")
stopifnot(length(samples) == length(imgs),
length(samples) == length(spotfiles),
length(samples) == length(json))
infoTable <- tibble(samples, imgs, spotfiles, json) |>
mutate(sample_id = basename(dirname(samples)),
slide_id = sample_id) |>
left_join(metadata, by = "sample_id") |>
arrange(pcw_age)
stopifnot(!any(is.na(infoTable$pcw_age)))
print(infoTable) # eyeball before loading# A tibble: 3 × 10
samples imgs spotfiles json sample_id slide_id organ organism pcw_age sex
<chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <int> <chr>
1 ./data/… ./da… ./data/s… ./da… V10A13-1… V10A13-… heart human 6 male
2 ./data/… ./da… ./data/s… ./da… V10A13-1… V10A13-… heart human 9 fema…
3 ./data/… ./da… ./data/s… ./da… V10F24-1… V10F24-… heart human 10 male Here you see in the column pcw_age (post conceptional week), that we have 3 sections of 3 distinct ages. We also know the sex of the embryos, let’s figure it out first!
2. Load data
This time we will call our object se, short for Seurat object. Note that now we are loading all the information of the 3 sections into a single object.
se <- ReadVisiumData(infoTable)
se <- LoadImages(se, image_height = 400)
se$sample_id <- factor(se$sample_id, levels = unique(se$sample_id))ImagePlot(se, ncol = 3)
3. Gene filtering
We repeat the 2 filtering steps we did earlier. Removing the MT, HB, RP and MALAT1, and annotating the genes in order to keep only the protein coding genes.
But first, let’s take advantage of some transcript to identify the sex of each section:
sex_genes <- c("XIST", "RPS4Y1", "DDX3Y", "UTY", "EIF1AY", "KDM5D")
FetchData(se, layer = "counts", vars = intersect(sex_genes, rownames(se))) |>
mutate(sample_id = se$sample_id) |>
group_by(sample_id) |>
summarise(across(everything(), mean))# A tibble: 3 × 7
sample_id XIST RPS4Y1 DDX3Y UTY EIF1AY KDM5D
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 V10A13-160_A1 0.00163 16.0 1.81 0.553 2.45 1.14
2 V10A13-157_C1 1.97 0 0 0 0 0
3 V10F24-105_C1 0 4.12 0.481 0.0868 1.09 0.226We used raw counts (before normalization), because we do the processing steps after filtering. Here the values are more of a qualitative measurement rather than quantitative. We want to know if the genes are there (male) or not (female).
infoTable |>
select(sample_id, pcw_age, sex)# A tibble: 3 × 3
sample_id pcw_age sex
<chr> <int> <chr>
1 V10A13-160_A1 6 male
2 V10A13-157_C1 9 female
3 V10F24-105_C1 10 male Now we can proceed with our filtering:
pattern <- paste0(
"^MT-",
"|^HB[ABDEGMQZ][0-9]?$",
"|^RP[LS][0-9]+[A-Z]?[0-9]?$|^RPL36AL$|^RPLP[0-9]$|^RPSA$",
"|^MALAT1$"
)
undesirable_genes <- setdiff(grep(pattern, rownames(se), value = TRUE), "RPL3L")
annot <- readRDS("./data/GRCh38-2020-A_gene_biotypes.rds")
protein_coding <- subset(annot, gene_type == "protein_coding")$gene_name
genes_to_keep <- setdiff(rownames(se), undesirable_genes) |> intersect(protein_coding)
se <- SubsetSTData(se, features = genes_to_keep)
# metadata is NOT refreshed by subsetting — recompute it
se$nCount_Spatial <- colSums(LayerData(se, "counts"))
se$nFeature_Spatial <- colSums(LayerData(se, "counts") > 0)
stopifnot(sum(se$nCount_Spatial == 0) == 0)
seAn object of class Seurat
19280 features across 4394 samples within 1 assay
Active assay: Spatial (19280 features, 0 variable features)
1 layer present: counts4. Per-section QC
p1 <- ggplot(se@meta.data, aes(nCount_Spatial, fill = sample_id)) +
geom_histogram(bins = 50) +
facet_wrap(~ sample_id, scales = "free_y") +
labs(title = "Library size per spot, by section")
p2 <- MapFeatures(se, features = "nCount_Spatial", pt_size = 1, ncol = 3, label_by = "sample_id")
p <- p1 / p2
p
5. Processing
se <- NormalizeData(se)
se <- FindVariableFeatures(se, selection.method = "vst", nfeatures = 3000)
se <- ScaleData(se)
se <- RunPCA(se, npcs = 50)
ElbowPlot(se, ndims = 50, reduction = "pca", plot_type = "variance")
6. Batch correction
Three sections were captured on different slides, so part of what separates them is technical. Harmony adjusts the PCA embedding so that spots group by biology rather than by slide. Run the clustering below BOTH ways and compare: if clusters split perfectly by sample_id on “pca”, that split is batch.
se <- harmony::RunHarmony(se, group.by.vars = "sample_id",
reduction = "pca", reduction.save = "harmony")Initializing centroids7. Clustering: uncorrected versus corrected
We now cluster twice, once on the raw PCA and once on the Harmony-corrected embedding, and keep both results so we can compare them side by side. Everything else — the number of dimensions, the algorithm, the resolution — stays identical, so any difference we see comes from the correction alone.
# 7.1 — clustering on the uncorrected PCA
se <- FindNeighbors(se, reduction = "pca", dims = 1:20)
se <- FindClusters(se, algorithm = 1, resolution = 0.8)Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 4394
Number of edges: 152102
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8819
Number of communities: 17
Elapsed time: 0 secondsse$clusters_pca <- se$seurat_clusters # save before it gets overwritten
se <- RunUMAP(se, reduction = "pca", dims = 1:20, seed.use = 42,
reduction.name = "umap_pca")
# 7.2 — same steps, this time on the Harmony embedding
se <- FindNeighbors(se, reduction = "harmony", dims = 1:20)
se <- FindClusters(se, algorithm = 1, resolution = 0.8)Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 4394
Number of edges: 163592
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8428
Number of communities: 16
Elapsed time: 0 secondsse$clusters_harmony <- se$seurat_clusters
se <- RunUMAP(se, reduction = "harmony", dims = 1:20, seed.use = 42,
reduction.name = "umap_harmony")Note that we copy seurat_clusters into a new column each time. FindClusters() always writes to that same column, so without the copy the first result would be lost. The same applies to RunUMAP(), which is why we give each embedding its own reduction.name.
Are the spots grouping by biology or by section?
The clearest diagnostic is to colour the UMAP by section rather than by cluster. If spots separate cleanly by sample_id, the dominant signal is where the spot came from, not what tissue it is.
p1 <- DimPlot(se, reduction = "umap_pca", group.by = "sample_id") + ggtitle("Before (PCA)")
p2 <- DimPlot(se, reduction = "umap_harmony", group.by = "sample_id") + ggtitle("After (Harmony)")
p <- p1 | p2
p
On the uncorrected embedding, the three sections form largely separate territories. After Harmony they overlap, meaning spots from different sections now sit next to each other when they are transcriptionally similar.
We can also ask how the two clusterings relate. Each row is an uncorrected cluster, each column a corrected one:
table(se$clusters_pca, se$clusters_harmony)
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
0 21 374 0 2 0 57 3 1 2 0 47 0 0 11 12 0
1 2 87 188 186 0 8 1 0 1 0 0 0 0 0 2 3
2 442 8 3 0 0 1 0 0 5 0 2 0 0 0 10 0
3 0 0 0 0 203 0 122 0 0 0 0 0 54 0 0 0
4 324 17 1 0 0 15 0 0 1 0 6 0 0 2 2 0
5 61 4 8 102 0 93 0 12 5 0 3 0 0 13 2 9
6 11 3 193 27 1 31 0 0 1 0 0 0 0 0 0 30
7 1 1 0 0 0 1 2 189 1 0 36 19 1 42 0 0
8 0 0 0 0 100 30 124 3 0 0 3 0 7 0 0 0
9 0 165 0 83 0 2 0 0 9 0 0 0 0 1 0 0
10 0 0 0 0 0 0 0 12 0 151 0 0 0 0 0 0
11 1 0 0 2 0 96 0 2 2 0 0 45 0 0 0 0
12 0 0 66 11 0 0 0 0 58 0 0 0 0 0 2 4
13 0 0 0 0 38 0 41 0 0 0 0 0 18 0 0 0
14 0 0 0 5 0 0 10 0 80 0 0 0 0 0 0 0
15 0 0 0 0 0 3 0 36 0 0 0 24 0 0 0 0
16 1 1 0 0 0 1 0 1 0 0 0 0 0 0 28 0A cluster whose spots spread across many columns was split or reassigned by the correction. A tight diagonal block means that cluster was robust either way.
Back to the tissue
The UMAP is a summary; the section is the ground truth. A good correction should make the same anatomical structure receive the same cluster label in all three sections.
MapLabels(se, column_name = "clusters_pca", pt_size = 1, ncol = 3) +
plot_layout(guides = "collect") &
guides(fill = guide_legend(override.aes = list(size = 3),
ncol = 2)) &
theme(legend.position = "right")
MapLabels(se, column_name = "clusters_harmony", pt_size = 1, ncol = 3) +
plot_layout(guides = "collect") &
guides(fill = guide_legend(override.aes = list(size = 3),
ncol = 2)) &
theme(legend.position = "right")
An important limitation of this dataset
We have one section per age, so sample_id and pcw_age carry exactly the same information:
table(se$sample_id, se$pcw_age)
6 9 10
V10A13-160_A1 615 0 0
V10A13-157_C1 0 1545 0
V10F24-105_C1 0 0 2234That means batch and biology are perfectly confounded here. When we tell Harmony to remove differences between sections, we are also telling it to remove differences between developmental stages — it has no way to tell them apart. The corrected clustering is therefore useful for finding shared anatomy across sections, but it cannot be used as evidence that a structure is or is not changing with age.
The published atlas avoided this by using at least two near-consecutive sections per heart as technical replicates, across 16 hearts spanning PCW 6 to 12. With replicates, the section-to-section variation can be estimated separately from the age effect.
Batch correction is as much judgement as method, and it depends on the question being asked. Over-correct and you erase real differences between samples; under-correct and the technical signal drives your clusters. Harmony exposes several parameters (notably theta, which controls how aggressively it mixes groups) that are worth reading about the day you need this on your own data.
8. Cluster markers
se <- SetIdent(se, value = "clusters_harmony")
detable <- FindAllMarkers(se,
only.pos = TRUE,
max.cells.per.ident = 300,
logfc.threshold = 0.1,
assay = "Spatial",
test.use = "wilcox",
min.pct = 0.05)
detable <- detable |>
filter(p_val_adj < 0.05) |>
mutate(pct.diff = pct.1 - pct.2)
# slice_* rather than top_n: p_val_adj has many exact ties at 0
top5 <- detable |>
group_by(cluster) |>
slice_min(p_val_adj, n = 20, with_ties = FALSE) |>
slice_max(pct.diff, n = 10, with_ties = FALSE) |>
slice_max(avg_log2FC, n = 5, with_ties = FALSE) |>
ungroup()
DotPlot(se, features = unique(top5$gene), assay = "Spatial") +
coord_flip() +
scale_color_gradientn(colours = rev(RColorBrewer::brewer.pal(11, "RdBu"))) +
scale_x_discrete(limits = rev) +
labs(x = "", y = "", title = "Cluster markers across 3 sections")
Here, each cluster shares spots coming from our 3 sections, which means that the markers we find in our DEG analysis, are markers that best describe these populations. This can be very useful to understand which markers drive these spots to cluster together and we expect to find markers that are significantly present in all sections analysed.
We can also plot a single marker on all the sections simultaneously:
gene <- "MYH6"
p <- MapFeatures(
object = se,
features = gene,
ncol = 3
)
p
9. Checkpoint
When analysis get long, it is common practice to save the status of our object. Many results such as plots and tables can be saved along the way too. Here we will simply see how it is done, but in order to not clutter the server at our disposal, we will not save it.
# dir.create("./checkpoints", showWarnings = FALSE) # create a directory to save results
# saveRDS(se, "./checkpoints/3sections_clustered.rds") # save the 3 sections object
# saveRDS(detable, "./checkpoints/3sections_markers.rds") # save the detable
# se <- readRDS("./checkpoints/3sections_clustered.rds") # load the saved object, when we want to continue working on the projectAdvanced analysis
10. Temporal DEG
We just did a DEG analysis on the 3 sections together, but an information we do not know yet is what markers are specific per clusters AND time point. For that we need to subset each clusters per time points individually and repeat FindAllMarkers() function to look at each of the time point separately.
Let’s try it out:
1. Select cluster of choice to split per time point
# pick a cluster that is well represented in all age
table(se$clusters_harmony, se$pcw_age)
6 9 10
0 62 338 464
1 4 189 467
2 74 194 191
3 113 117 188
4 38 110 194
5 96 174 68
6 43 130 130
7 48 35 173
8 62 88 15
9 0 0 151
10 3 24 70
11 24 45 19
12 18 17 45
13 13 27 29
14 4 30 24
15 13 27 6cl <- "3"
sub <- SubsetSTData(se, spots = colnames(se)[which(se$clusters_harmony == cl)]) # subset the se object to only the spots from our cluster of choice
sub$pcw_age <- factor(sub$pcw_age, levels = sort(unique(as.numeric(se$pcw_age)))) # to get ages in the right order
Idents(sub) <- "pcw_age"
res <- FindAllMarkers(sub,
only.pos = TRUE,
max.cells.per.ident = 300,
logfc.threshold = 0.1,
assay = "Spatial",
test.use = "wilcox",
min.pct = 0.05) |>
filter(p_val_adj < 0.05)
res <- as_tibble(res)
res_top <- res |> dplyr::rename(age = cluster) |> group_by(age) |> slice_max(avg_log2FC, n = 10) |> ungroup()By doing this, we isolate all the spots corresponding to a cluster of interest and ask “Within this cluster, which are the markers that differ the most across age/section?”. The way it is calculated is, taking the spots of the clusters of choice of 1 age (e.g. pcw 6) and compare it to the spots of age 2 and 3 together (e.g. pcw 9 and 10 pooled together). This step is repeated for all the age categories (pcw 9 against 6 & 10, and then pcw 10 against 6 and 9 together).
2. Visualisation
top2_per_age <- res_top |>
group_by(age) |>
slice_min(p_val_adj, n = 2, with_ties = FALSE) |>
ungroup()
top_2_markers <- unique(top2_per_age$gene)
VlnPlot(sub, features = top_2_markers, group.by = "pcw_age", ncol = 2)
# Plot top2 markers per category
DotPlot(sub, features = top_2_markers, assay = "Spatial") +
coord_flip() +
scale_color_gradientn(colours = rev(RColorBrewer::brewer.pal(11, "RdBu"))) +
scale_x_discrete(limits = rev) +
labs(x = "", y = "", title = "Age-associated genes within a cluster")
But maybe, the actual spatial plot might be even easier. Let’s plot it on the 3 sections.
MapFeatures(
object = se,
features = top_2_markers,
ncol = 3,
pt_size = 0.8
)
11. How to take advantage of the spatial information bioinformatically
We covered many steps until here. Many of the steps that are covered in most Visium data analysis. But it is important to note that most of these steps are borrowed from single-cell RNA seq. But we can do way more with spatial data by taking advantage of the spots coordinates as an actual spatial information, instead of just as a barcode.
A question that often comes out is:
“Are the genes we called variable actually spatially structured?”
A gene might have an intrinsic high variability across spots, but its expression be completely scattered across the tissue section. We can rely on some spatial statistics to answer that question.
# select the top 500 variable features
test_features <- VariableFeatures(se)[1:500]
# evaluate their spatial autocorrelation
sp_cor <- CorSpatialFeatures(object = se, features = test_features)
# Check the list on section 1
# genes with the most spatially structured expression
head(sp_cor[[1]], 20) # A tibble: 20 × 2
gene cor
<chr> <dbl>
1 NPPA 0.970
2 MYH6 0.958
3 FBLN5 0.915
4 ELN 0.895
5 COL1A1 0.893
6 SULT1E1 0.880
7 OGN 0.866
8 BMP10 0.850
9 CRABP2 0.847
10 SLN 0.847
11 PAM 0.844
12 COL1A2 0.843
13 PI15 0.835
14 GPC3 0.834
15 HTRA1 0.829
16 CXCL12 0.829
17 ANKRD1 0.815
18 ACTA1 0.812
19 MFAP4 0.810
20 IRX3 0.805# genes with the least spatially structured expression
tail(sp_cor[[1]], 20)# A tibble: 20 × 2
gene cor
<chr> <dbl>
1 ONECUT2 -0.00546
2 CELF3 -0.00951
3 CPA3 -0.0101
4 PHOX2B -0.0104
5 SIGLEC15 -0.0108
6 DUSP9 -0.0112
7 PRG4 -0.0113
8 NEFM -0.0115
9 MYO3B -0.0116
10 NECAB2 -0.0150
11 CXCL13 -0.0169
12 DBH -0.0197
13 TREML1 -0.0259
14 TAC1 -0.0270
15 NKG7 -0.0328
16 DEFA3 NaN
17 KRT13 NaN
18 RHEX NaN
19 CHL1 NaN
20 OPRM1 NaN Again, it is the best to visualize them spatially:
Section 1
section_number <- 1
positive_corr <- sp_cor[[section_number]]$gene[1:6]
MapFeatures(
object = se,
features = positive_corr,
section_number = section_number,
ncol = 3
)
negative_corr <- tail(subset(sp_cor[[section_number]], !is.nan(cor))$gene, 6)
MapFeatures(
object = se,
features = negative_corr,
section_number = section_number,
ncol = 3
)
Section 2
section_number <- 2
positive_corr <- sp_cor[[section_number]]$gene[1:6]
MapFeatures(
object = se,
features = positive_corr,
section_number = section_number,
ncol = 3
)
negative_corr <- tail(subset(sp_cor[[section_number]], !is.nan(cor))$gene, 6)
MapFeatures(
object = se,
features = negative_corr,
section_number = section_number,
ncol = 3
)
Section 3
section_number <- 3
positive_corr <- sp_cor[[section_number]]$gene[1:6]
MapFeatures(
object = se,
features = positive_corr,
section_number = section_number,
ncol = 3
)
negative_corr <- tail(subset(sp_cor[[section_number]], !is.nan(cor))$gene, 6)
MapFeatures(
object = se,
features = negative_corr,
section_number = section_number,
ncol = 3
)
Closing note
In this workshop we went through the most important pre-processing, processing, clustering steps for a single section and for a couple of sections together. We explored some quality metrics, tried to correct for batch effect using Harmony, and explored a way to identify temporal changes across clusters of interest. Finally we understood that this is just the start, and from there much more spatially aware analyses can be performed.
Due to lack of time, we have not covered deconvolution, niche analysis, pathways, spatially aware clustering, spatial interactions, pseudotime analysis, and more.
For further explorations, here are some suggested resources:
- semla documentation
- Seurat spatial vignette
- Seurat guided clustering
- Seurat integration & label transfer
- Single-cell QC & Clustering
- NBIS detailed slides
- Deconvolution with Stereoscope
- Spatially aware clustering with BANKSY
- Spatial interactions with SpatialDM
sessionInfo()R version 4.3.3 (2024-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 24.04.4 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
[4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
time zone: Etc/UTC
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] future_1.75.0 harmony_2.0.5 Rcpp_1.1.2 RColorBrewer_1.1-3
[5] patchwork_1.3.2 semla_1.4.2 Seurat_5.5.1 SeuratObject_5.4.0
[9] sp_2.2-3 lubridate_1.9.5 forcats_1.0.1 stringr_1.6.0
[13] dplyr_1.2.1 purrr_1.2.2 readr_2.2.0 tidyr_1.3.2
[17] tibble_3.3.1 ggplot2_4.0.3 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] rstudioapi_0.19.0 jsonlite_2.0.0 magrittr_2.0.5
[4] spatstat.utils_3.2-4 magick_2.9.1 farver_2.1.2
[7] rmarkdown_2.31 ragg_1.5.2 vctrs_0.7.3
[10] ROCR_1.0-12 spatstat.explore_3.8-2 htmltools_0.5.9
[13] sctransform_0.4.3 parallelly_1.48.0 KernSmooth_2.23-22
[16] htmlwidgets_1.6.4 ica_1.0-3 plyr_1.8.9
[19] plotly_4.12.1 zoo_1.9-0 igraph_2.3.3
[22] mime_0.13 lifecycle_1.0.5 pkgconfig_2.0.3
[25] Matrix_1.6-5 R6_2.6.1 fastmap_1.2.0
[28] fitdistrplus_1.2-6 shiny_1.14.0 digest_0.6.39
[31] tensor_1.5.1 RSpectra_0.16-2 irlba_2.3.7
[34] textshaping_1.0.5 labeling_0.4.3 progressr_1.0.0
[37] spatstat.sparse_3.2-0 timechange_0.4.0 httr_1.4.8
[40] polyclip_1.10-7 abind_1.4-8 compiler_4.3.3
[43] bit64_4.8.4 withr_3.0.3 S7_0.2.2
[46] fastDummies_1.7.6 MASS_7.3-60.0.1 tools_4.3.3
[49] lmtest_0.9-40 otel_0.2.0 httpuv_1.6.17
[52] future.apply_1.20.2 goftest_1.2-3 glue_1.8.1
[55] dbscan_1.2.6 nlme_3.1-164 promises_1.5.0
[58] grid_4.3.3 Rtsne_0.17 cluster_2.1.6
[61] reshape2_1.4.5 generics_0.1.4 hdf5r_1.3.12
[64] gtable_0.3.6 spatstat.data_3.1-9 tzdb_0.5.0
[67] data.table_1.18.6.1 hms_1.1.4 utf8_1.2.6
[70] spatstat.geom_3.8-2 RcppAnnoy_0.0.23 ggrepel_0.9.8
[73] RANN_2.6.3 pillar_1.11.1 spam_2.11-4
[76] RcppHNSW_0.7.0 later_1.4.8 splines_4.3.3
[79] lattice_0.22-5 bit_4.6.0 survival_3.5-8
[82] deldir_2.0-4 tidyselect_1.2.1 miniUI_0.1.2
[85] pbapply_1.7-4 knitr_1.51 gridExtra_2.3.1
[88] scattermore_1.2 RhpcBLASctl_0.23-42 xfun_0.60
[91] matrixStats_1.5.0 stringi_1.8.9 yaml_2.3.12
[94] evaluate_1.0.5 codetools_0.2-19 cli_3.6.6
[97] uwot_0.2.4 xtable_1.8-8 reticulate_1.46.0
[100] systemfonts_1.3.2 globals_0.19.1 spatstat.random_3.5-1
[103] zeallot_0.2.0 png_0.1-9 spatstat.univar_3.2-0
[106] parallel_4.3.3 dotCall64_1.2 listenv_1.0.0
[109] viridisLite_0.4.3 scales_1.4.0 ggridges_0.5.7
[112] rlang_1.3.0 cowplot_1.2.0 shinyjs_2.1.1