Digital unrolling
Last compiled: 26 September 2023
digital_unrolling.rmdIn this tutorial, we will look at how to perform digital unrolling of folded “swiss roll” tissue sections.
The data set comes from a “swiss roll” of a mouse colon which was published by M. Parigi et al. (including yours truly) The spatial transcriptomic landscape of the healing mouse intestine following damage.
Mouse colon data
Let’s load and have a quick look at the data.
se_mcolon <- readRDS(system.file("extdata/mousecolon",
"se_mcolon",
package = "semla"))
se_mcolon <- se_mcolon |>
LoadImages()## ## ── Loading H&E images ──## ## ℹ Loading image from /private/var/folders/91/twz8ld_x3f98sr9yc2hq9xpn47blx2/T/RtmpQDm3p7/temp_libpathed1b712617f8/semla/extdata/mousecolon/spatial/tissue_lowres_image.jpg## ℹ Scaled image from 541x600 to 400x444 pixels## ℹ Saving loaded H&E images as 'rasters' in Seurat object
ImagePlot(se_mcolon)
In the H&E image, we can clearly see the organization of large part of the lower intestine. The outermost part of the roll is the ascending colon and the center of the roll is the rectum.
We can identify genes with a regionalized pattern pretty easily but it would we great if we could somehow unroll the tissue so that we can investigate gene expression patterns along the proximal-distal axis.
Below is a spatial map of the expression of Fxyd4 which is up-regulated closer to the distal part of the colon.
MapFeatures(se_mcolon, features = "Fxyd4",
image_use = "raw", override_plot_dims = TRUE,
colors = c("darkblue", "cyan", "yellow", "red", "darkred"))
Pre-processing
The goal of this tutorial is to show how semla can be
used to find an “unrolled” coordinate system. To achieve this, we can
start a shiny application with CutSpatialNetwork() which
will allow us to do some interesting stuff. More about that later.
Before we get started, we need to do a little bit of pre-processing.
First, we will tile the H&E image using TileImage().
Tiling is necessary to make it possible to visualize the H&E images
in a more interactive way. The path to the raw H&E image can be
found with GetStaffli(se_mcolon)@imgs.
NB: If the image file has been moved, TileImage the path
might no longer be valid. Make sure that the path is valid before
proceeding. To get decent resolution, we you should use the
“tissue_hires_image.png” output by spaceranger.
im <- GetStaffli(se_mcolon)@imgs[1]
sprintf("File exists: %s", file.exists(im))## [1] "File exists: TRUE"Now we can load the image with magick, tile the image
and export it to a local directory. The function returns the path to the
directory the tiles were exported to. We will save this to a variable
called tilepath
library(magick)
im <- image_read(im)
tilepath <- TileImage(im = im, outpath = "~/Downloads/")The CutSpatialNetwork() function requires another file
to work properly. This file holds information about the “spatial
network” (see GetSpatialNetwork() for details) which will
be used in the shiny application to define how spots in the Visium data
are connected to each other.
To demonstrate what this means, we can generate a spatial network and visualize it.
spatnet <- GetSpatialNetwork(se_mcolon)[[1]]
library(ggplot2)
ggplot() +
geom_segment(data = spatnet, aes(x = x_start, xend = x_end, y = y_start, yend = y_end)) +
scale_y_reverse() +
coord_fixed()
In the plot above, you can see the edges that connects neighboring spots to each other.
Now let’s export this spatial network to make it available for our shiny application.
export_graph(se_mcolon, sampleID = 1, outdir = tilepath$datapath)Prepare network
Now we are ready to run the app! Once the app loads, you should see
the H&E image with spots and edges. The goal is to disconnect the
layer of the roll by cutting the edges between spots in different
layers. Cutting can be done by holding the SHIFT key while
move the cursor across edges which will make them red. You can also
“repair” edges by holding the CTRL key.
The image below illustrates how edges between two layers have been sliced.
Sliced edges are colored in red
Once all edges has been cut, you have to press the
Quit & Save button. If you just close the app without
pressing “Quit & Save”, the results will not be returned back to
R.
NB: We need to make sure to save the output to a variable. Here, we
will save the output to tidy_network.
tidy_network <- CutSpatialNetwork(se_mcolon, datadir = tilepath$datapath)For downstream steps, it is crucial that the layers have been separated appropriately. Even a single edge between two layer will make it impossible to run down stream steps reliably. The image below demonstrates a situation where the separation will fail.
The arrow points at two edges (colored in black) which still connects two separate layers
Luckily, any changes you made in the app are by default written to the spatial network file, so if you miss something, you can simply run the app again and fix the issue.
Let’s plot our full network and our filtered network:
library(tidygraph)
library(patchwork)
library(dplyr)
edges_full <- tidy_network |>
activate(edges) |>
as_tibble()
edges_filtered <- tidy_network |>
activate(edges) |>
as_tibble() |>
filter(keep)
p1 <- ggplot(edges_full,
aes(x, xend = x_end, y, yend = y_end)) +
geom_segment() +
scale_y_reverse() +
coord_fixed() +
ggtitle("Full network")
p2 <- ggplot(edges_filtered,
aes(x, xend = x_end, y, yend = y_end)) +
geom_segment() +
scale_y_reverse() +
coord_fixed() +
ggtitle("Filtered network")
p1 + p2
And just like that, we have separated the the layers!
Unrolling the network
Now comes the part where we do the actual “unrolling” using
AdjustTissueCoordinates(). There are a few important things
to note here. AdjustTissueCoordinates applies an algorithm
to sort nodes in the spatial graph between the end points. Distances
between nodes are used as a proxy for the actual distances in the
tissue. These distances between nodes are the shortest paths, i.e. the
minimum number of edges needed to visit (geodesic). For obvious reasons,
distances cannot be calculated between nodes that are not connected, and
for this reason, the function only handles connected graphs.
Unfortunately, our roll is not fully connected, so the function will only deal with the largest sub graph.
all_networks <- to_components(tidy_network |>
activate(edges) |>
filter(keep))
edges_labeled <- do.call(bind_rows, lapply(seq_along(all_networks), function(i) {
all_networks[[i]] |>
activate(edges) |>
as_tibble() |>
mutate(graph = paste0("sub graph", i))
}))
ggplot(edges_labeled, aes(x, xend = x_end, y, yend = y_end)) +
geom_segment() +
scale_y_reverse() +
coord_fixed() +
facet_grid(~graph)
The second, smaller sub graph will be discarded.
Another thing to note is that this function relies on a few assumptions about the shape of the folded tissue. If you were to try it on somethings completely different, let’s say a tissue folded in a zig-zag shape, it might not work as well.
tidy_network_adjusted <- AdjustTissueCoordinates(full_graph = tidy_network)## ℹ Removing 537 edges## ℹ Checking for disconnected graphs## ℹ More than 1 subgraph identified. Keeping the largest subgraph## ℹ Calculating pairwise geodesics between nodes in graph## ℹ Identifying end points## ℹ Finding shortest path beetween end points## ℹ Checking location of nodes relative to shortest path nodes## ✔ Finished!Our tidy_network_adjusted variable holds a
tibble with the node coordinates and the new
x_dist and y_dist coordinates. As a sanity
check, we can color our original spatial map by the new distance values
to see if they make sense.
p1 <- ggplot(tidy_network_adjusted, aes(x, y, color = x_dist)) +
geom_point() +
scale_color_gradientn(colours = RColorBrewer::brewer.pal(n = 9, name = "Spectral") |>
rev()) +
coord_fixed() +
scale_y_reverse()
p2 <- ggplot(tidy_network_adjusted, aes(x, y, color = y_dist)) +
geom_point() +
scale_color_gradientn(colours = RColorBrewer::brewer.pal(n = 9, name = "Spectral") |>
rev()) +
coord_fixed() +
scale_y_reverse()
p1 + p2
The x_dist values, which is arguably the most relevant
for regionalization patterns, looks quite good. y_dist
corresponds to the geodesic distances from the nodes along the shortest
path, where negative values indicate that the spots are located outside
the shortest path and positive values indicate that the spots are
located inside the shortest path.
Now that we have our new coordinate system, we can start looking at regionalization of gene expression along the proximal-distal axis.
library(tibble)
library(tidyr)
selected_genes <- c("Fxyd4", "Slc37a2", "Hmgcs2", "Cyp4b1")
se_mcolon <- ScaleData(se_mcolon)
gg <- FetchData(se_mcolon, vars = selected_genes, slot = "scale.data") |>
rownames_to_column(var = "name") |>
left_join(y = tidy_network_adjusted, by = "name") |>
pivot_longer(all_of(selected_genes), names_to = "variable", values_to = "value")
ggplot(gg, aes(x_dist, value, color = variable)) +
geom_smooth()
With these coordinates in our hands, it is quite straightforward to extract genes whose expression depend on distance along the proximal-distal axis. This will not be part of this tutorial, but below is a rather simple way of sorting genes by expression in a heatmap. Here, it becomes quite clear what genes are regionalized. There is one part of the heatmap that stands out, just on the right side of the center, which are the genes that are expressed in a lymphoid structure located in that part of the colon.
sorted_spots <- tidy_network_adjusted |>
arrange(x_dist) |>
pull(name)
# Get scale data
scale_data <- GetAssayData(se_mcolon,
slot = "scale.data",
assay = "Spatial")[VariableFeatures(se_mcolon), sorted_spots] |>
as.matrix()
# Bin data
bins <- cut(1:ncol(scale_data), breaks = 50)
binMat <- do.call(cbind, lapply(levels(bins), function(lvl) {
rowMeans(scale_data[, bins == lvl])
}))
pheatmap::pheatmap(binMat, cluster_cols = FALSE, border_color = FALSE,
breaks = seq(-1, 1, length.out = 50),
color = scico::scico(palette = "bam", n = 51))
Other applications
The tools described here could also be used for other applications.
Maybe we are interested in a specific part of a tissue where we want to
define an axis. Letäs make a manual selection of the outermost layer of
the mouse brain and subset the Seurat object to contain
spots from this layer.
se_mbrain <- readRDS(system.file("extdata/mousebrain", "se_mbrain", package = "semla"))
se_mbrain <- LoadImages(se_mbrain)
se_mbrain <- FeatureViewer(se_mbrain)
se_mbrain_subset <- SubsetSTData(se_mbrain, expression = region == "outer_layer")
MapLabels(se_mbrain_subset, column_name = "region", image_use = "raw")
Then we can create a spatial network from this layer and modify it a
bit to make it compatible with AdjustTissueCoordinates:
spatial_network <- GetSpatialNetwork(se_mbrain_subset)[[1]]
nodes <- spatial_network |> select(from, x_start, y_start) |> group_by(from) |>
slice_head(n = 1) |> rename(name = from, x = x_start, y = y_start)
edges <- spatial_network |> select(from, to) |> mutate(keep = TRUE)
spatial_network <- tidygraph::tbl_graph(nodes = nodes, edges = edges, directed = FALSE)Now that we have the spatial network, we can run
AdjustTissueCoordinates to compute distances between the
end points of the selected region:
# Run adjust tissue coordinates
test <- AdjustTissueCoordinates(full_graph = spatial_network)
se_mbrain_subset$x_dist <- test[match(test$name, colnames(se_mbrain_subset)), "x_dist", drop = TRUE]
# Plot x_dist
MapFeatures(se_mbrain_subset, features = "x_dist", image_use = "raw",
colors = RColorBrewer::brewer.pal(n = 9, name = "Spectral"),
override_plot_dims = TRUE, pt_size = 2.5) &
theme(plot.title = element_blank())
Package version
-
semla: 1.1.6
Session info
## R version 4.2.3 (2023-03-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices datasets utils methods base
##
## other attached packages:
## [1] tidyr_1.3.0 tibble_3.2.1 patchwork_1.1.3 tidygraph_1.2.3
## [5] semla_1.1.6 ggplot2_3.4.3 dplyr_1.1.3 SeuratObject_4.1.3
## [9] Seurat_4.3.0.1
##
## loaded via a namespace (and not attached):
## [1] systemfonts_1.0.4 plyr_1.8.8 igraph_1.5.1
## [4] lazyeval_0.2.2 sp_2.0-0 splines_4.2.3
## [7] listenv_0.9.0 scattermore_1.2 digest_0.6.33
## [10] htmltools_0.5.6 magick_2.7.4 fansi_1.0.4
## [13] magrittr_2.0.3 memoise_2.0.1 tensor_1.5
## [16] cluster_2.1.4 ROCR_1.0-11 globals_0.16.2
## [19] matrixStats_1.0.0 pkgdown_2.0.7 spatstat.sparse_3.0-2
## [22] colorspace_2.1-0 ggrepel_0.9.3 textshaping_0.3.6
## [25] xfun_0.40 jsonlite_1.8.7 progressr_0.14.0
## [28] spatstat.data_3.0-1 zeallot_0.1.0 survival_3.5-7
## [31] zoo_1.8-12 glue_1.6.2 polyclip_1.10-4
## [34] gtable_0.3.4 leiden_0.4.3 scico_1.5.0
## [37] future.apply_1.11.0 abind_1.4-5 scales_1.2.1
## [40] pheatmap_1.0.12 spatstat.random_3.1-5 miniUI_0.1.1.1
## [43] Rcpp_1.0.11 viridisLite_0.4.2 xtable_1.8-4
## [46] reticulate_1.31 htmlwidgets_1.6.2 httr_1.4.7
## [49] RColorBrewer_1.1-3 ellipsis_0.3.2 ica_1.0-3
## [52] pkgconfig_2.0.3 farver_2.1.1 sass_0.4.7
## [55] uwot_0.1.16 deldir_1.0-9 utf8_1.2.3
## [58] tidyselect_1.2.0 labeling_0.4.3 rlang_1.1.1
## [61] reshape2_1.4.4 later_1.3.1 munsell_0.5.0
## [64] tools_4.2.3 cachem_1.0.8 cli_3.6.1
## [67] dbscan_1.1-11 generics_0.1.3 ggridges_0.5.4
## [70] evaluate_0.21 stringr_1.5.0 fastmap_1.1.1
## [73] yaml_2.3.7 ragg_1.2.5 goftest_1.2-3
## [76] knitr_1.43 fs_1.6.3 fitdistrplus_1.1-11
## [79] purrr_1.0.2 RANN_2.6.1 pbapply_1.7-2
## [82] future_1.33.0 nlme_3.1-163 mime_0.12
## [85] compiler_4.2.3 plotly_4.10.2 png_0.1-8
## [88] spatstat.utils_3.0-3 bslib_0.5.1 stringi_1.7.12
## [91] highr_0.10 desc_1.4.2 forcats_1.0.0
## [94] lattice_0.21-8 Matrix_1.5-1 shinyjs_2.1.0
## [97] vctrs_0.6.3 pillar_1.9.0 lifecycle_1.0.3
## [100] BiocManager_1.30.22 spatstat.geom_3.2-4 lmtest_0.9-40
## [103] jquerylib_0.1.4 RcppAnnoy_0.0.21 data.table_1.14.8
## [106] cowplot_1.1.1 irlba_2.3.5.1 httpuv_1.6.11
## [109] R6_2.5.1 promises_1.2.1 renv_1.0.2
## [112] KernSmooth_2.23-22 gridExtra_2.3 parallelly_1.36.0
## [115] codetools_0.2-19 MASS_7.3-60 rprojroot_2.0.3
## [118] withr_2.5.0 sctransform_0.3.5 mgcv_1.9-0
## [121] parallel_4.2.3 grid_4.2.3 rmarkdown_2.24
## [124] Rtsne_0.16 spatstat.explore_3.2-1 shiny_1.7.5