SpatialDE by Svensson et al., 2018, is a method to identify spatially variable genes (SVGs) in spatially resolved transcriptomics data.
You can install spatialDE from Bioconductor with the following code:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("spatialDE")Reproducing the example analysis from the original SpatialDE Python package.
library(spatialDE)
library(ggplot2)Files originally retrieved from SpatialDE GitHub repository from the following links: https://github.com/Teichlab/SpatialDE/blob/master/Analysis/MouseOB/data/Rep11_MOB_0.csv https://github.com/Teichlab/SpatialDE/blob/master/Analysis/MouseOB/MOB_sample_info.csv
# Expression file used in python SpatialDE.
data("Rep11_MOB_0")
# Sample Info file used in python SpatialDE
data("MOB_sample_info")The Rep11_MOB_0 object contains spatial expression data for
16218 genes on 262 spots, with genes as rows
and spots as columns.
Rep11_MOB_0[1:5, 1:5]
#> 16.92x9.015 16.945x11.075 16.97x10.118 16.939x12.132 16.949x13.055
#> Nrf1 1 0 0 1 0
#> Zbtb5 1 0 1 0 0
#> Ccnl1 1 3 1 1 0
#> Lrrfip1 2 2 0 0 3
#> Bbs1 1 2 0 4 0
dim(Rep11_MOB_0)
#> [1] 16218 262The MOB_sample_info object contains a data.frame with coordinates for each
spot.
head(MOB_sample_info)Rep11_MOB_0 <- Rep11_MOB_0[rowSums(Rep11_MOB_0) >= 3, ]Rep11_MOB_0 <- Rep11_MOB_0[, row.names(MOB_sample_info)]
MOB_sample_info$total_counts <- colSums(Rep11_MOB_0)
head(MOB_sample_info)MOB_sample_infoX <- MOB_sample_info[, c("x", "y")]
head(X)stabilizeThe SpatialDE method assumes normally distributed data, so we stabilize the variance of the negative binomial distributed counts data using Anscombe’s approximation.
The stabilize() function takes as input a data.frame of expression values with samples in columns and genes in rows. Thus, in this case, we have to transpose the data.
norm_expr <- stabilize(Rep11_MOB_0)
#> + /home/biocbuild/.cache/R/basilisk/1.18.0/0/bin/conda create --yes --prefix /home/biocbuild/.cache/R/basilisk/1.18.0/spatialDE/1.12.0/env 'python=3.10.14' --quiet -c conda-forge --override-channels
#> + /home/biocbuild/.cache/R/basilisk/1.18.0/0/bin/conda install --yes --prefix /home/biocbuild/.cache/R/basilisk/1.18.0/spatialDE/1.12.0/env 'python=3.10.14' -c conda-forge --override-channels
#> + /home/biocbuild/.cache/R/basilisk/1.18.0/0/bin/conda install --yes --prefix /home/biocbuild/.cache/R/basilisk/1.18.0/spatialDE/1.12.0/env -c conda-forge 'python=3.10.14' 'numpy=1.23.5' 'scipy=1.9.3' 'patsy=0.5.3' 'pandas=1.5.2' --override-channels
norm_expr[1:5, 1:5]
#> 16.92x9.015 16.945x11.075 16.97x10.118 16.939x12.132 16.949x13.055
#> Nrf1 1.227749 0.8810934 0.8810934 1.2277491 0.8810934
#> Zbtb5 1.227749 0.8810934 1.2277491 0.8810934 0.8810934
#> Ccnl1 1.227749 1.6889027 1.2277491 1.2277491 0.8810934
#> Lrrfip1 1.484676 1.4846765 0.8810934 0.8810934 1.6889027
#> Bbs1 1.227749 1.4846765 0.8810934 1.8584110 0.8810934regress_outNext, we account for differences in library size between the samples by regressing out the effect of the total counts for each gene using linear regression.
resid_expr <- regress_out(norm_expr, sample_info = MOB_sample_info)
resid_expr[1:5, 1:5]
#> 16.92x9.015 16.945x11.075 16.97x10.118 16.939x12.132 16.949x13.055
#> Nrf1 -0.75226761 -1.2352000 -1.0164479 -0.7903289 -1.0973214
#> Zbtb5 0.09242373 -0.3323719 0.1397144 -0.2760560 -0.2533134
#> Ccnl1 -2.77597164 -2.5903783 -2.6092013 -2.8529340 -3.1193883
#> Lrrfip1 -1.92331333 -2.1578718 -2.3849405 -2.5924072 -1.7163300
#> Bbs1 -1.12186064 -1.0266476 -1.3706460 -0.5363646 -1.4666155runTo reduce running time, the SpatialDE test is run on a subset of 1000 genes. Running it on the complete data set takes about 10 minutes.
# For this example, run spatialDE on the first 1000 genes
sample_resid_expr <- head(resid_expr, 1000)
results <- spatialDE::run(sample_resid_expr, coordinates = X)
head(results[order(results$qval), ])model_searchFinally, we can classify the DE genes to interpetable DE classes using the model_search function.
We apply the model search on filtered DE results, using a threshold of 0.05 for the Q-values.
de_results <- results[results$qval < 0.05, ]
ms_results <- model_search(
sample_resid_expr,
coordinates = X,
de_results = de_results
)
# To show ms_results sorted on qvalue, uncomment the following line
# head(ms_results[order(ms_results$qval), ])
head(ms_results)spatial_patternsFurthermore, we can group spatially variable genes (SVGs) into spatial patterns using automatic expression histology (AEH).
sp <- spatial_patterns(
sample_resid_expr,
coordinates = X,
de_results = de_results,
n_patterns = 4L, length = 1.5
)
sp$pattern_resultsVisualizing one of the most significant genes.
gene <- "Pcp4"
ggplot(data = MOB_sample_info, aes(x = x, y = y, color = norm_expr[gene, ])) +
geom_point(size = 7) +
ggtitle(gene) +
scale_color_viridis_c() +
labs(color = gene)