Spatial Transcriptomics Deconvolution with SPOTlight (2023)

1.1 What is SPOTlight?

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

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

1.2 Starting point

The minimal unit of data required to run SPOTlight are:

• ST (sparse) matrix with the expression, raw or normalized, where rows = genesand 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 cellexpression matrix.

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

2.1 Data description

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

SCE data comes from the The Tabula MurisConsortium 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 kidneysubset to map it to the Visium slide.

2.2 Loading the data

Both datasets are available through Biocondcutor packages and can be loaded into R as follows.`Load the spatial data:

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

Load the single cell data. Since our data comes from theTabula Muris Sensisdataset, 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 agesexcept at 24m. In order to reduce the potential noise introduced by age andbatch effects we are going to select cells all coming from the same age.

# Keep cells from 18m micesce <- sce[, sce$age == "18m"]# Keep cells with clear cell type annotationssce <- sce[, !sce$free_annotation %in% c("nan", "CD45")]

3.1 Preprocessing

If the dataset is very large we want to downsample it to train the model, bothin of number of cells and number of genes. To do this, we want to keep arepresentative amount of cells per cluster and the most biologically relevantgenes.

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 percell identity provides marginal improvement while greatly increasingcomputational time and resources. Furthermore, restricting the gene set to themarker genes for each cell type along with up to 3.000 highly variable genesfurther optimizes performance and computational resources. You can find a moredetailed explanation in the originalpaper.

sce <- logNormCounts(sce)

# Get vector indicating which genes are neither ribosomal or mitochondrialgenes <- !grepl(pattern = "^Rp[l|s]|Mt", x = rownames(sce))dec <- modelGeneVar(sce, subset.row = genes)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)

colLabels(sce) <- colData(sce)$free_annotation# Compute marker genesmgs <- scoreMarkers(sce, subset.row = genes)

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)

# split cell indices by identityidx <- split(seq(ncol(sce)), sce$free_annotation)# downsample to at most 20 per identity & subset# We are using 5 here to speed up the process but set to 75-100 for your real# life analysisn_cells <- 5cs_keep <- lapply(idx, function(i) { n <- length(i) if (n < n_cells) n_cells <- n sample(i, n_cells)})sce <- sce[, unlist(cs_keep)]

3.2 Deconvolution

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

Briefly, here is how it works:

1. NMF is used to factorize a matrix into two lower dimensionality matriceswithout 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 performedin 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 isintended to return sparser results during the factorization in W and H, thuspromoting cell-type-specific topic profile and reducing overfitting duringtraining. Before running factorization, we initialize each topic, column,of W with the unique marker genes for each cell type with weights. In turn, eachtopic of H in SPOTlight is initialized with the corresponding belongance of eachcell 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 alsoaims at reducing variability and improving the consistency between runs.
   

   (Video) Fast and Easy Single Nucleus and Spatial Transcriptomics | PanHunter Hands-on Demo

2. 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 profiledistribution over each capture location which we can use to determine itscomposition.
   

3. we obtain Q, cell-type specific topic profiles, from H. We select all cellsfrom the same cell type and compute the median of each topic for a consensuscell-type-specific topic signature. We then use NNLS to find the weights of eachcell 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 = as.character(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: 1.04min

## Deconvoluting mixture data

Alternatively you can run SPOTlight in two steps so that you can have the trained model. Having the trained model allows you to reuse with other datasets you also want to deconvolute with the same reference. This allows you to skip the training step, the most time consuming and computationally expensive.

mod_ls <- trainNMF( x = sce, y = spe, groups = sce$type, mgs = mgs, weight_id = "weight", group_id = "type", gene_id = "gene") # Run deconvolutionres <- runDeconvolution( x = spe, mod = mod_ls[["mod"]], ref = mod_ls[["topic"]])

Extract data from SPOTlight:

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

## CD45 B cell CD45 NK cell CD45 T cell## AAACCGTTCGTCCAGG-1 0.02581825 0.02962119 0.05933737## AAACCTAAGCAGCCGG-1 0.01310012 0.03845660 0.00000000## AAACGAGACGGTTGAT-1 0.02743782 0.01962589 0.02788517## AAACGGTTGCGAACTG-1 0.01229016 0.02945113 0.04741127## AAACTCGGTTCGCAAT-1 0.00000000 0.17266965 0.09930546## AAACTGCTGGCTCCAA-1 0.04551176 0.04695416 0.02766609

# Extract NMF model fitmod <- res$NMF

4.1 Topic profiles

We first take a look at the Topic profiles. By setting facet = FALSE we want toevaluate how specific each topic signature is for each cell identity.Ideally each cell identity will have a unique topic profile associated to it asseen 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 sharea similar topic profile since this will mean that SPOTlight has learned aconsistent 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)