Analysis of spatial transcriptomics with singleCellHaystack

knitr::opts_chunk$set(
    fig.align = "center",
    message = FALSE,
    warning = TRUE
)
library(Seurat)
library(SeuratData)
library(singleCellHaystack)

set.seed(1)

We can apply singleCellHaystack to spatial transcriptomics data as well. Here we use Seurat (v3.2 or higher) and the spatial transcriptomics data available in the SeuratData package. For this example we use 10x Genomics Visium platform brain data. For more details about analyzing spatial transcriptomics with Seurat take a look at their spatial transcriptomics vignette here.

Preparing input data

library(Seurat)
library(SeuratData)
library(singleCellHaystack)

set.seed(1)

We focus on the anterior1 slice.

if (!"stxBrain" %in% SeuratData::InstalledData()[["Dataset"]]) {
  SeuratData::InstallData("stxBrain")
}

anterior1 <- LoadData("stxBrain", type = "anterior1")
anterior1

We filter genes with less 10 cells with non-zero counts. This reduces the computational time by eliminating very lowly expressed genes.

counts <- GetAssayData(anterior1, slot = "counts")
sel.ok <- Matrix::rowSums(counts > 1) > 10

anterior1 <- anterior1[sel.ok, ]
anterior1

We can plot the total number of counts per bead, superimposed on the image of the brain.

SpatialFeaturePlot(anterior1, features = "nCount_Spatial")

We normalize the data we use log normalization.

anterior1 <- NormalizeData(anterior1)

Running haystack on the spatial coordinates

At the moment there is no direct interface to apply haystack on spatial transcriptomics in Seurat. Once this version of Seurat is released we will add it. For now, we can obtain the two pieces of information required. One is the normalized counts and the other is the spatial coordinates. Then we pass to haystack_2D the coordinates of the beads and the detection matrix. Here we use a naive approach an define detected genes as those having non-zero counts.

coord <- GetTissueCoordinates(anterior1, "anterior1")
counts <- GetAssayData(anterior1, slot = "data")

res <- haystack_2D(coord$imagecol, coord$imagerow, detection = as.matrix(counts > 1))

We can check the top genes with spatial biased distribution.

sum <- show_result_haystack(res)
head(sum, n = 20)

And we can visualize the expression of the 6 top-scoring genes in the spatial plot.

SpatialFeaturePlot(anterior1, features = rownames(sum)[1:6])


Try the singleCellHaystack package in your browser

Any scripts or data that you put into this service are public.

singleCellHaystack documentation built on March 28, 2021, 9:12 a.m.