Spatial Transcriptomics Deconvolution with `SPOTlight`

Marc Elosua-Bayes

National Center for Genomic Analysis - Center for Genomic RegulationUniversity Pompeu Fabra
marc.elosua@cnag.crg.eu

Helena L. Crowell

Institute for Molecular Life Sciences, University of Zurich, SwitzerlandSIB Swiss Institute of Bioinformatics, University of Zurich, Switzerland
helena.crowell@uzh.ch

18 enero 2022

Source: vignettes/SPOTlight_kidney.Rmd

SPOTlight_kidney.Rmd

Abstract

Spatially resolved gene expression profiles are key to understand tissue organization and function. However, novel spatial transcriptomics (ST) profiling techniques lack single-cell resolution and require a combination with single-cell RNA sequencing (scRNA-seq) information to deconvolute the spatially indexed datasets. Leveraging the strengths of both data types, we developed SPOTlight, a computational tool that enables the integration of ST with scRNA-seq data to infer the location of cell types and states within a complex tissue. SPOTlight is centered around a seeded non-negative matrix factorization (NMF) regression, initialized using cell-type marker genes and non-negative least squares (NNLS) to subsequently deconvolute ST capture locations (spots).

For a more detailed explanation of SPOTlight consider looking at our manuscript: > Elosua-Bayes M, Nieto P, Mereu E, Gut I, Heyn H.
SPOTlight: seeded NMF regression to deconvolute
spatial transcriptomics spots with single-cell transcriptomes.
Nucleic Acids Res. 2021;49(9):e50. doi: 10.1093

Load packages

library(ggplot2)
library(SPOTlight)
library(SingleCellExperiment)
library(SpatialExperiment)
library(scater)
library(scran)

Introduction

What is `SPOTlight`?

SPOTlight is a tool that enables the deconvolution of cell types and cell type proportions present within each capture location comprising mixtures of cells. Originally developed for 10X’s Visium - spatial transcriptomics - technology, it can be used for all technologies returning mixtures of cells.

SPOTlight is based on learning topic profile signatures, by means of an NMFreg model, for each cell type and finding which combination of cell types fits best the spot we want to deconvolute. Find below a graphical abstract visually summarizing the key steps.

Starting point

The minimal unit of data required to run SPOTlight are:

ST (sparse) matrix with the expression, raw or normalized, where rows = genes and columns = capture locations.
Single cell (sparse) matrix with the expression, raw or normalized,
where rows = genes and columns = cells.
Vector indicating the cell identity for each column in the single cell expression matrix.

Data inputs can also be objects of class SpatialExperiment (SE), SingleCellExperiment (SCE) or Seurat objects from which the minimal data will be extracted.

Getting started

Data description

For this vignette, we will use a SE put out by 10X Genomics containing a Visium kidney slide. The raw data can be accessed here.

SCE data comes from the The Tabula Muris Consortium which contains >350,000 cells from from male and female mice belonging to six age groups, ranging from 1 to 30 months. From this dataset we will only load the kidney subset to map it to the Visium slide.

Loading the data

Both datasets are available through Biocondcutor’s ExperimentHub and can be loaded into R as follows:

library(ExperimentHub)

# initialize a Hub instance which stores a complete set of recordd
eh <- ExperimentHub()

# retrieve any records that match our keyword(s) of interest
query(eh, "Tabula Muris Senis droplet Kidney")

## ExperimentHub with 7 records
## # snapshotDate(): 2021-10-19
## # $dataprovider: The Tabula Muris Consortium
## # $species: Mus musculus
## # $rdataclass: matrix, H5File, DFrame
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
## #   rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["EH6386"]]' 
## 
##            title                                             
##   EH6386 | Tabula Muris Senis droplet Kidney colData         
##   EH6387 | Tabula Muris Senis droplet Kidney counts          
##   EH6388 | Tabula Muris Senis droplet Kidney processed counts
##   EH6389 | Tabula Muris Senis droplet Kidney rowData         
##   EH6390 | Tabula Muris Senis droplet Kidney PCA             
##   EH6391 | Tabula Muris Senis droplet Kidney tSNE            
##   EH6392 | Tabula Muris Senis droplet Kidney UMAP

query(eh, "MouseKidneyCoronal")

## ExperimentHub with 2 records
## # snapshotDate(): 2021-10-19
## # $dataprovider: 10X Genomics
## # $species: Mus musculus
## # $rdataclass: SpatialExperiment
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
## #   rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["EH6707"]]' 
## 
##            title                   
##   EH6707 | MouseKidneyCoronal      
##   EH6743 | MouseKidneyCoronal_v3.13

Load the spatial data:

library(TENxVisiumData)
spe <- MouseKidneyCoronal()
# Use symbols instead of Ensembl IDs as feature names
rownames(spe) <- rowData(spe)$symbol

Load the single cell data. Since our data comes from the Tabula Muris Sensis dataset, we can directly load the SCE object as follows:

library(TabulaMurisSenisData)
sce <- TabulaMurisSenisDroplet(tissues = "Kidney")$Kidney

Quick data exploration:

table(sce$free_annotation, sce$age)

##                                                             
##                                                                1m   3m  18m
##   CD45                                                         11   32   76
##   CD45    B cell                                                7    5   45
##   CD45    NK cell                                               1    4    8
##   CD45    T cell                                                8   18   48
##   CD45    macrophage                                           59  132  254
##   CD45    plasma cell                                           0    0    9
##   Epcam     kidney distal convoluted tubule epithelial cell    78  126  169
##   Epcam    brush cell                                          52   15   78
##   Epcam    kidney collecting duct principal cell               77  110  132
##   Epcam    kidney proximal convoluted tubule epithelial cell  945  684 1120
##   Epcam    podocyte                                            92   94   64
##   Epcam    proximal tube epithelial cell                       25  195  296
##   Epcam    thick ascending tube S epithelial cell             465  312  248
##   Pecam    Kidney cortex artery cell                           75   78   91
##   Pecam    fenestrated capillary endothelial                  164  182  164
##   Pecam    kidney capillary endothelial cell                   49   38   25
##   Stroma    fibroblast                                         15   16   37
##   Stroma    kidney mesangial cell                              80   51   18
##   nan                                                         285  238  256
##                                                             
##                                                               21m  24m  30m
##   CD45                                                         54 1010  106
##   CD45    B cell                                               54 2322   62
##   CD45    NK cell                                               2   47    4
##   CD45    T cell                                               42  846  177
##   CD45    macrophage                                          101  259  514
##   CD45    plasma cell                                          12  169   42
##   Epcam     kidney distal convoluted tubule epithelial cell   153    0  131
##   Epcam    brush cell                                         169    0   31
##   Epcam    kidney collecting duct principal cell               58    0  370
##   Epcam    kidney proximal convoluted tubule epithelial cell  868    0  817
##   Epcam    podocyte                                            66    0  170
##   Epcam    proximal tube epithelial cell                        5    0 1977
##   Epcam    thick ascending tube S epithelial cell             228    0  313
##   Pecam    Kidney cortex artery cell                           69    0  115
##   Pecam    fenestrated capillary endothelial                  134    0  211
##   Pecam    kidney capillary endothelial cell                   18    0    7
##   Stroma    fibroblast                                         13    0   80
##   Stroma    kidney mesangial cell                              22    0    7
##   nan                                                         189 1068  579

We see how there is a good representation of all the cell types across ages except at 24m. In order to reduce the potential noise introduced by age and batch effects we are going to select cells all coming from the same age.

# Keep cells from 18m mice
sce <- subset(sce, , age == "18m")
# Keep cells with clear cell type annotations
sce <- subset(sce, , !free_annotation %in% c("nan", "CD45"))

Workflow

Preprocessing

If the dataset is very large we want to downsample it to train the model, both in of number of cells and number of genes. To do this, we want to keep a representative amount of cells per cluster and the most biologically relevant genes.

In the paper we show how downsampling the number of cells per cell identity to ~100 doesn’t affect the performance of the model. Including >100 cells per cell identity provides marginal improvement while greatly increasing computational time and resources. Furthermore, restricting the gene set to the marker genes for each cell type along with up to 3.000 highly variable genes further optimizes performance and computational resources. You can find a more detailed explanation in the original paper.

Feature selection

Our first step is to get the marker genes for each cell type. We follow the Normalization procedure as described in OSCA. We first carry out library size normalization to correct for cell-specific biases:

sce <- logNormCounts(sce)

Variance modelling

We aim to identify highly variable genes that drive biological heterogeneity. By feeding these genes to the model we improve the resolution of the biological structure and reduce the technical noise.

dec <- modelGeneVar(sce)
plot(dec$mean, dec$total, xlab = "Mean log-expression", ylab = "Variance")
curve(metadata(dec)$trend(x), col = "blue", add = TRUE)

# Get the top 3000 genes.
hvg <- getTopHVGs(dec, n = 3000)

Next we obtain the marker genes for each cell identity. You can use whichever method you want as long as it returns a weight indicating the importance of that gene for that cell type. Examples include avgLogFC, AUC, pct.expressed, p-value…

colLabels(sce) <- colData(sce)$free_annotation

# Get vector indicating which genes are neither ribosomal or mitochondrial
genes <- !grepl(pattern = "^Rp[l|s]|Mt", x = rownames(sce))

# Compute marker genes
mgs <- scoreMarkers(sce, subset.row = genes)

Then we want to keep only those genes that are relevant for each cell identity:

mgs_fil <- lapply(names(mgs), function(i) {
    x <- mgs[[i]]
    # Filter and keep relevant marker genes, those with AUC > 0.8
    x <- x[x$mean.AUC > 0.8, ]
    # Sort the genes from highest to lowest weight
    x <- x[order(x$mean.AUC, decreasing = TRUE), ]
    # Add gene and cluster id to the dataframe
    x$gene <- rownames(x)
    x$cluster <- i
    data.frame(x)
})
mgs_df <- do.call(rbind, mgs_fil)

Cell Downsampling

Next, we randomly select at most 100 cells per cell identity. If a cell type is comprised of <100 cells, all the cells will be used. If we have very biologically different cell identities (B cells vs. T cells vs. Macrophages vs. Epithelial) we can use fewer cells since their transcriptional profiles will be very different. In cases when we have more transcriptionally similar cell identities we need to increase our N to capture the biological heterogeneity between them.

In our experience we have found that for this step it is better to select the cells from each cell type from the same batch if you have a joint dataset from multiple runs. This will ensure that the model removes as much signal from the batch as possible and actually learns the biological signal.

For the purpose of this vignette and to speed up the analysis, we are going to use 20 cells per cell identity:

# split cell indices by identity
idx <- split(seq(ncol(sce)), sce$free_annotation)
# downsample to at most 20 per identity & subset
n_cells <- 2
cs_keep <- lapply(idx, function(i) {
    n <- length(i)
    if (n < n_cells)
        n_cells <- n
    sample(i, n_cells)
})
sce <- sce[, unlist(cs_keep)]

Deconvolution

You are now set to run SPOTlight to deconvolute the spots!

Briefly, here is how it works:

NMF is used to factorize a matrix into two lower dimensionality matrices without negative elements. We first have an initial matrix V (SCE count matrix), which is factored into W and H. Unit variance normalization by gene is performed in V and in order to standardize discretized gene expression levels, ‘counts-umi’. Factorization is then carried out using the non-smooth NMF method, implemented in the R package NMF NMF (14). This method is intended to return sparser results during the factorization in W and H, thus promoting cell-type-specific topic profile and reducing overfitting during training. Before running factorization, we initialize each topic, column, of W with the unique marker genes for each cell type with weights. In turn, each topic of H in SPOTlight is initialized with the corresponding belongance of each cell for each topic, 1 or 0. This way, we seed the model with prior information, thus guiding it towards a biologically relevant result. This initialization also aims at reducing variability and improving the consistency between runs.
NNLS regression is used to map each capture location’s transcriptome in V’ (SE count matrix) to H’ using W as the basis. We obtain a topic profile distribution over each capture location which we can use to determine its composition.
we obtain Q, cell-type specific topic profiles, from H. We select all cells from the same cell type and compute the median of each topic for a consensus cell-type-specific topic signature. We then use NNLS to find the weights of each cell type that best fit H’ minimizing the residuals.

You can visualize the above explanation in the following workflow scheme:

res <- SPOTlight(
    x = sce,
    y = spe,
    groups = sce$free_annotation,
    mgs = mgs_df,
    hvg = hvg,
    weight_id = "mean.AUC",
    group_id = "cluster",
    gene_id = "gene")

## Scaling count matrix

## Seeding initial matrices

## Training NMF model

## Time for training: 0.2min

## Deconvoluting mixture data

Extract data from SPOTlight:

# Extract deconvolution matrix
head(mat <- res$mat)[, seq_len(3)]

##                    CD45    B cell CD45    NK cell CD45    T cell
## AAACCGTTCGTCCAGG-1     0.04085780      0.03745327     0.09125076
## AAACCTAAGCAGCCGG-1     0.00000000      0.00000000     0.09927724
## AAACGAGACGGTTGAT-1     0.05342870      0.01022877     0.11173765
## AAACGGTTGCGAACTG-1     0.02799133      0.03798122     0.07084192
## AAACTCGGTTCGCAAT-1     0.01718571      0.14946677     0.08317662
## AAACTGCTGGCTCCAA-1     0.07110740      0.05131965     0.08557376

# Extract NMF model fit
mod <- res$NMF

Visualization

In the next section we show how to visualize the data and interpret SPOTlight’s results.

Topic profiles

We first take a look at the Topic profiles. By setting facet = FALSE we want to evaluate how specific each topic signature is for each cell identity. Ideally each cell identity will have a unique topic profile associated to it as seen below.

plotTopicProfiles(
    x = mod,
    y = sce$free_annotation,
    facet = FALSE,
    min_prop = 0.01,
    ncol = 1) +
    theme(aspect.ratio = 1)

Next we also want to ensure that all the cells from the same cell identity share a similar topic profile since this will mean that SPOTlight has learned a consistent signature for all the cells from the same cell identity.

plotTopicProfiles(
    x = mod,
    y = sce$free_annotation,
    facet = TRUE,
    min_prop = 0.01,
    ncol = 6)

Lastly we can take a look at which genes the model learned for each topic. Higher values indicate that the gene is more relevant for that topic. In the below table we can see how the top genes for Topic1 are characteristic for B cells (i.e. Cd79a, Cd79b, Ms4a1…).

library(NMF)
sign <- basis(mod)
colnames(sign) <- paste0("Topic", seq_len(ncol(sign)))
head(sign)

##            Topic1        Topic2        Topic3        Topic4      Topic5
## Cd79a 0.002810294 4.234637e-319  0.000000e+00 6.705534e-221 0.001914301
## Ly6d  0.004164431 2.064875e-248 2.478523e-285  0.000000e+00 0.001618620
## Fau   0.006895937  2.490729e-03  1.553858e-03  1.942712e-03 0.001695927
## Cd37  0.002678810 3.230637e-220  1.606590e-03  1.225580e-48 0.001991599
## Cd79b 0.003905927 2.728267e-301  0.000000e+00 1.105812e-137 0.001339975
## Cd74  0.001593506 1.968397e-175 2.540083e-226  1.486929e-03 0.001602849
##              Topic6        Topic7        Topic8        Topic9       Topic10
## Cd79a  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## Ly6d  8.975658e-316  0.000000e+00 1.729897e-307 5.433815e-314  2.696091e-04
## Fau    1.692851e-60  5.041989e-16  1.945795e-66  1.568102e-11  6.454253e-04
## Cd37  1.404302e-287 4.579565e-150 2.706540e-279 8.501571e-286  0.000000e+00
## Cd79b  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## Cd74  8.529272e-230 7.334888e-239 1.694522e-234 1.525215e-204 6.501910e-292
##             Topic11       Topic12       Topic13       Topic14       Topic15
## Cd79a  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## Ly6d   0.000000e+00 1.704725e-302 1.477727e-300  0.000000e+00 4.943387e-298
## Fau    1.391518e-15  4.520660e-26  1.159283e-04  7.740333e-18  7.263718e-30
## Cd37   0.000000e+00 2.667158e-274 2.312004e-272 3.517747e-321 7.734263e-270
## Cd79b  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## Cd74  2.976077e-301 1.747066e-243 8.821218e-161  4.980218e-20 4.804072e-225
##             Topic16       Topic17
## Cd79a  0.000000e+00  0.000000e+00
## Ly6d   0.000000e+00  0.000000e+00
## Fau    2.607292e-04  3.976858e-67
## Cd37   0.000000e+00 9.653549e-320
## Cd79b  0.000000e+00  0.000000e+00
## Cd74  1.765619e-282 1.843635e-238

# This can be dynamically visualized with DT as shown below
# DT::datatable(sign, fillContainer = TRUE, filter = "top")

Spatial Correlation Matrix

See here for additional graphical parameters.

plotCorrelationMatrix(mat)

Co-localization

Now that we know which cell types are found within each spot we can make a graph representing spatial interactions where cell types will have stronger edges between them the more often we find them within the same spot.