Spatial Transcriptomics with Visium

Analysing sections of the developing human heart

Raphaël Mauron

V4SDB Student Summer School · 2 September 2026

Welcome

Today

You are going to analyse spatial transcriptomics data from a human embryonic heart.

  • All the code is written for you
  • Your job is to read it, run it, and understand why each step exists
  • The biology is yours — you know hearts better than the code does

Note

Everything we do today is a simplified version of a real, published pipeline.

Plan for the session

Wednesday 2 September 2026 · 09:00 – 13:00

Time Topic
09:00 Introduction, setup, first look at the tissue
09:30 Filtering and quality control one section
10:15 Normalization, scaling, variable genes
10:40 Break
10:55 PCA, clustering, mapping clusters onto anatomy
11:35 Cluster markers, comparison with the published atlas
12:10 Three sections together, batch correction three sections
12:45 Temporal changes, spatial statistics, wrap-up
13:00 End

Parts 10 and 11 are written so you can finish them on your own.

Two rules for the day

  1. Ask immediately. If a chunk errors, do not silently fall behind — we fix it together.
  2. Run in order, top to bottom. Several steps overwrite objects. Jumping around produces errors that look mysterious but are not.

Red text in the console is often just a warning. Not everything red is broken.

Part 0 · What is Visium?

The problem

Bulk RNA-seq — what genes, but not where. Single-cell RNA-seq — what genes in which cell, but the tissue is dissociated and position is destroyed.

Spatial transcriptomics — keep the tissue intact, and record where each transcript came from. A poly-T oligo binds the poly-A tail in situ; the barcode underneath records the position.

How the array works

A glass slide carries a grid of spots:

  • ~5,000 spots per capture area
  • 55 µm diameter · 100 µm centre-to-centre

Every spot carries millions of capture oligos sharing a spatial barcode unique to that position.

Section onto the slide → permeabilise → mRNA diffuses down → poly-T captures poly-A → sequence.

The barcode says where. The read says which gene.

Two consequences that shape everything

1. A spot is a position, not a cell. At 55 µm it covers a handful of cells. Every measurement is a local mixture.

2. Poly-A capture is indiscriminate. Anything polyadenylated gets sequenced, whether or not you wanted it.

A large part of today is deciding what to keep.

Where the data comes from

Lázár, Mauron, Andrusivová et al.
Spatiotemporal gene expression and cellular dynamics of the developing human heart.
Nature Genetics 57, 2756–2771 (2025)

The full study

  • 16 hearts · 38 Visium sections
  • 69,114 spots · PCW 6–12
  • scRNA-seq: 76,991 cells
  • ISS: 150 transcripts

Today

  • 3 sections
  • PCW 6, 9, 10
  • ~2,200 spots for part 1

Interactive atlas → https://hdcaheart.serve.scilifelab.se

Tissue donated after elective medical abortion with written informed consent, Swedish ethical permit 2018/769-31.

Part 1 · One section

What Space Ranger gives us

File What it is
filtered_feature_bc_matrix.h5 the counts — genes × spots
tissue_hires_image.png the H&E image
tissue_positions_list.csv which pixel each barcode sits on
scalefactors_json.json how to scale image to spots 
object <- ReadVisiumData(infoTable)
object <- LoadImages(object)
object
An object of class Seurat 
36601 features across 2234 samples within 1 assay 
Active assay: Spatial (36601 features, 0 variable features)
 1 layer present: counts

How many features? how many spots? which assay? which layers exist so far?

Before anything else: look at the tissue

At PCW 10 you should find the compact and trabeculated myocardium, the atrial wall, the interventricular septum, the valve region, the chamber lumens and the epicardium.

Important

Every plot today puts coloured dots on this image. They only mean something if you know what is underneath.

First orientation: MYH6 vs MYH7

Atrial and ventricular myosin heavy chains — they should split the section in two.

Does the split match where you expected atria and ventricles? Find a region where neither is expressed — what is there instead?

We have ~36,000 genes. Humans have ~19,000.

Four families dominate the counts and vary for reasons unrelated to anatomy:

  • Hemoglobin — the heart is full of blood
  • Mitochondrial — 13 genes; tracks muscle density and RNA degradation
  • Ribosomal proteins — ~85 genes, move as a block, 10–25% of counts
  • MALAT1 — nuclear lncRNA; reads out permeabilisation more than biology

Together, about a third of all counts in this dataset.

The pattern we actually use

pattern <- paste0(
  "^MT-",                      # 13 mitochondrial genes
  "|^HB[ABDEGMQZ][0-9]?$",     # HBA1/2, HBB, HBD, HBE1, HBG1/2, HBM, HBQ1, HBZ
  "|^RP[LS][0-9]+[A-Z]?[0-9]?$|^RPL36AL$|^RPLP[0-9]$|^RPSA$",
  "|^MALAT1$"
)
undesirable_genes <- setdiff(grep(pattern, rownames(object), value = TRUE), "RPL3L")
length(undesirable_genes); head(undesirable_genes, 18)
[1] 110
 [1] "RPL22"   "RPL11"   "RPS8"    "RPL5"    "RPS27"   "RPS7"    "RPS27A" 
 [8] "RPL31"   "RPL37A"  "RPL32"   "RPL15"   "RPSA"    "RPL14"   "RPL29"  
[15] "RPL24"   "RPL22L1" "RPL39L"  "RPL35A"

Those $ anchors are load-bearing — a looser pattern silently swallowed eleven RPS6K kinases. Always check the other direction too.

Second filter: keep protein-coding only

Annotate the remainder against the GRCh38-2020-A reference and keep protein_coding. Project-dependent — studying immunoglobulin loci or lncRNAs? Different call.

A trap

Subsetting genes does not update the metadata. Recompute it, or your QC describes a matrix you no longer have.

object <- SubsetSTData(object, features = genes_to_keep_2)

object$nCount_Spatial   <- colSums(LayerData(object, "counts"))
object$nFeature_Spatial <- colSums(LayerData(object, "counts") > 0)

all(object$nCount_Spatial == colSums(LayerData(object, "counts")))  # must be TRUE

Quality control — but for spatial

In single-cell, low counts means a dying cell — filter it. In Visium a spot is a fixed position, so low counts usually means low cellularity: a lumen, cardiac jelly, a section edge.

Tip

So we keep every spot — the low-count regions are anatomy.

Normalization and scaling — different axes

NormalizeData() · per spot

log1p(count / total × 10000)

ScaleData() · per gene

(x − µ) / σ — a z-score

Same gene, three layers. Scaling is applied to the normalized values.

Variable features

3,000 genes instead of 19,000 — better signal-to-noise, much faster.

vst fits log(variance) against log(mean) by loess, then standardises each gene against its expected variance:

\[z_{ij}=\frac{x_{ij}-\bar{x}_i}{\sigma_i^{\text{exp}}}\]

Genes are ranked by the variance of those standardised values.

High variance ≠ spatially organised. We return to this at the very end.

Clustering, in five steps

PCA — 3,000 genes condensed into a few dozen “master trends”

Neighbours — build a shared-nearest-neighbour graph

Clusters — Louvain community detection on that graph

UMAP — a 2D picture of that graph

Map it back onto the tissue ← the important one

The UMAP is not the result

Important

A cluster tracing a recognisable structure is believable. A cluster scattered across unrelated regions is usually a technical gradient.

Cluster markers

FindAllMarkers() — Wilcoxon test, each cluster against all others.

detable <- FindAllMarkers(object, only.pos = TRUE,
                          max.cells.per.ident = 300, ...)

We add pct.diff = pct.1 − pct.2: how specific a gene is to its cluster.

Do you recognise any genes?

Why does MYH6 appear in several clusters at once?

Comparing with the published atlas

Us

  • 1 section · ~2,200 spots
  • default parameters
  • 20 minutes

The atlas

  • 17 sections · 25,208 spots
  • 23 spatial clusters
  • deconvolved to 72 cell states

Not a fair fight — and that is the point. What can a quick analysis recover, and what does it miss?

Why we will not find 23 clusters

Power — the SA node is a handful of spots. No resolution will find it.

Sampling — whether a structure appears depends on the cutting plane.

Parameters — we picked resolution = 1 and 20 PCs without looking at alternatives.

Coarse is not wrong. Knowing which results survive these choices is most of the skill.

Part 2 · Three sections

Why more than one section?

  • Is your cluster reproducible, or a one-slide artefact?
  • Does a structure change between PCW 6 and PCW 10?
  • n = 1 tells you nothing about variability

Before we look

What should look the same across ages? What should differ biologically? What might differ for purely technical reasons?

Can we recover the sex from the data?

Mean raw counts per section, next to the reported metadata:

sample_id pcw_age sex XIST DDX3Y
V10A13-160_A1 6 male 0.00 1.81
V10A13-157_C1 9 female 1.97 0.00
V10F24-105_C1 10 male 0.00 0.48

Must be done before filtering — XIST is a lncRNA, and RPS4Y1 looks ribosomal to our regex.

Batch effect

Three sections, three slides, three capture events. Some of what separates them is technical.

Harmony adjusts the PCA embedding so spots group by biology rather than by slide.

se <- RunHarmony(se, group.by.vars = "sample_id")

We cluster both ways and compare. If clusters split perfectly by sample_id, that split was batch.

The diagnostic: colour by section, not by cluster

Separate territories → the dominant signal is which slide the spot came from. Overlapping → spots sit together when transcriptionally similar.

But: our design is confounded

One section per age, so sample_id and pcw_age carry identical information.

               
                   6    9   10
  V10A13-160_A1  615    0    0
  V10A13-157_C1    0 1545    0
  V10F24-105_C1    0    0 2234

Telling Harmony to remove differences between sections also tells it to remove differences between developmental stages. It cannot tell them apart.

With at least 2 sections per heart they stop being the same variable. You could correct on sample_id while pcw_age still varies within each corrected group — the slide effect goes, the age signal survives. Exactly why the atlas took near-consecutive replicate sections from all 16 hearts.

Over- and under-correction

Over-correct

You erase real differences between samples

Under-correct

Technical signal drives your clusters

There is no setting correct in general. There is a setting appropriate to your question, and you have to be able to defend it.

Part 3 · Going further

Does anything change with age?

Take one cluster, split it by age. Every spot comes from the same cluster, so the comparison is about age, not cluster identity. Genes ranked by pct.diff, not fold-change — specificity, not loudness.

age gene avg_log2FC pct.diff
6 KRT19 4.69 0.90
6 HMGA2 2.67 0.73
6 CLDN6 2.46 0.71
9 AHSP 2.96 0.51
9 HSPA1B 2.94 0.50
9 ALAS2 3.15 0.48
10 EIF1AY 2.23 0.48
10 TMEM176B 3.31 0.47
10 TMEM176A 3.43 0.36

Each age is tested against the other two pooled — and with n = 1 per age, every “temporal” gene is equally a “section” gene. Hypothesis-generating only.

Everything so far ignored the coordinates

Almost every step today came from single-cell RNA-seq. The coordinates were used only for plotting.

The question that needs space

“Are the genes we called variable actually spatially structured?”

A gene can vary enormously across spots and still be scattered randomly across the tissue.

sp_cor <- CorSpatialFeatures(se, features = VariableFeatures(se)[1:500])

For each gene: how similar is its expression between neighbouring spots?

Proof — same 500 genes, two ends of the ranking

Left two: highest spatial autocorrelation — coherent domains.

Right two: lowest — scattered, despite all four being in the top 500 variable genes.

Wrap-up

What we did

  • Loaded Visium data and plotted genes on tissue
  • Made and justified filtering decisions
  • Normalized, scaled, reduced, clustered
  • Mapped clusters back onto anatomy and named them
  • Compared against a published atlas
  • Scaled to three sections and confronted batch effects
  • Asked a temporal question, and found its limits
  • Used the coordinates as actual spatial information

What we skipped

  • Deconvolution — mapping single-cell states into spots
  • Niche analysis — which structures neighbour which
  • Spatially aware clustering — BANKSY and similar
  • Cell–cell interaction inference
  • Pathway and pseudotime analysis

All of these are in the paper, and the atlas viewer lets you explore the results without writing any code.

Where to go next

The atlas → https://hdcaheart.serve.scilifelab.se

The code → https://github.com/rmauron/HDCA_heart_dev

The one thing to take away

Every step was a choice

Which genes to remove. Whether to filter spots. How many PCs. What resolution. Whether to correct for batch, and how hard.

The software will never refuse. You have to be able to say why.

Thank you

Questions welcome — today, and afterwards.

Raphaël Mauron