16 MERFISH: Human Heart

Learning Objectives

Go through standard analysis steps for preprocessed MERFISH data
Understand how to visualise and interpret MERFISH data
Basic interpretation of markers and spatial patterns

For this example, we will be using the preprocessed MERFISH dataset of the human heart again. The dataset contains spatially resolved gene expression data for various cell types in the heart tissue. As there is no image associated with the data, this is basically a spatially resolved single-cell dataset. We will be using single cell functions in Seurat to analyse and visualise the data. The spatial coordinates are defined by dimensional reduction “spatial”, which we will use for visualisation instead of UMAP.

# Load the required libraries
library(Seurat)
library(paletteer)

# Load the preprocessed MERFISH dataset
merfish <- readRDS("data/human_heart_merfish/overall_merfish.rds")

# Visualise the structure of the heart in the MERFISH data
DimPlot(
  merfish,
  reduction = "spatial",
  group.by = "communities",
  cols = paletteer_d("ggthemes::Tableau_20")
)

As you can see, the data is already annotated with cluster identities and tissue types. We can use this information to explore the spatial distribution of different cell types in the heart tissue. However, there are three full slices of the heart in the dataset, so we will focus on the middle slice for more detailed visualisation.

# Subset to one slice for better visualisation
# Check the subset with a zoomed-in DimPlot first
DimPlot(merfish, reduction = "spatial", raster = FALSE) +
  coord_cartesian(xlim = c(-5000, 500))

# Subset the object by spatial reduction coordinates
# Grab the coordinates for the spatial reduction and convert to a data frame
spatial_coordinates <- merfish[["spatial"]] |>
  Embeddings() |>
  as.data.frame()

# Filter to keep only the region of interest
middle_heart_coordinates <- spatial_coordinates |>
  filter(SPATIAL_1 > -5000 & SPATIAL_1 < 500)

# Subset the Seurat object
cells_to_keep <- rownames(middle_heart_coordinates)
merfish_heart <- subset(merfish, cells = cells_to_keep)

We can visualise the expression of known marker genes to identify different cell types and tissue architecture. We will look at some markers for heart tissue, including cardiac fibroblasts, atrial cardiomyocytes, and ventricular cardiomyocytes.

# Visualise MYH6 and MYH7 in the selected slice
FeaturePlot(
  merfish_heart,
  features = c("MYH6", "MYH7"),
  reduction = "spatial",
  raster = FALSE
)

# Visualise fibroblast markers in the selected slice
FeaturePlot(
  merfish_heart,
  features = c("PDGFRA", "PDGFRB", "DCN", "LUM"),
  reduction = "spatial",
  raster = FALSE
)

To identify other markers for the different clusters, we can use the FindAllMarkers function on the existing cluster identities. This dataset already has metadata assigned, including Leiden clusters, communities and populations. First we will visualise each of these annotations.

# Visualise the existing annotations on the spatial slice
DimPlot(
  merfish_heart,
  reduction = "spatial",
  raster = FALSE,
  group.by = "leiden"
) +
  ggtitle("Leiden Clusters")
DimPlot(
  merfish_heart,
  reduction = "spatial",
  raster = FALSE,
  group.by = "populations"
) +
  ggtitle("Populations")
DimPlot(
  merfish_heart,
  reduction = "spatial",
  raster = FALSE,
  group.by = "communities",
  label = TRUE
) +
  ggtitle("Communities")

Of course we can decide to find all markers for any of these annotations, but for this example we will use the “communities” annotation and use the FindMarkers function to identify differentially expressed genes between two communities of interest. Particularly interesting cases here are comparing the left and right atrium and the left and right ventricles.

# Find markers for left atrium versus right atrium
markers_atria <- FindMarkers(
  merfish_heart,
  ident.1 = "Left Atria",
  ident.2 = "Right Atria",
  group.by = "communities",
  only.pos = TRUE,
  min.pct = 0.25,
  logfc.threshold = 0.25,
  pvalue.cutoff = 0.05
)

# View the top five markers for the left atrium
head(markers_atria[order(markers_atria$p_val_adj, decreasing = FALSE), ], 5)

Let’s look at one of the top genes more highly expressed in the left atrium, COL2A1, which is a collagen gene. It’s IIA isoform is transiently expressed in the developing heart.

# Visualise COL2A1 in the left atrium
FeaturePlot(
  merfish_heart,
  features = "COL2A1",
  reduction = "spatial",
  raster = FALSE
)

We can do the same for the ventricles. They are divided into communities for the outer and innerventricular regions, so we will compare these two separately for left and right ventricles and also compare outter to inner regions.

# Find markers for outer left ventricle versus outer right ventricle
markers_outer_ventricles <- FindMarkers(
  merfish_heart,
  ident.1 = "Outer-LV",
  ident.2 = "Outer-RV",
  group.by = "communities",
  only.pos = TRUE,
  min.pct = 0.25,
  logfc.threshold = 0.25,
  pvalue.cutoff = 0.05
)

# Find markers for inner left ventricle versus inner right ventricle
markers_inner_ventricles <- FindMarkers(
  merfish_heart,
  ident.1 = "Inner-LV",
  ident.2 = "Inner-RV",
  group.by = "communities",
  only.pos = TRUE,
  min.pct = 0.25,
  logfc.threshold = 0.25,
  pvalue.cutoff = 0.05
)

# Find markers for outer versus inner left ventricle
markers_left_ventricle <- FindMarkers(
  merfish_heart,
  ident.1 = "Outer-LV",
  ident.2 = "Inner-LV",
  group.by = "communities",
  only.pos = TRUE,
  min.pct = 0.25,
  logfc.threshold = 0.25,
  pvalue.cutoff = 0.05
)

# Find markers for outer versus inner right ventricle
markers_right_ventricle <- FindMarkers(
  merfish_heart,
  ident.1 = "Outer-RV",
  ident.2 = "Inner-RV",
  group.by = "communities",
  only.pos = TRUE,
  min.pct = 0.25,
  logfc.threshold = 0.25,
  pvalue.cutoff = 0.05
)

From the top results we can try looking at the FeaturePlots to identify markers for left vs right ventricles and inner vs outer regions.

Answer

# Visualise markers that distinguish the ventricles
FeaturePlot(
  merfish_heart,
  features = "IRX4",
  reduction = "spatial",
  raster = FALSE
)

# Visualise a marker for outer versus inner regions
FeaturePlot(
  merfish_heart,
  features = "HEY2",
  reduction = "spatial",
  raster = FALSE
)

As this is not a proper spatial dataset (for Seurat), we cannot use spatially aware functions like FindSpatiallyVariableFeatures or packages like BANKSY here. However, we can use the spatial coordinates to visualise the expression of genes of interest in the tissue context, and the full dataset of heart samples from different developmental stages can be used to track gene expression changes during development.

16.1 Conclusion

MERFISH is a powerful technique for spatially resolved transcriptomics, allowing for the visualisation and analysis of gene expression patterns in tissue sections. In this example, we have explored the spatial distribution of different cell types in the human heart tissue using MERFISH data. We have also identified marker genes for different communities in the heart tissue, including atrial and ventricular regions. This analysis provides insights into the cellular composition and organisation of the heart tissue, which can be further explored to understand the underlying biological processes and disease mechanisms.

16.2 Summary

Key Points

MERFISH data can be analysed using Seurat functions for single-cell data, with spatial coordinates used for visualisation.
Marker genes can be identified for different communities in the tissue, providing insights into cellular composition and organisation.
Visualisation of gene expression patterns in the tissue context can reveal spatial distribution of cell types and regions.
MERFISH can be applied in a larger scale to compare gene expression patterns across different developmental stages or disease conditions, as can be more affordable compared to other spatial transcriptomics technologies.