SPOTlight
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
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.
The minimal unit of data required to run SPOTlight
are:
Data inputs can also be objects of class SpatialExperiment (SE), or SingleCellExperiment (SCE).
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.
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 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:
Quick data exploration:
##
## 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.
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.
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:
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.
# Get vector indicating which genes are neither ribosomal or mitochondrial
genes <- !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)
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
# 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)
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
# We are using 5 here to speed up the process but set to 75-100 for your real
# life analysis
n_cells <- 5
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)]
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 membership 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 = 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: 0.35min
## 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 deconvolution
res <- runDeconvolution(
x = spe,
mod = mod_ls[["mod"]],
ref = mod_ls[["topic"]])
Extract data from SPOTlight
:
## CD45 B cell CD45 NK cell CD45 T cell
## AAACCGTTCGTCCAGG-1 0.03066541 0.02995862 0.04840226
## AAACCTAAGCAGCCGG-1 0.02012244 0.03629587 0.00000000
## AAACGAGACGGTTGAT-1 0.05022077 0.02101514 0.00000000
## AAACGGTTGCGAACTG-1 0.02669496 0.05262647 0.01841882
## AAACTCGGTTCGCAAT-1 0.04535987 0.17023944 0.08191323
## AAACTGCTGGCTCCAA-1 0.04455906 0.01397359 0.07082672
In the next section we show how to visualize the data and interpret
SPOTlight
’s results.
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.
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…).
## Topic1 Topic2 Topic3 Topic4 Topic5
## Cd79a 0.003478566 0.0000000000 0.000000000 0.000000e+00 2.334614e-03
## Ly6d 0.004845656 0.0000000000 0.000000000 0.000000e+00 1.059204e-03
## Fau 0.009465876 0.0022512272 0.002793259 3.524951e-03 1.242161e-03
## Cd37 0.004784734 0.0007358089 0.001869904 5.947212e-163 8.101215e-93
## Cd79b 0.005371728 0.0000000000 0.000000000 0.000000e+00 3.840430e-04
## Cd74 0.001280021 0.0000000000 0.000000000 3.441755e-03 5.045690e-04
## Topic6 Topic7 Topic8 Topic9 Topic10
## Cd79a 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## Ly6d 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## Fau 2.857013e-57 4.631754e-33 3.190149e-32 1.368013e-102 2.028879e-04
## Cd37 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## Cd79b 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.895718e-251
## Cd74 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## Topic11 Topic12 Topic13 Topic14 Topic15
## Cd79a 0.000000e+00 0.000000e+00 0.000000000 0.0000000000 0.000000e+00
## Ly6d 0.000000e+00 7.417576e-145 0.000000000 0.0000000000 0.000000e+00
## Fau 3.151606e-19 2.327201e-45 0.000627946 0.0003681967 2.271439e-04
## Cd37 0.000000e+00 0.000000e+00 0.000000000 0.0000000000 8.081291e-273
## Cd79b 0.000000e+00 1.433358e-312 0.000000000 0.0000000000 0.000000e+00
## Cd74 0.000000e+00 0.000000e+00 0.000000000 0.0000000000 0.000000e+00
## Topic16 Topic17
## Cd79a 0.000000e+00 0.000000e+00
## Ly6d 0.000000e+00 0.000000e+00
## Fau 3.772372e-05 1.326872e-53
## Cd37 0.000000e+00 0.000000e+00
## Cd79b 0.000000e+00 0.000000e+00
## Cd74 0.000000e+00 0.000000e+00
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.
See here for additional graphical parameters.
We can also visualize the cell type proportions as sections of a pie chart for each spot. You can modify the colors as you would a standard ggplot2.
ct <- colnames(mat)
mat[mat < 0.1] <- 0
# Define color palette
# (here we use 'paletteMartin' from the 'colorBlindness' package)
paletteMartin <- c(
"#000000", "#004949", "#009292", "#ff6db6", "#ffb6db",
"#490092", "#006ddb", "#b66dff", "#6db6ff", "#b6dbff",
"#920000", "#924900", "#db6d00", "#24ff24", "#ffff6d")
pal <- colorRampPalette(paletteMartin)(length(ct))
names(pal) <- ct
plotSpatialScatterpie(
x = spe,
y = mat,
cell_types = colnames(mat),
img = FALSE,
scatterpie_alpha = 1,
pie_scale = 0.4) +
scale_fill_manual(
values = pal,
breaks = names(pal))
With the image underneath - we are rotating it 90 degrees counterclockwise and mirroring across the horizontal axis to show how to align if the spots don’t overlay the image.
plotSpatialScatterpie(
x = spe,
y = mat,
cell_types = colnames(mat),
img = FALSE,
scatterpie_alpha = 1,
pie_scale = 0.4,
# Rotate the image 90 degrees counterclockwise
degrees = -90,
# Pivot the image on its x axis
axis = "h") +
scale_fill_manual(
values = pal,
breaks = names(pal))
Lastly we can also take a look at how well the model predicted the proportions for each spot. We do this by looking at the residuals of the sum of squares for each spot which indicates the amount of biological signal not explained by the model.
spe$res_ss <- res[[2]][colnames(spe)]
xy <- spatialCoords(spe)
spe$x <- xy[, 1]
spe$y <- xy[, 2]
ggcells(spe, aes(x, y, color = res_ss)) +
geom_point() +
scale_color_viridis_c() +
coord_fixed() +
theme_bw()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] NMF_0.27 cluster_2.1.6
## [3] rngtools_1.5.2 registry_0.5-1
## [5] rhdf5_2.49.0 TabulaMurisSenisData_1.11.0
## [7] TENxVisiumData_1.13.0 ExperimentHub_2.13.1
## [9] AnnotationHub_3.13.1 BiocFileCache_2.13.0
## [11] dbplyr_2.5.0 scran_1.33.1
## [13] scater_1.33.4 scuttle_1.15.4
## [15] SpatialExperiment_1.15.1 SingleCellExperiment_1.27.2
## [17] SummarizedExperiment_1.35.1 GenomicRanges_1.57.1
## [19] GenomeInfoDb_1.41.1 IRanges_2.39.2
## [21] S4Vectors_0.43.2 MatrixGenerics_1.17.0
## [23] matrixStats_1.3.0 SPOTlight_1.9.0
## [25] Biobase_2.65.0 BiocGenerics_0.51.0
## [27] ggplot2_3.5.1 BiocStyle_2.33.1
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 sys_3.4.2 jsonlite_1.8.8
## [4] magrittr_2.0.3 ggbeeswarm_0.7.2 magick_2.8.4
## [7] farver_2.1.2 rmarkdown_2.27 zlibbioc_1.51.1
## [10] vctrs_0.6.5 memoise_2.0.1 htmltools_0.5.8.1
## [13] S4Arrays_1.5.7 curl_5.2.1 BiocNeighbors_1.23.0
## [16] Rhdf5lib_1.27.0 SparseArray_1.5.31 sass_0.4.9
## [19] bslib_0.8.0 plyr_1.8.9 cachem_1.1.0
## [22] buildtools_1.0.0 igraph_2.0.3 mime_0.12
## [25] lifecycle_1.0.4 iterators_1.0.14 pkgconfig_2.0.3
## [28] rsvd_1.0.5 Matrix_1.7-0 R6_2.5.1
## [31] fastmap_1.2.0 GenomeInfoDbData_1.2.12 digest_0.6.36
## [34] colorspace_2.1-1 AnnotationDbi_1.67.0 dqrng_0.4.1
## [37] irlba_2.3.5.1 RSQLite_2.3.7 beachmat_2.21.5
## [40] labeling_0.4.3 filelock_1.0.3 gdata_3.0.0
## [43] fansi_1.0.6 nnls_1.5 polyclip_1.10-7
## [46] httr_1.4.7 abind_1.4-5 compiler_4.4.1
## [49] bit64_4.0.5 withr_3.0.1 doParallel_1.0.17
## [52] ggcorrplot_0.1.4.1 BiocParallel_1.39.0 viridis_0.6.5
## [55] DBI_1.2.3 highr_0.11 ggforce_0.4.2
## [58] HDF5Array_1.33.6 MASS_7.3-61 rappdirs_0.3.3
## [61] DelayedArray_0.31.11 rjson_0.2.21 bluster_1.15.0
## [64] gtools_3.9.5 tools_4.4.1 vipor_0.4.7
## [67] scatterpie_0.2.3 beeswarm_0.4.0 glue_1.7.0
## [70] rhdf5filters_1.17.0 grid_4.4.1 gridBase_0.4-7
## [73] reshape2_1.4.4 generics_0.1.3 gtable_0.3.5
## [76] tidyr_1.3.1 BiocSingular_1.21.2 ScaledMatrix_1.13.0
## [79] metapod_1.13.0 utf8_1.2.4 XVector_0.45.0
## [82] ggrepel_0.9.5 BiocVersion_3.20.0 foreach_1.5.2
## [85] pillar_1.9.0 stringr_1.5.1 limma_3.61.9
## [88] tweenr_2.0.3 dplyr_1.1.4 lattice_0.22-6
## [91] bit_4.0.5 tidyselect_1.2.1 locfit_1.5-9.10
## [94] maketools_1.3.0 Biostrings_2.73.1 knitr_1.48
## [97] gridExtra_2.3 edgeR_4.3.8 xfun_0.46
## [100] statmod_1.5.0 stringi_1.8.4 UCSC.utils_1.1.0
## [103] ggfun_0.1.5 yaml_2.3.10 evaluate_0.24.0
## [106] codetools_0.2-20 tibble_3.2.1 BiocManager_1.30.23
## [109] cli_3.6.3 munsell_0.5.1 jquerylib_0.1.4
## [112] Rcpp_1.0.13 png_0.1-8 parallel_4.4.1
## [115] blob_1.2.4 sparseMatrixStats_1.17.2 viridisLite_0.4.2
## [118] scales_1.3.0 purrr_1.0.2 crayon_1.5.3
## [121] rlang_1.1.4 KEGGREST_1.45.1