Spatial Transcriptomics
CSAMA 2024 Version
11 June 2024
Contents
1 Introduction
In this tutorial we walk through a typical spatial transcriptomics analysis using Bioconductor packages.
Spatial transcriptomics is a fast evolving set of technologies and we cannot cover all protocols here. However, we will show an example of spot-based protocols, namely 10X Genomics Visium, and an example of imaging-based methods, namely Nanostring CosMX. Note that this is an area of current development and some methods will likely be improved in the next months.
Most of the steps covered here, especially those for spot-based data, are described in the Best Practices ST book.
The Voyager Bioconductor package (Moses et al. 2023) also has an extensive set of tutorials that cover most of the currently available spatial transcriptomics technologies.
While not covered in this tutorial, there are packages and software tools for the analysis of spatial transcriptomics data outside of Bioconductor too. Popular tools include the Seurat R package, the Giotto R package and the SpatialData python package.
2 Visium Data
We will initially focus on 10X Genomics Visium. We start from the output of the Space Ranger preprocessing software. This is the 10X Genomics software suite that allows to pre-process the FASTQ files generated by the sequencing platform and perform alignment and quantification. We will perform exploratory data analysis (EDA) and quality control (QC). We will then cover normalization, the identification of spatially variable genes dimensionality reduction and cell type identification and .
2.1 Experimental data
We will use one sample of human brain from the dorsolateral prefrontal cortex (DLPFC) region, measured using the 10x Genomics Visium platform.
In the full dataset, there are 12 samples in total, from 3 individuals, with 2 pairs of spatially adjacent replicates (serial sections) per individual (4 samples per individual). The individuals and spatially adjacent replicates can be used as blocking factors. Each sample spans the six layers of the cortex plus white matter in a perpendicular tissue section. For the examples in this workflow we use a single sample from this dataset (sample 151673).
For more details on the dataset, see Maynard et al. (2021). The full dataset is publicly available through the spatialLIBD Bioconductor package.
2.2 The SpatialExperiment class
SpatialExperiment is a S4 class that extends SingleCellExperiment and can be used for efficiently storing and working with spatial data in R/Bioconductor.
This class is itself extended by MoleculeExperiment and SpatialFeatureExperiment, which allow to more easily work with imaging-based data. For the moment, we will use SpatialExperiment.
A more thorough overview of SingleCellExperiment can be found in Righelli et al. (2022).
We start by loading the data.
library(SpatialExperiment)
library(STexampleData)
spe <- Visium_humanDLPFC()
spe
## class: SpatialExperiment
## dim: 33538 4992
## metadata(0):
## assays(1): counts
## rownames(33538): ENSG00000243485 ENSG00000237613 ... ENSG00000277475
## ENSG00000268674
## rowData names(3): gene_id gene_name feature_type
## colnames(4992): AAACAACGAATAGTTC-1 AAACAAGTATCTCCCA-1 ...
## TTGTTTGTATTACACG-1 TTGTTTGTGTAAATTC-1
## colData names(8): barcode_id sample_id ... reference cell_count
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## spatialCoords names(2) : pxl_col_in_fullres pxl_row_in_fullres
## imgData names(4): sample_id image_id data scaleFactorThe SpatialExperiment class should be fairly familiar, since it is heavily based on the SingleCellExperiment class. In addition to the slots that you already know, a SpatialExperiment object includes the spatialCoords and the imgData slots.
head(spatialCoords(spe))
## pxl_col_in_fullres pxl_row_in_fullres
## AAACAACGAATAGTTC-1 3913 2435
## AAACAAGTATCTCCCA-1 9791 8468
## AAACAATCTACTAGCA-1 5769 2807
## AAACACCAATAACTGC-1 4068 9505
## AAACAGAGCGACTCCT-1 9271 4151
## AAACAGCTTTCAGAAG-1 3393 7583imgData(spe)
## DataFrame with 2 rows and 4 columns
## sample_id image_id data scaleFactor
## <character> <character> <list> <numeric>
## 1 sample_151673 lowres #### 0.0450045
## 2 sample_151673 hires #### 0.1500150Note that the colData DataFrame is still where most of the useful information about the spots are available.
colData(spe)
## DataFrame with 4992 rows and 8 columns
## barcode_id sample_id in_tissue array_row
## <character> <character> <integer> <integer>
## AAACAACGAATAGTTC-1 AAACAACGAATAGTTC-1 sample_151673 0 0
## AAACAAGTATCTCCCA-1 AAACAAGTATCTCCCA-1 sample_151673 1 50
## AAACAATCTACTAGCA-1 AAACAATCTACTAGCA-1 sample_151673 1 3
## AAACACCAATAACTGC-1 AAACACCAATAACTGC-1 sample_151673 1 59
## AAACAGAGCGACTCCT-1 AAACAGAGCGACTCCT-1 sample_151673 1 14
## ... ... ... ... ...
## TTGTTTCACATCCAGG-1 TTGTTTCACATCCAGG-1 sample_151673 1 58
## TTGTTTCATTAGTCTA-1 TTGTTTCATTAGTCTA-1 sample_151673 1 60
## TTGTTTCCATACAACT-1 TTGTTTCCATACAACT-1 sample_151673 1 45
## TTGTTTGTATTACACG-1 TTGTTTGTATTACACG-1 sample_151673 1 73
## TTGTTTGTGTAAATTC-1 TTGTTTGTGTAAATTC-1 sample_151673 1 7
## array_col ground_truth reference cell_count
## <integer> <character> <character> <integer>
## AAACAACGAATAGTTC-1 16 NA NA NA
## AAACAAGTATCTCCCA-1 102 Layer3 Layer3 6
## AAACAATCTACTAGCA-1 43 Layer1 Layer1 16
## AAACACCAATAACTGC-1 19 WM WM 5
## AAACAGAGCGACTCCT-1 94 Layer3 Layer3 2
## ... ... ... ... ...
## TTGTTTCACATCCAGG-1 42 WM WM 3
## TTGTTTCATTAGTCTA-1 30 WM WM 4
## TTGTTTCCATACAACT-1 27 Layer6 Layer6 3
## TTGTTTGTATTACACG-1 41 WM WM 16
## TTGTTTGTGTAAATTC-1 51 Layer2 Layer2 5In this case, we have a ground truth annotation of the spots, as well as an indication of the number of cells covered by each spots.
2.3 Visualization
One nice feature of spatial transcriptomics is that we can “see” the tissue under study and visualize the spots on top of the histology image. We will use the `ggspavis package to do so.
library(ggspavis)
plotVisium(spe)
We can also visualize the ground truth annotation of the spots.
plotSpots(spe, annotate = "ground_truth",
pal = "libd_layer_colors")
We can see that, while the Visium platform consists of a grid of spots, only some of them will be covered by the tissue: these are called “in_tissue” in SpatialExperiment.
spe$in_tissue <- as.factor(spe$in_tissue)
plotSpots(spe, in_tissue = NULL, annotate = "in_tissue")
For the rest of the workflow, we will only keep the in-tissue spots.
spe <- spe[, colData(spe)$in_tissue == 1]2.4 Quality control and normalization
Spot-based spatial transcriptomics is based on short-read sequencing of RNA molecules and as such many methods originally developed for single-cell RNA-seq can be used to analyze Visium data.
In particular, scater’s perCellQCMetrics function can be used to compute a set of metrics useful to evaluate the quality of the samples. The isOutlier function uses a data driven threshold to define cells of lower quality compared to the rest of the dataset.
Similarly to what we have done in the single-cell workflow, it is useful to identify spots with high percentages of mitochondrial reads, which may be a symptom of poor sample quality.
library(EnsDb.Hsapiens.v86)
rowData(spe)$location <- mapIds(EnsDb.Hsapiens.v86,
keys=rowData(spe)$gene_id,
column="SEQNAME", keytype="GENEID")
table(rowData(spe)$location=="MT")
##
## FALSE TRUE
## 33213 13We can now compute the QC metrics with scater.
library(scater)
spe <- addPerCellQC(spe, subsets = list(mito = which(rowData(spe)$location=="MT")))
colData(spe)
## DataFrame with 3639 rows and 14 columns
## barcode_id sample_id in_tissue array_row
## <character> <character> <factor> <integer>
## AAACAAGTATCTCCCA-1 AAACAAGTATCTCCCA-1 sample_151673 1 50
## AAACAATCTACTAGCA-1 AAACAATCTACTAGCA-1 sample_151673 1 3
## AAACACCAATAACTGC-1 AAACACCAATAACTGC-1 sample_151673 1 59
## AAACAGAGCGACTCCT-1 AAACAGAGCGACTCCT-1 sample_151673 1 14
## AAACAGCTTTCAGAAG-1 AAACAGCTTTCAGAAG-1 sample_151673 1 43
## ... ... ... ... ...
## TTGTTTCACATCCAGG-1 TTGTTTCACATCCAGG-1 sample_151673 1 58
## TTGTTTCATTAGTCTA-1 TTGTTTCATTAGTCTA-1 sample_151673 1 60
## TTGTTTCCATACAACT-1 TTGTTTCCATACAACT-1 sample_151673 1 45
## TTGTTTGTATTACACG-1 TTGTTTGTATTACACG-1 sample_151673 1 73
## TTGTTTGTGTAAATTC-1 TTGTTTGTGTAAATTC-1 sample_151673 1 7
## array_col ground_truth reference cell_count sum
## <integer> <character> <character> <integer> <numeric>
## AAACAAGTATCTCCCA-1 102 Layer3 Layer3 6 8458
## AAACAATCTACTAGCA-1 43 Layer1 Layer1 16 1667
## AAACACCAATAACTGC-1 19 WM WM 5 3769
## AAACAGAGCGACTCCT-1 94 Layer3 Layer3 2 5433
## AAACAGCTTTCAGAAG-1 9 Layer5 Layer5 4 4278
## ... ... ... ... ... ...
## TTGTTTCACATCCAGG-1 42 WM WM 3 4324
## TTGTTTCATTAGTCTA-1 30 WM WM 4 2761
## TTGTTTCCATACAACT-1 27 Layer6 Layer6 3 2322
## TTGTTTGTATTACACG-1 41 WM WM 16 2331
## TTGTTTGTGTAAATTC-1 51 Layer2 Layer2 5 6281
## detected subsets_mito_sum subsets_mito_detected
## <numeric> <numeric> <numeric>
## AAACAAGTATCTCCCA-1 3586 1407 13
## AAACAATCTACTAGCA-1 1150 204 11
## AAACACCAATAACTGC-1 1960 430 13
## AAACAGAGCGACTCCT-1 2424 1316 13
## AAACAGCTTTCAGAAG-1 2264 651 12
## ... ... ... ...
## TTGTTTCACATCCAGG-1 2170 370 12
## TTGTTTCATTAGTCTA-1 1560 314 12
## TTGTTTCCATACAACT-1 1343 476 13
## TTGTTTGTATTACACG-1 1420 308 12
## TTGTTTGTGTAAATTC-1 2927 991 13
## subsets_mito_percent total
## <numeric> <numeric>
## AAACAAGTATCTCCCA-1 16.6351 8458
## AAACAATCTACTAGCA-1 12.2376 1667
## AAACACCAATAACTGC-1 11.4089 3769
## AAACAGAGCGACTCCT-1 24.2223 5433
## AAACAGCTTTCAGAAG-1 15.2174 4278
## ... ... ...
## TTGTTTCACATCCAGG-1 8.55689 4324
## TTGTTTCATTAGTCTA-1 11.37269 2761
## TTGTTTCCATACAACT-1 20.49957 2322
## TTGTTTGTATTACACG-1 13.21321 2331
## TTGTTTGTGTAAATTC-1 15.77774 6281We can visualize the metrics and compute data-driven thresholds to flag low-quality cells.
hist(spe$subsets_mito_percent)
high.mito <- isOutlier(spe$subsets_mito_percent, type="higher")
table(high.mito)
## high.mito
## FALSE TRUE
## 3627 12hist(spe$sum)
low.umi <- isOutlier(spe$sum, type="lower")
table(low.umi)
## low.umi
## FALSE
## 3639In this case, only a few spots are discarded for high mitochondrial content. We can check whether there are spatial patterns of low quality.
colData(spe)$high.mito <- high.mito
plotSpotQC(spe, plot_type = "spot",
annotate = "high.mito")
We can plot the library sizes against the number of cells per spot (which is available for this dataset). We expect that spots with more cells will tend to have a higher number of UMIs.
plotSpotQC(spe, plot_type = "scatter",
x_metric = "cell_count", y_metric = "sum")
We see that the spots with very high cell counts also have low numbers of expressed genes. This indicates that the experiments have failed for these spots, and they should be removed. We select a threshold of 10 cells per spot. The number of spots above this threshold is relatively small, and there is a clear downward trend in the number of expressed genes above this threshold.
high.cellcount <- colData(spe)$cell_count > 10
table(high.cellcount)
## high.cellcount
## FALSE TRUE
## 3549 90spe$high.cellcount <- high.cellcount
plotSpotQC(spe, plot_type = "spot",
annotate = "high.cellcount")
The discarded spots are all on the edges of the tissue. It seems plausible that something has gone wrong with the cell segmentation on the edges of the images, so it makes sense to remove these spots.
Since only a small number of spots are flagged as low quality, we can leave them in the object and simply flag them as low-quality.
table(high.cellcount | high.mito)
##
## FALSE TRUE
## 3537 102low_quality <- which(high.cellcount | high.mito)As for quality control, normalization can be borrowed from the single-cell literature. Published work use simple methods, such as the library size normalization. We use here this simple strategy, but more advanced methods can be used (see the single-cell workflow).
spe <- logNormCounts(spe)
spe
## class: SpatialExperiment
## dim: 33538 3639
## metadata(0):
## assays(2): counts logcounts
## rownames(33538): ENSG00000243485 ENSG00000237613 ... ENSG00000277475
## ENSG00000268674
## rowData names(4): gene_id gene_name feature_type location
## colnames(3639): AAACAAGTATCTCCCA-1 AAACAATCTACTAGCA-1 ...
## TTGTTTGTATTACACG-1 TTGTTTGTGTAAATTC-1
## colData names(17): barcode_id sample_id ... high.cellcount sizeFactor
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## spatialCoords names(2) : pxl_col_in_fullres pxl_row_in_fullres
## imgData names(4): sample_id image_id data scaleFactorExercise: Use the scran package to normalize the data with the pool deconvolution method and compare the size factors with the library sizes.
3 Spatially variable genes
In a typical single-cell analysis, we focus on the highly variable genes in the hope that they will carry most of the biological information. Spatial transcriptomics allows us to identify spatially variable genes (SVGs), i.e., genes with spatially correlated patterns of expression across the tissue area.
Several methods to identify SVGs in ST data have recently been developed; here we focus on the nnSVG method (Weber et al. 2023).
In this example, we run nnSVG using a small subset of the dataset for faster runtime. We select a subset by subsampling on the set of spots and including stringent filtering for low-expressed genes. A full analysis using all spots for this dataset and default filtering parameters for Visium data from human brain tissue takes around 45 minutes for one Visium sample on a standard laptop.
library(nnSVG)
n <- 100
set.seed(123)
ix <- sample(seq_len(n), n)
spe_nnSVG <- spe[, ix]
spe_nnSVG <- filter_genes(
spe_nnSVG, filter_genes_ncounts = 10, filter_genes_pcspots = 3
)
spe_nnSVG <- logNormCounts(spe_nnSVG)
spe_nnSVG <- nnSVG(spe_nnSVG)
head(rowData(spe_nnSVG), 3)
## DataFrame with 3 rows and 18 columns
## gene_id gene_name feature_type location
## <character> <character> <character> <character>
## ENSG00000074800 ENSG00000074800 ENO1 Gene Expression 1
## ENSG00000171603 ENSG00000171603 CLSTN1 Gene Expression 1
## ENSG00000162545 ENSG00000162545 CAMK2N1 Gene Expression 1
## sigma.sq tau.sq phi loglik runtime mean
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSG00000074800 0.00689397 0.537541 21.01207 -111.492 0.009 1.71503
## ENSG00000171603 0.57866902 0.125757 23.32804 -120.270 0.010 2.09564
## ENSG00000162545 0.18362513 0.562470 4.97238 -123.776 0.011 2.58247
## var spcov prop_sv loglik_lm LR_stat rank
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSG00000074800 0.549317 0.048413 0.0126626 -111.437 -0.109614 161
## ENSG00000171603 0.721747 0.362993 0.8214763 -125.087 9.634798 59
## ENSG00000162545 0.759200 0.165932 0.2461151 -127.617 7.682376 75
## pval padj
## <numeric> <numeric>
## ENSG00000074800 1.0000000 1.0000000
## ENSG00000171603 0.0080878 0.0234409
## ENSG00000162545 0.0214681 0.0489472We can for instance inspect the top SVGs by ordering by the rank column.
rowData(spe_nnSVG)[order(rowData(spe_nnSVG)$rank),]
## DataFrame with 171 rows and 18 columns
## gene_id gene_name feature_type location
## <character> <character> <character> <character>
## ENSG00000197971 ENSG00000197971 MBP Gene Expression 18
## ENSG00000123560 ENSG00000123560 PLP1 Gene Expression X
## ENSG00000173786 ENSG00000173786 CNP Gene Expression 17
## ENSG00000109846 ENSG00000109846 CRYAB Gene Expression 11
## ENSG00000131095 ENSG00000131095 GFAP Gene Expression 17
## ... ... ... ... ...
## ENSG00000173812 ENSG00000173812 EIF1 Gene Expression 17
## ENSG00000128989 ENSG00000128989 ARPP19 Gene Expression 15
## ENSG00000150991 ENSG00000150991 UBC Gene Expression 12
## ENSG00000131143 ENSG00000131143 COX4I1 Gene Expression 16
## ENSG00000138326 ENSG00000138326 RPS24 Gene Expression 10
## sigma.sq tau.sq phi loglik runtime mean
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSG00000197971 2.83483 0.173900 1.550925 -120.999 0.011 3.74783
## ENSG00000123560 3.80992 0.513980 0.942131 -143.802 0.010 2.78344
## ENSG00000173786 2.02307 0.390224 1.151300 -126.786 0.009 1.73816
## ENSG00000109846 2.00210 0.303676 1.793466 -127.424 0.008 1.80283
## ENSG00000131095 2.49706 0.513030 2.687935 -154.650 0.009 2.01780
## ... ... ... ... ... ... ...
## ENSG00000173812 0.01120531 0.552721 11.405284 -113.2418 0.009 2.20198
## ENSG00000128989 0.01913710 0.724128 11.921848 -127.0426 0.008 2.17656
## ENSG00000150991 0.00563829 0.332138 23.772113 -87.6236 0.007 2.69460
## ENSG00000131143 0.01643595 0.654543 15.903211 -121.9345 0.008 2.69960
## ENSG00000138326 0.00746273 0.580789 0.227044 -115.3014 0.009 2.44050
## var spcov prop_sv loglik_lm LR_stat rank
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSG00000197971 2.74088 0.449245 0.942201 -191.805 141.6122 1
## ENSG00000123560 3.16595 0.701256 0.881130 -199.014 110.4231 2
## ENSG00000173786 1.92044 0.818305 0.838303 -174.019 94.4664 3
## ENSG00000109846 1.90257 0.784854 0.868298 -173.552 92.2545 4
## ENSG00000131095 2.77227 0.783134 0.829563 -192.375 75.4484 5
## ... ... ... ... ... ... ...
## ENSG00000173812 0.568611 0.0480727 0.0198702 -113.1634 -0.156703 167
## ENSG00000128989 0.749152 0.0635576 0.0257474 -126.9506 -0.183908 168
## ENSG00000150991 0.340548 0.0278663 0.0166924 -87.5314 -0.184430 169
## ENSG00000131143 0.676280 0.0474896 0.0244955 -121.8339 -0.201117 170
## ENSG00000138326 0.592135 0.0353973 0.0126863 -115.1903 -0.222210 171
## pval padj
## <numeric> <numeric>
## ENSG00000197971 0 0
## ENSG00000123560 0 0
## ENSG00000173786 0 0
## ENSG00000109846 0 0
## ENSG00000131095 0 0
## ... ... ...
## ENSG00000173812 1 1
## ENSG00000128989 1 1
## ENSG00000150991 1 1
## ENSG00000131143 1 1
## ENSG00000138326 1 1These results already offer a way to explore potentially interesting spatial signal, e.g., by plotting the expression pattern of the top ranked genes.
plotSpots(spe, annotate = "ENSG00000197971")