Vesalius: MSB Manuscript Analysis

IMPORTANT NOTICE

The code provide here shows the work flow using vesalius 1.0.1.

The forzen source code can be found (here)[https://github.com/patrickCNMartin/Vesalius/releases/tag/Release_1.0.1]

Preface

The following file contains all code related to Spatial Transcriptomics analysis using Vesalius. While this file contains the code used to produce all analysis figures in the manuscript, the order of the code is somewhat different to that of the paper.

We hope that this can also serve as a vignette describing various analysis possibilities provided by Vesalius. Each section or sub-section will indicate in which figure their respective output can be found.

Please note that the code related to Vesalius benchmarking and comparison to other methods can be found here.

NOTE change path to files accordingly

Set up

First, we will set up the R environment to run Vesalius. Note that some of the packages are Vesalius dependencies but for sake of clarity, they are also listed here. We also set a seed. Some functions (such as kmeans) requires random start points. This seed ensure that you will (hopefully) recover the same territories as we had in the paper.

Libraries

library(vesalius)
library(imager)
library(imagerExtra)
library(Seurat)
library(ggplot2)
library(patchwork)
library(viridis)
library(RColorBrewer)
library(dplyr)
library(tvR)
library(sp)
library(grid)


# Standard libraries - just in case you are a command line fiend
library(utils)
library(stats)
library(graphics)
library(grDevices)
### Set seed
set.seed(1)

Vesalius Images

RGB embeddings

The first step of the Vesalius algorithm is to load and pre-process spatial transcriptomic data. For this manuscript, we used slide-seqV2 data provided on the Single Cell Portal. Our analysis mainly focused on the mouse brain and embryo (E12.5). We also used SeqScope data available here. Please note that the SeqScope dataset have already be preprocess and do not require further processing.

Once the data has been loaded and pre-processed (log normalisation, scaling, Finding variable features - handled by Seurat), the data can be parsed to Vesalius for embedding into the RGB colour space.

#------------------------------------------------------------------------------#
# Only 2 pucks were used - but ultimately you can add as many as you want
# The first one is the mouse Hippocampus, the second is mouse embryo
#------------------------------------------------------------------------------#
slideTag <- c("Puck_200115_08","Puck_190926_03")

slideBeads <-c("~/group/slide_seqV2/Puck_200115_08_bead_locations.csv",
               "~/group/slide_seqV2/Puck_190926_03_bead_locations.csv")

slideCounts <- c("~/group/slide_seqV2/Puck_200115_08.digital_expression.txt.gz",
                 "~/group/slide_seqV2/Puck_190926_03.digital_expression.txt.gz")



#------------------------------------------------------------------------------#
# Next we can set a few parameters for the pre-processing
#------------------------------------------------------------------------------#
## number of variable features
nfeatures <- 2000

## Plot output directory
plots <- "~/group/slide_seqV2/plots/"


## Output directory
IP <- "~/group/slide_seqV2/IP/"

#------------------------------------------------------------------------------#
# Looping over files and converting UMAP projection values to RGB colours
#------------------------------------------------------------------------------#
countList <- list()
for(i in seq_along(slideTag)){
  message(paste("Loading Bead file",slideTag[i]))
  #----------------------------------------------------------------------------#
  # Loading coordinates
  #----------------------------------------------------------------------------#
  bead <- ReadSlideSeq(slideBeads[i])


  message(paste("Loading Count file",slideTag[i]))
  #----------------------------------------------------------------------------#
  # Unconventional loading - however required as some data sets
  # Fail to load  - This ensures all data sets can be loaded
  # Count matrix
  #----------------------------------------------------------------------------#
  count <- read.table(slideCounts[i], header = TRUE )
  rownames(count) <- count[,1]
  count <- count[,-1]


  #----------------------------------------------------------------------------#
  # Creating seurat spatial object
  # NOTE this code is taken from the Seurat source code as it does not seem that
  # Slide seq loading function are all exported
  # If this has been updated - this section can be changed accordingly
  #----------------------------------------------------------------------------#
  count <- CreateSeuratObject(count, assay ="Spatial")
  bead <- bead[Cells(x = count)]
  DefaultAssay(object = bead) <- "Spatial"
  count[["slice1"]] <- bead

  #----------------------------------------------------------------------------#
  # Seuart pre-processing steps
  #----------------------------------------------------------------------------#
  count <- NormalizeData(count)
  count <- FindVariableFeatures(count, nfeatures = nfeatures)
  count <- ScaleData(count)

  #----------------------------------------------------------------------------#
  # Embed PCA loading into the RGB colour space.
  # NOTE that the output of this function only assigns colours to each
  # coordinates we do not have an image yet.
  #----------------------------------------------------------------------------#
  count <- rgbUMAP(count, pcs =30, conserveSparse = FALSE)
  countList[[i]] <- count
  #----------------------------------------------------------------------------#
  # exporting embeddings
  #----------------------------------------------------------------------------#

  filenames <- paste0(IP,"UMAP_to_RGB_log_",nfeatures,"_",slideTag[i],".csv")
  exportRGB.csv(count,file = filenames)


}

Image Array Creation

The next step uses the data frames generated above to build image arrays. If you have saved the data frame above, you can directly load that data frame. This function converts punctual coordinate locations into tiles. Each tile is then rasterised to include all possible pixel locations and assigned its colour code to it. We will only use one of the arrays as an example. The analysis of both Slide-seq data is shown below.

You will notice below that two of the arguments contain the key word “filter”. filterGrid removes stray beads and filterThreshold removes any tiles that might be to large (i.e tiles that fill gaps in the ST assay).

The resolution argument changes the resolution of the output image. Smaller Images decrease run time. The image resolution may be decreased but all beads are still retained.

#------------------------------------------------------------------------------#
# Loading saved data frame after PCA to RGB conversion
# Replace countList[[i]] with this data frame
#------------------------------------------------------------------------------#
#imageBrain <- utils::read.csv("path/to/file.csv",header=T, stringsAsFactors=F)

#------------------------------------------------------------------------------#
# Next we build the image
# We recommend using more than one core.
# Using first list element
#------------------------------------------------------------------------------#
imageBrain <- buildImageArray(countList[[1]],
                              filterGrid = 0.01,
                              filterThreshold = 0.9975,
                              resolution = 40,
                              cores = 10)

Image processing

The first step in this section of the analysis will be to regularise the image via Total Variance Regularisation. Vesalius also offers the possibility to equalise color histograms to increase contrast between territories. This process is not always required and is better suited for PCA embeddings. See section on PCA embeddings below.

The higher the lambda the smoother the output.

imageBrain <- regulariseImage(imageBrain,
                              lambda = 5,
                              niter = 50,
                              normalise=T)

Finally, we will iteratively segment the image. This is achieved by the iterativeSegmentation.array that offers quite a few different options. As this function uses kmeans clustering to generate color segments, the colDepth argument define the number of final colour you wish to have. This function can take multiple values for colDepth. Essentially, the image will first be smoothed and then segmented into kcolors. Then, another round of smoothing will be applied followed by another round of segmentation. This will be repeated for every value of k. In this context, we highly recommend using a colDepth with decreasing values (e.g. seq(12,6, by = -2)).

Image smoothing is handle within the segmentation function. If you wish to smooth independently, please use the smoothArrayfunction.

#------------------------------------------------------------------------------#
# colDepth = number of colour segments
# smoothIter = number of smoothing rounds PRIOR to segmentation
# methods = smoothing methods (they will be applied in the order they are given)
# sigma = sigma for gaussian blur
# box = size of median kernel
# useCenter = we use only the centre pixel
# invert = if colour should be inverted (if true then background is "white")
#------------------------------------------------------------------------------#

imageBrain <- iterativeSegmentation.array(imageBrain, colDepth =11,
                                          smoothIter = 20,
                                          method = c("iso","box"),
                                          sigma=1,box = 15,
                                          useCenter = T)

The parameters selected during this process will affect the number final territories. While case by case parameter selection would optimised the final number of territories, we have noticed that these parameters are more closely tied to ST method. We used the same parameters for all data sets produced by the same ST method.

Territory Isolation

Vesalius takes an image containing colour segments and further separates them into distinct spatial territories. For each color segment, beads are pooled into a territory as long as they are within a certain captureRadius of other beads belonging to the same colour segment. If there are still beads remaining in that color segment, the process will be repeated until all beads have been pooled into a territory. Vesalius also considers that territories that fall under a certain threshold (defined by minBar) are isolated territories. All Isolated territories will be pooled together.

#------------------------------------------------------------------------------#
# captureRadius = proportion of maximum distance between two beads.
# minBar = minimum number of cells that should be in a territory
#------------------------------------------------------------------------------#
imageBrain <- isolateTerritories.array(imageBrain,
                                        captureRadius = 0.008,
                                        minBar = 50)

Vesalius Image and Territory visualisation

Both image arrays and territories can be visualised using Vesalius plotting functions. Please note that these functions are minimalistic: we assume that users might want to use their own plotting style and this is straight forward as the output of all Vesalius functions are simple data frames. Feel free to use which ever tool you prefer to visualise your data.

This section produces:

Figure 2a

#------------------------------------------------------------------------------#
# IMPORTANT: the "split" show all territories separately - it just makes
# visualisation and territory selection easier. Please refer to this plot
# to select territories if in doubt.
# If desired, each function will return a ggplot object that can be further
# customised as demonstrated below.
#------------------------------------------------------------------------------#
imgSmoothed <- imagePlot(imageBrain, as.cimg = F,cex = 15) + theme_void()
imgSmoothed

imgTerritory <- territoryPlot(imageBrain, randomise = TRUE,cex =15 , cex.pt=0.25) +
   theme_void()+
   theme(legend.text = element_text(size = 12),
         legend.title = element_text(size=12),
         plot.title = element_text(size =15),
         legend.position = "right")

pdf("vesalius_Territory.pdf", width = 8, height =6)
imgTerritory
dev.off()

Isolating Vesalius Territories for in depth analysis

In this section, we will focus on the creation and manipulation of Vesalius images in order to produce territories used for in depth analysis. We are assuming that you have already run the embedding functions and potentially saved the embeddings intermediate file.

Slide-SeqV2 - Mouse Hippocampus

First, we load and process the mouse hippocampus data set taken from Slide-Seq. We are using Puck_200115_08 for this section of the work.

#------------------------------------------------------------------------------#
# Loading raw data into a Seurat object
# Pre-processing count data
# We keep a raw count seurat object to parse to DEG functions later on
#------------------------------------------------------------------------------#
bead <- ReadSlideSeq("~/group/slide_seqV2/Puck_200115_08_bead_locations.csv")
brainRaw<- read.table("~/group/slide_seqV2/Puck_200115_08.digital_expression.txt.gz", header = TRUE )
rownames(brainRaw) <- brainRaw[,1]
brainRaw <- brainRaw[,-1]
brainRaw <- CreateSeuratObject(brainRaw, assay ="Spatial")
bead <- bead[Cells(x = brainRaw)]
DefaultAssay(object = bead) <- "Spatial"
brainRaw[["slice1"]] <- bead

brainNorm <- NormalizeData(brainRaw)
brainNorm <- FindVariableFeatures(brainNorm, nfeatures = 2000)
brainNorm <- ScaleData(brainNorm)


#------------------------------------------------------------------------------#
# Loading and processing mouse hippocampus Slide-seqV2 data
# If you have saved the intermediate data frame
# you can load that file and parse that instead of countList[[1]]
#------------------------------------------------------------------------------#
#imageBrain <- utils::read.csv("~/group/slide_seqV2/IP/PCA_to_RGB_log_2000_slice1_Puck_200115_08.csv",
#                              header=T, stringsAsFactors=F)

imageBrain <- buildImageArray(countList[[1]],
                              filterThreshold = 0.9975,
                              resolution = 40,
                              cores = 10)




imageBrain <- regulariseImage(imageBrain, lambda = 5,
                              niter = 50, normalise=T)


imageBrain <- iterativeSegmentation.array(imageBrain, colDepth =11,
                                          smoothIter = 20,
                                          method = c("iso","box"),
                                          sigma=1,box = 15,
                                          useCenter = T)
imageBrain <- isolateTerritories.array(imageBrain,
                                        captureRadius = 0.008,
                                        minBar = 50)

Once the data has been processed, we can select territories for in depth analysis. We recommend using the territoryPlot function with the split argument set to TRUE to better visualise and select the desired territories. Here we will use both raw counts and processed counts.

This section produces:

Figure 3a
Figure 3d
DEG are included in supplementary material of the manuscript.

#------------------------------------------------------------------------------#
# seedID = numeric vector containing territories of interest
# morphologyFactor = numeric vector used to apply any morpholgical operators
# NOTE that positive value with inflate while negative values will erode.
# CACluster = CA Pyramidal layer field
# medCluster = Medial habenula and thrid ventricle

#------------------------------------------------------------------------------#

CACluster <- extractTerritories(imageBrain,brainRaw,
                                seedID = c(15),morphologyFactor=0)

medCluster <- extractTerritories(imageBrain,brainRaw,
                                 seedID = c(40),morphologyFactor=15)


#------------------------------------------------------------------------------#
# Clustering analysis carried out by Seurat
# Here we use the recommended pipeline as we are clustering barcodes and not
# producing images.
#------------------------------------------------------------------------------#
CACluster <- CACluster %>%
             SCTransform(assay = "Spatial") %>%
             RunPCA(dims =1:30) %>%
             RunUMAP(dims =1:30) %>%
             FindNeighbors() %>%
             FindClusters(resolution = 0.3)


#------------------------------------------------------------------------------#
# Finding and extracting markers for each cluster (top 15 based on log FC)
# This compares clusters between each other and does not consider all other
# barcodes. Used for annotation.
#------------------------------------------------------------------------------#
CAMarkers <- FindAllMarkers(CACluster) %>% group_by(cluster) %>%
             filter(avg_log2FC >0)

#------------------------------------------------------------------------------#
# Extracting genes that are differentially expressed in the cluster
# after comparing those barcodes to all other barcodes. This differs from
# the section above as we are not only considering barcodes in the isolated
# territory but all all other territories
# As such, we use normalised count values (brainCounts was previously preprocessed)
# Used for annotation.
#------------------------------------------------------------------------------#
CAClusterMarkers <- extractClusterMarkers(CACluster,brainNorm) %>%
                    group_by(seedTerritory) %>% top_n(20,wt = logFC)


#------------------------------------------------------------------------------#
# Custom plotting  with cell type annotation
# Barcodes were annotated manually - see supplementary table 1 in manuscript
# Colour comes from a colour blind friendly palette
#------------------------------------------------------------------------------#
CA <- FetchData(CACluster, c("UMAP_1","UMAP_2","seurat_clusters"))


CA_ISH <- c("CA1 Pyramidal Layer",
            "Astrocyte",
            "CA3 Pyramidal Layer",
            "CA1 Pyramidal Layer",
            "Microglia",
            "Neuron",
            "CA2 Pyramidal Layer",
             "Oligodendrocyte")

coordCA <- GetTissueCoordinates(CACluster)
CA <- cbind(CA,coordCA[,c("x","y")])

lvls <- CA_ISH[as.numeric(as.character(CA$seurat_clusters))+1]
CA$seurat_clusters <- lvls

CA_col <- length(unique(CA$seurat_clusters))
#cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73",
#          "#F0E442", "#0072B2", "#D55E00", "#CC79A7")#[c(2,6,8,5,4)]
ca_pal <- colorRampPalette(c("#999999", "#E69F00", "#56B4E9", "#009E73",
                                "#F0E442", "#0072B2", "#D55E00", "#CC79A7"))
cols <- ca_pal(CA_col)[sample(1:CA_col)]

coord_CA <- ggplot(CA, aes(x,y,col = as.factor(seurat_clusters))) +
            geom_point(size = 0.4, alpha = 1) +
            theme_void() +
            scale_color_manual(values = cols) +
            theme(legend.text = element_text(size = 20),
                  #axis.text = element_text(size = 20),
                  #axis.title = element_text(size = 20),
                  legend.position = "right",
                  plot.title = element_text(size = 20),
                  plot.tag = element_text(size=20),
                  plot.margin = margin(0, 0, 0.5, 0, "cm")) +
            guides(colour = guide_legend(override.aes = list(size=9)))+
            labs(colour = " ", title = " ",
                x = "X coordinates", y = "Y coordinates")


#------------------------------------------------------------------------------#
# The same process is applied to the third ventricle and medial habenula
#------------------------------------------------------------------------------#
medCluster <- medCluster %>%
             SCTransform(assay = "Spatial") %>%
             RunPCA(dims =1:30) %>%
             RunUMAP(dims =1:30) %>%
             FindNeighbors()
medCluster <- medCluster %>%FindClusters(resolution = 0.85)

medMarkers <- FindAllMarkers(medCluster) %>% group_by(cluster) %>%
             filter(avg_log2FC > 0)
medClusterMarkers <- extractClusterMarkers(medCluster,brainNorm) %>%
                     group_by(seedTerritory) %>% top_n(25,wt = logFC)

#------------------------------------------------------------------------------#
# Here we compare cluster between each other and more specifically clusters
# associated with the medial habenula and the third ventricle.
# This is where we get compartment specific gene expression
# medComp = medial habenula
# tvCompt = third ventricle
#------------------------------------------------------------------------------#
medComp <- FindMarkers(medCluster,ident.1 = 0, ident.2 = 5)
tvCompt <- FindMarkers(medCluster,ident.1 =c(2,4,8),ident.2=1)
#EpendymalComp <- FindMarkers(medCluster,ident.1 =1,ident.2=6)

#------------------------------------------------------------------------------#
# Plotting
#------------------------------------------------------------------------------#
med <- FetchData(medCluster, c("UMAP_1","UMAP_2","seurat_clusters"))
med_ISH <-c("Medial Habenula - Low",
            "Third Ventricle - Low",
             "Third Ventricle - High",
              "Endothelial Cells",
              "Third Ventricle - High",
              "Medial Habenula - High",
              "Ependymal",
              "Microglia",
              "Oligodendrocyte",
              "Oligodendrocyte")
coordmed <- GetTissueCoordinates(medCluster)
med <- cbind(med,coordmed[,c("x","y")])
lvls <- med_ISH[as.numeric(as.character(med$seurat_clusters))+1]
med$seurat_clusters <- lvls

med_col <- length(unique(med$seurat_clusters))
med_pal <- colorRampPalette(c("#999999", "#E69F00", "#56B4E9", "#009E73",
                                "#F0E442", "#0072B2", "#D55E00", "#CC79A7"))
cols <- med_pal(med_col)[sample(1:med_col)]


coord_med <- ggplot(med, aes(x,y,col = as.factor(seurat_clusters))) +
            geom_point(size = 1.5, alpha = 1) +
            theme_void() +
            scale_color_manual(values = cols) +
            theme(legend.text = element_text(size = 20),
                  legend.position = "right",
                  #axis.text = element_text(size = 20),
                  #axis.title = element_text(size = 20),,
                  plot.title = element_text(size = 20),
                  plot.tag = element_text(size=20),
                  plot.margin = margin(0, 0, 0.5, 0, "cm")) +
            guides(colour = guide_legend(override.aes = list(size=9)))+
            labs(colour = " ", title = "",
                x = "X coordinates", y = "Y coordinates")

Finally, we can visualise cell clustering in isolated territories

Figure 3a
Figure 3d

#pdf("SSV2_Brain_IsolatedTerritory.pdf", width = 10, height= 5)
coord_CA / coord_med
#dev.off()

Slide-SeqV2 - Mouse Embryo (E12.5)

We applied the same principles in the mouse embryo as in the mouse hippocampus. Here, we are using Puck_190926_03. If you have not saved the RGB embedding file please run the embedding section first and substitute object/files accordingly.

The following section produces:

Figure 2c

bead <- ReadSlideSeq("~/group/slide_seqV2/Puck_190926_03_bead_locations.csv")


embryoCounts <- read.table("~/group/slide_seqV2/Puck_190926_03.digital_expression.txt.gz", header = TRUE )
rownames(embryoCounts) <- embryoCounts[,1]
embryoCounts <- embryoCounts[,-1]

embryoCounts <- CreateSeuratObject(embryoCounts, assay ="Spatial")
bead <- bead[Cells(x = embryoCounts)]
DefaultAssay(object = bead) <- "Spatial"
embryoCounts[["slice1"]] <- bead
embryoNorm <- NormalizeData(embryoCounts)
embryoNorm <- FindVariableFeatures(embryoNorm, nfeatures = 2000)
embryoNorm <- ScaleData(embryoNorm)

#imageEmbryo <- utils::read.csv("~/group/slide_seqV2/IP/PCA_to_RGB_log_2000_slice1_Puck_190926_03.csv",
#                                 header=T, stringsAsFactors=F)

imageEmbryo <- buildImageArray(countList[[2]],
                               filterThreshold = 0.999,
                               resolution = 40,
                               cores = 5)


#imageEmbryo <- equalizeHistogram(imageEmbryo,
#                                 sleft = 2.5, sright = 2.5,invert = T)


imageEmbryo <- regulariseImage(imageEmbryo,lambda =10)


imageEmbryo <- iterativeSegmentation.array(imageEmbryo,
                                           colDepth = 12,
                                           smoothIter = 20,
                                           method = c("iso","box"),
                                           sigma = 1,
                                           box = 10,
                                           useCenter=T)


imageEmbryo <- isolateTerritories.array(imageEmbryo,
                                        captureRadius = 0.025,minBar = 40)



#----------------------------------------------------------------------------#
# Finally some plotting - Figure 2c
#----------------------------------------------------------------------------#
imgSmoothedEm <- imagePlot(imageEmbryo, as.cimg = F,cex = 15) + theme_void()
imgTerritoryEm <- territoryPlot(imageEmbryo, randomise = TRUE,cex =15 , cex.pt=0.25)+
                    theme_void()+
                    theme(legend.text = element_text(size = 12),
                          legend.title = element_text(size=12),
                          plot.title = element_text(size=15))+
                    labs(title = "Slide-seq V2 - Embryo (E12.5)")
pdf(paste0("Embryo.pdf"), width = 7, height=5.5)
imgTerritoryEm
dev.off()

In the embryo, we isolated the eye for in depth analysis. The principles applied in hippocampus still apply here. The DEG are included in the supplementary material of the manuscript.

#------------------------------------------------------------------------------#
# Extract territories
#------------------------------------------------------------------------------#

eyeCluster <- extractTerritories(imageEmbryo,embryoCounts,seedID = 17,morphologyFactor=15)

#------------------------------------------------------------------------------#
# Eye analysis
#------------------------------------------------------------------------------#
eyeCluster <- eyeCluster %>%
              SCTransform(assay = "Spatial") %>%
              RunPCA(dims =1:30) %>%
              RunUMAP(dims =1:30) %>%
              FindNeighbors()
eyeCluster <- eyeCluster %>% FindClusters(resolution = 0.6)


#------------------------------------------------------------------------------#
# Embryonic eye markers
#------------------------------------------------------------------------------#
eyeMarkers <- FindAllMarkers(eyeCluster) %>% group_by(cluster) %>%
                             filter(avg_log2FC >0)

eyeClusterMarkers <- extractClusterMarkers(eyeCluster,embryoNorm) %>%
                     group_by(seedTerritory) %>% top_n(30,wt = logFC)

#------------------------------------------------------------------------------#
# Here we compare cluster between each other and more specifically clusters
# associated with different cell type layers in the embryonic eye.
#------------------------------------------------------------------------------#
epiComp <- FindMarkers(eyeCluster,ident.1 = 1, ident.2 = 3)
lensComp <- FindMarkers(eyeCluster,ident.1 =4,ident.2=6)
#------------------------------------------------------------------------------#
# Plotting and Assigning cell types
#------------------------------------------------------------------------------#
eye <- FetchData(eyeCluster, c("UMAP_1","UMAP_2","seurat_clusters"))
eye_ISH <- c("Periocular Mesenchyme",
             "Anterior Lens Epithelial Cells",
             "Primary Lens Fiber Cells",
             "Anterior Lens Epithelial Cells - L2",
             "Lens Vesicle Cells",
             "Lens Vesicle Cells - L2",
             "Preplacodal Lens Ectoderm",
             "Primary Lens Fiber Cells",
             "Macrophage")

coordeye <- GetTissueCoordinates(eyeCluster)
eye <- cbind(eye,coordeye[,c("x","y")])

lvls <- eye_ISH[as.numeric(as.character(eye$seurat_clusters))+1]
eye$seurat_clusters <- lvls

eye_col <- length(unique(eye$seurat_clusters))

eye_pal <- colorRampPalette(c("#999999", "#E69F00", "#56B4E9", "#009E73",
                                "#F0E442", "#0072B2", "#D55E00", "#CC79A7"))
cols <- eye_pal(eye_col)[sample(1:eye_col)]




coord_eye <- ggplot(eye, aes(x,y,col = as.factor(seurat_clusters))) +
             geom_point(size = 1, alpha = 1) +
             theme_void() +
             scale_color_manual(values = cols) +
             theme(legend.text = element_text(size = 20),
                   legend.position ="right",
                   plot.title = element_text(size = 20),
                   plot.tag = element_text(size=20),
                   plot.margin = margin(0.5, 0, 0, 0, "cm")) +
             guides(colour = guide_legend(override.aes = list(size=9)))+
            labs(colour = "", title = "",
                                x = "X coordinates", y = "Y coordinates")

We can visualise the isolated eye.

Figure 3g

coord_eye

SeqScope - Liver

Vesalius is able to handle any type of spatial transcriptomic data as long as there are coordinates and gene counts provided. Here, we demonstrate Vesalius ability to recover tissue territories in another high resolution Spatial transcriptomics assay SeqScope. SeqScope data is available here. The SeqScope team provide process spatial data sets in the form of .rdsfiles containing Seurat objects as well as H&E images. An important aspect to consider when working with these datasets, is that they may contain more than one “slice”. As each data set has already been pre-processed, we do not need to pre-process the data again.

This section produces:

Figure 2e

#------------------------------------------------------------------------------#
# Load data set
#------------------------------------------------------------------------------#
seq <-readRDS("~/group/seqScope/Liver_TD_10um_annotated.rds")

#------------------------------------------------------------------------------#
# Converting to RGB
#------------------------------------------------------------------------------#
seq <- rgbUMAP(SO = seq, conserveSparse=F)

#------------------------------------------------------------------------------#
# Reducing coordinates size
# There is a lot of space between coordinates that is not required
#------------------------------------------------------------------------------#
seq$x <- seq$x/10
seq$y <- seq$y/10

#------------------------------------------------------------------------------#
# Build Image and process
#------------------------------------------------------------------------------#
image <- buildImageArray(seq,sliceID=1,resolution=30,filterThreshold=0.99, cores= 10)
tmpImage <- equalizeHistogram(image, sleft = 2.5 , sright = 2.5)
tmpImage <- regulariseImage(tmpImage, lambda = 10,invert =F)
tmpImage <- iterativeSegmentation.array(tmpImage,colDepth = 9,
  smoothIter = 10,
  method = c("box","iso"),
  acrossLevels = "mean",
  sigma = seq(1,3),
  box = 10,
  useCenter = T)
tmpImage <- isolateTerritories.array(tmpImage,minBar=5,captureRadius = 0.035)

#------------------------------------------------------------------------------#
# Get single slice for plotting example
#------------------------------------------------------------------------------#
plotTmp <- filter(tmpImage, x > 800 & x < 1700)

for(i in seq_len(length(unique(plotTmp$territory)))){
  print(i)
    plotTmp$territory[plotTmp$territory == unique(plotTmp$territory)[i]] <- i
}

pdf("seqScopeTerLiver.pdf", width = 7, height=6)
g <-territoryPlot(plotTmp,cex = 15,cex.pt=2)+
   theme_void()+
   theme(legend.text = element_text(size = 12),
         legend.title = element_text(size=12),
         plot.title = element_text(size =15),
         legend.position = "right")+
  labs(colour= "Territory nr.")
print(g)
dev.off()

SeqScope - Colon

We can apply the same process as described above but for the SeqScope Colon data set. The .rds files also contains more than one slice.

This section produces:

Figure 2f

#------------------------------------------------------------------------------#
# Load data set
#------------------------------------------------------------------------------#
seq <-readRDS("~/group/seqScope/Colon_10um_annotated.rds")

#------------------------------------------------------------------------------#
# Converting to RGB
#------------------------------------------------------------------------------#
seq <- rgbUMAP(SO = seq, conserveSparse=F)

#------------------------------------------------------------------------------#
# Reducing coordinates size
# There is a lot of space between coordinates that is not required
#------------------------------------------------------------------------------#
seq$x <- seq$x/10
seq$y <- seq$y/10

#------------------------------------------------------------------------------#
# Build Image and process
#------------------------------------------------------------------------------#
image <- buildImageArray(seq,sliceID=1,resolution=30,filterThreshold=0.99, cores= 10)
tmpImage <- equalizeHistogram(image, sleft = 2.5 , sright = 2.5)
tmpImage <- regulariseImage(tmpImage, lambda = 10,invert =F)
tmpImage <- iterativeSegmentation.array(tmpImage,colDepth = 9,
  smoothIter = 10,
  method = c("box","iso"),
  acrossLevels = "mean",
  sigma = seq(1,3),
  box = 10,
  useCenter = T)
tmpImage <- isolateTerritories.array(tmpImage,minBar=5,captureRadius = 0.035)

#------------------------------------------------------------------------------#
# Get single slice for plotting example
#------------------------------------------------------------------------------#

plotTmp <- filter(tmpImage, x > 4500 & y >900)
for(i in seq_len(length(unique(plotTmp$territory)))){
  print(i)
    plotTmp$territory[plotTmp$territory == unique(plotTmp$territory)[i]] <- i
}

pdf("seqScopeTerColon.pdf", width = 7, height=6)
g1 <-territoryPlot(plotTmp,cex = 15,cex.pt=2)+
   theme_void()+
   theme(legend.text = element_text(size = 12),
         legend.title = element_text(size=12),
         plot.title = element_text(size =15),
         legend.position = "right")+
  labs(colour= "Territory nr.")
print(g1)
dev.off()

Visualising Gene expression patterns.

Expression Pattern in the isolated CA Field

We showed that we recovered the CA2 field in isolated CA territory. This was due to the fact that we could recover subtle gene expression patterns that were otherwise lost or overshadowed by stronger patterns. In this section, we show how we can visualise these expression patterns.

This section produces:

Figure 3c
Figure S19

#------------------------------------------------------------------------------#
# First we show the overall and isolated pattern for pcp4 (Fig 3)
#------------------------------------------------------------------------------#

pcp4All <- viewGeneExpression(imageBrain,brainNorm,
                              ter = NULL, genes = "Pcp4",cex=20)+
           theme_void() +
           theme(legend.position='bottom',
                 plot.title = element_text(hjust=0.5))


pcp4 <- viewGeneExpression(imageBrain,CACluster, genes = "Pcp4",cex=20)+
        theme_void()+
        theme(legend.position='bottom',
              plot.title = element_text(hjust=0.5))
pdf("pcp4.pdf",width=5,height=6)
pcp4All
pcp4
dev.off()

#------------------------------------------------------------------------------#
# Next we do the same thing for other CA2 marker genes (Fig SX)
#------------------------------------------------------------------------------#
rgs14All <- viewGeneExpression(imageBrain,brainCounts,
                              ter = NULL, genes = "Rgs14",cex=20)+
            theme_void() +
            theme(legend.position='left',
            plot.title = element_text(hjust=0.5))
rgs14 <- viewGeneExpression(imageBrain,CACluster,
                            genes = "Rgs14",cex=20)+
         theme_void()+
         theme(legend.position='left',
               plot.title = element_text(hjust=0.5))


necab2All <- viewGeneExpression(imageBrain,brainCounts,
                                ter = NULL, genes = "Necab2",cex=20)+
             theme_void() +
             theme(legend.position='right',
                   plot.title = element_text(hjust=0.5))
necab2 <- viewGeneExpression(imageBrain,CACluster,
                             genes = "Necab2",cex=20)+
          theme_void()+
          theme(legend.position='right',
                plot.title = element_text(hjust=0.5))



pdf("CA2_inClust_overall.pdf",width = 10,height = 4)
rgs14All + necab2All + plot_layout(ncol=2,nrow = 1)
dev.off()

pdf("CA2_inClust.pdf",width = 10,height = 4)
rgs14 + necab2 + plot_layout(ncol=2,nrow =1)
dev.off()

Expression Patterns in the Isolated Eye

We can also visualise differentially expressed genes between layers in the eye.

This section produces:

Figure 3h

Cryba4 <- viewGeneExpression(imageEmbryo,eyeCluster,
                             genes = "Cryba4",cex=20)+
          theme_void()+
          theme(legend.position='none',
                plot.title = element_text(hjust=0.5))



Ccnd2 <- viewGeneExpression(imageEmbryo,eyeCluster,
                             genes = "Ccnd2",cex=20)+
          theme_void()+
          theme(legend.position='right',
                plot.title = element_text(hjust=0.5))



Pmel <- viewGeneExpression(imageEmbryo,eyeCluster,
                             genes = "Pmel",cex=20)+
          theme_void()+
          theme(legend.position='none',
                plot.title = element_text(hjust=0.5))


Aldh1a1 <- viewGeneExpression(imageEmbryo,eyeCluster,
                             genes = "Aldh1a1",cex=20)+
          theme_void()+
          theme(legend.position='right',
                plot.title = element_text(hjust=0.5))
pdf("Eye_patterns.pdf", width = 5,height =4)

Cryba4+Ccnd2+Pmel +Aldh1a1 +plot_layout(ncol=2)
dev.off()

Territory specific expression patterns within cell types

Thanks to easy territory isolation, Vesalius can easily compare the expression of cells between different anatomical structures. Here, we present how we compared the expression of different cell types between anatomical structures in the mouse hippocampus. The data is taken from Puck_200115_08. First, we need to assign cell types to each barcode. In this instance, we used RCTD to decompose cell types contained within each barcode and store annotations into a Seurat object.

You do not need to use RCTD - you only need to assign cell types to the barcodes contained in your data. Feel free to use any type of annotation tool you prefer. The code shown here applies to Seurat objects

First, we load our Seurat Object and extract homotypic barcodes. Homotypic barcodes represent barcodes only containing a single cell type according to RCTD.

ref <- readRDS("/isilonsund/NextGenSeqData/project-data/hkim/ST_project/Slide-seq/Slide-seq_hippocampus.rds")
cells <- unique(unlist(strsplit(ref@meta.data$celltype, "\\+")))
homo <- paste0(cells,"+",cells)
homo_cells <- ref@meta.data[ref@meta.data$celltype %in% homo,]
homo_cells$barcodes <- rownames(homo_cells)
homo_cells <- split(homo_cells$barcodes,homo_cells$celltype)

Next, we can select territories of interest. For this example, we compared the Cortex and the Thalamus. Note that the numerical values are the Vesalius territories you wish to select. We recommend checking the territoryPlot of your ST assay in split mode to select your territories.

We will run the comparison on all homotypic beads. Not all cell type combinations will yield differential gene expression.

This section produces:

Figure 4a
DEGs are provided in the supplementary material of the manuscript

cortex <- c(1:6,8)
thalamus <- c(28,33)
#------------------------------------------------------------------------------#
# looping over all cells types in homotypic beads
# Keeping everything in a neat little data frame
#------------------------------------------------------------------------------#
by_cell <- vector("list",length(homo_cells))
names(by_cell) <- names(homo_cells)
for(i in seq_along(homo_cells)){
    print(length(homo_cells[[i]]))
    if(length(homo_cells[[i]]) <30){
      next()
    } else {
      deg <- compareCells(imageBrain,
                          brainNorm,
                          homo_cells[[i]],
                          method ="wilcox",
                          seed=cortex,
                          logFC = 0.25,
                          pval = 0.01,
                          query=thalamus,
                          minBar = 30)
    }
    if(is.null(deg)){next()}
    deg$celltype <- rep(names(homo_cells)[i], nrow(deg))
    deg$group_1 <- rep("Cortex", nrow(deg))
    deg$group_2 <- rep("Thalamus", nrow(deg))
    by_cell[[i]] <- deg
}
by_cell <- do.call("rbind",by_cell)
write.csv(by_cell, file = "vesalius_terComp_CortexVSThalamus.csv")

Finally, we can plot the results. First, we plot all differentially expressed genes and from these we will select two examples.

#------------------------------------------------------------------------------#
# Plot everything
#------------------------------------------------------------------------------#
up <- 1
#pdf("All_genes_Vesalius_diff_territory.pdf",width =20, height=20)
for(j in seq(1,4)){
  png(paste0("diffTer/All_genes_Vesalius_diff_territory",j,".png"),width = 2000,height=2000)
  # Make the panel
  plotCols <- 4
  plotRows <- 5

  # Set up the page
  grid.newpage()
  pushViewport(viewport(layout = grid.layout(plotRows, plotCols)))
  vplayout <- function(x, y)
      viewport(layout.pos.row = x, layout.pos.col = y)

  # Make each plot, in the correct location
  for (i in 1:20) {
      curRow = ceiling(i/plotCols)
      curCol = (i-1) %% plotCols + 1
      bar <- homo_cells[[by_cell$celltype[up]]]
      gene <- by_cell$genes[up]
      p <- viewCellExpression(image = imageBrain ,
          counts = brainNorm,
          cells = bar,
          gene = gene,
          ter1 =cortex,ter2=thalamus,
          cex =17) + labs(title = paste(by_cell$celltype[up],"-",gene))
      print(p, vp = vplayout(curRow, curCol ))
      up <- up +1
  }
  dev.off()
}


#------------------------------------------------------------------------------#
# Selecting examples to plot Figure 4a
#------------------------------------------------------------------------------#
pdf("vesalius_DEG_between_Ter.pdf",height = 5,width = 12)

### Cpe in astrocytes
local <- by_cell %>% filter(genes == "Cpe" & celltype == "Astrocyte+Astrocyte")
bar <- homo_cells[["Astrocyte+Astrocyte"]]
p <- viewCellExpression(image = imageBrain ,
        counts = brainNorm,
        cells = bar,
        gene = local$genes,
        ter1 =cortex,ter2=thalamus,
        normalise =TRUE,
        cex =10) +
        labs(title = "Astrocyte - Cpe")+
        scale_shape_discrete(labels = c("Cortex", "Thalamus"))

### Nrgn in Entorhinal
local <- by_cell %>% filter(genes == "Nrgn" & celltype == "Entorhinal+Entorhinal")
bar <- homo_cells[["Entorhinal+Entorhinal"]]
p1 <- viewCellExpression(image = imageBrain ,
        counts = brainNorm,
        cells = bar,
        gene = local$genes,
        ter1 =cortex,ter2=thalamus,
        normalise =TRUE,
        cex =10) +
        labs(title = "Entorhinal - Nrgn")+
        scale_shape_discrete(labels = c("Cortex", "Thalamus"))

print(p + p1)
dev.off()

Tissue border expression

The isolation of territories in combination with the image representation of territories makes for an easy way to analyse gene expression pattern at the border between different tissues. For example, we analysed the expression patterns at the border between the Dentate Gyrus Granule Cell Layer and the Dentate Gyrus Sub Granular Zone.

This section produces:

Figure 4b
DEGs are provided in the supplementary material of the manuscript

#------------------------------------------------------------------------------#
# Get territory and run clustering
#------------------------------------------------------------------------------#
DenCluster <- extractTerritories(imageBrain, brainRaw,
                                 seedID = 26,morphologyFactor = 7)

DenCluster <- DenCluster %>% SCTransform(assay = "Spatial") %>%
              RunPCA(dims = 1:30) %>%
              RunUMAP(dims = 1:30) %>%
              FindNeighbors() %>%
              FindClusters(resolution = 0.3)

#------------------------------------------------------------------------------#
# We can extract marker gene from each cluster and compare each cluster
# to all other territories
# NOTE: for cell type annotation we used anatomical position as we just wanted
# overall territory not specific cell types
#------------------------------------------------------------------------------#
denMarkers <- FindAllMarkers(DenCluster) %>%
              group_by(cluster) %>%
              top_n(25,wt=avg_log2FC)

denClusterMarkers <- extractClusterMarkers(DenCluster,brainCounts)%>%
                     group_by(seedTerritory) %>% top_n(15,wt = logFC)




#------------------------------------------------------------------------------#
# We now have all the data we need to do some plotting
# We start off with the clustered dentate gyrus
# We only show the coordinates
#------------------------------------------------------------------------------#

outer <- FetchData(DenCluster, c("UMAP_1","UMAP_2","seurat_clusters"))

coordouter <- GetTissueCoordinates(DenCluster)
outer <- cbind(outer,coordouter[,c("x","y")])

outer_ISH <- c("DG - Granule Cell Layer",
               "DG - Granule Cell Layer",
               "DG - Sub-Granular Zone",
                "DG - Molecular Layer",
              "DG - Granule Cell Layer")

lvls <- outer_ISH[as.numeric(as.character(outer$seurat_clusters))+1]
outer$seurat_clusters <- lvls
outer_col <- length(unique(outer$seurat_clusters))
#outer_pal <- colorRampPalette(c("#999999", "#E69F00", "#56B4E9", "#009E73",
                                #"#F0E442", "#0072B2", "#D55E00", "#CC79A7"))
#cols <- outer_pal(outer_col)[sample(1:outer_col)]
cols <-c("#999999", "#E69F00", "#56B4E9", "#009E73",
         "#F0E442", "#0072B2", "#D55E00", "#CC79A7")[c(2,3,4)]



coord_outer <- ggplot(outer, aes(x,y,col = as.factor(seurat_clusters))) +
               geom_point(size = 1, alpha = 1) +
               theme_void() +
               scale_color_manual(values = cols) +
               theme(legend.text = element_text(size = 16),
                     legend.position = "right",
                     legend.title = element_text(size = 15),
                     plot.title = element_text(size = 15,hjust =0),
                     plot.tag = element_text(size=15)) +
               guides(colour = guide_legend(override.aes = list(size=9)))+
               labs(colour = " ", title = " ",
                           x = " ", y = " ")

We can visualise the clustering results

Figure 4b

#pdf("DenMap.pdf", width =9, height=9)
coord_outer
#dev.off()

Tissue Layering

Vesalius provides a way to layer a territory and investigate gene expression difference between layers. In the manuscript, we isolated the Corpus Callosum and applied morphological operators to fill in the holes and gaps. Then, we layered this isolated territory and after extracting layer specific genes, we found the Kif5a and Stmn4 genes that were highly expressed in the inner layers of the Corpus Callosum.

This sections produces:

Figure 4d
DEGs are provided in the supplementary material

#------------------------------------------------------------------------------#
# First step is to dilate the territory to fill in gaps and consider
# neighbouring territories
#------------------------------------------------------------------------------#
sc <- layerTerritory.edge(imageBrain,7,morphologyFactor=c(-35,80),layerDepth =5)


#------------------------------------------------------------------------------#
# Now we can find layer specific gene expresison
#------------------------------------------------------------------------------#
totalLayerssc <- unique(sc$layer)

Layeredsc <- lapply(seq_along(totalLayerssc), function(idx,layers,ter,counts){
                    query <- layers[!layers %in% idx]
                    res <- compareLayers(ter,counts,idx,query,logFC =0.25)
                    if(length(res) ==0) return(NULL)
                    return(res)
},totalLayerssc,sc,brainNorm)
Layeredsc <- do.call("rbind",Layeredsc)
Layeredsc <- filter(Layeredsc ,logFC >0)


#------------------------------------------------------------------------------#
# Plotting Figure 4d
#------------------------------------------------------------------------------#
pdf("Kif5a.pdf",width=9,height=7)
scPlot <- viewLayerExpression(sc,brainNorm,genes = "Kif5a",cex=20)+
          theme_void()+
          theme(plot.margin = margin(0.5,0.5,0.5,0.5,"cm"),
                plot.title = element_text(hjust = 0.5, size = 25),
                legend.text = element_text(size = 10),
                legend.title = element_text(size = 20),
                legend.position = "bottom") +
          labs(title ="Kif5a")
scPlot
dev.off()

pdf("Stmn4.pdf",width=9,height=7)
scPlot <- viewLayerExpression(sc,brainNorm,genes = "Stmn4",cex=20)+
          theme_void()+
          theme(plot.margin = margin(0.5,0.5,0.5,0.5,"cm"),
                plot.title = element_text(hjust = 0.5, size = 25),
                legend.text = element_text(size = 10),
                legend.title = element_text(size = 20),
                legend.position = "bottom") +
          labs(title ="Stmn4")
scPlot
dev.off()

UMAP vs PCA

In the section, we aim to demonstrate the differences between using UMAP or PCA embeddings in your color images.

The main aspect to consider is how much data is being included in the image. UMAP projections can contain an arbitrary number of PCs and in turn this means that more information can be included into the final images. While this does sound like a good thing, it also means that more subtle patterns might be merged together in color space.

PCA slices - on the other hand - can only show 3 PCs at a time to create RGB images. Obviously, in this case there is less information included in the final image. However, this generates images that emphasise different part of the ST assay. We can use these different slices and images to explores different sections of the ST assay and focus our attention on different territories.

This section produces:

Figure S2

#------------------------------------------------------------------------------#
# Using slide-seq data sets
#------------------------------------------------------------------------------#
slideTag <- c("Puck_200115_08","Puck_190926_03")

slideBeads <-c("~/group/slide_seqV2/Puck_200115_08_bead_locations.csv",
               "~/group/slide_seqV2/Puck_190926_03_bead_locations.csv")

slideCounts <- c("~/group/slide_seqV2/Puck_200115_08.digital_expression.txt.gz",
                 "~/group/slide_seqV2/Puck_190926_03.digital_expression.txt.gz")



#------------------------------------------------------------------------------#
# Next we can set a few parameters for the pre-processing
#------------------------------------------------------------------------------#
## number of variable features
nfeatures <- 2000

## Output directory
IP <- "~/group/slide_seqV2/IP/"

#------------------------------------------------------------------------------#
# We will convert to UMAP and PCA Images
#------------------------------------------------------------------------------#

for(i in seq_along(slideTag)){
  #----------------------------------------------------------------------------#
  # storing plots for output
  #----------------------------------------------------------------------------#
  plots <- list()

  #----------------------------------------------------------------------------#
  # Loading coordinates
  #----------------------------------------------------------------------------#
  bead <- ReadSlideSeq(slideBeads[i])

  #----------------------------------------------------------------------------#
  # Unconventional loading - however required as some data sets
  # Fail to load  - This ensures all data sets can be loaded
  # Count matrix
  #----------------------------------------------------------------------------#
  count <- read.table(slideCounts[i], header = TRUE )
  rownames(count) <- count[,1]
  count <- count[,-1]


  #----------------------------------------------------------------------------#
  # Creating seurat spatial object
  # NOTE this code is taken from the Seurat source code as it does not seem that
  # Slide seq loading function are all exported
  # If this has been updated - this section can be changed accordingly
  #----------------------------------------------------------------------------#
  count <- CreateSeuratObject(count, assay ="Spatial")
  bead <- bead[Cells(x = count)]
  DefaultAssay(object = bead) <- "Spatial"
  count[["slice1"]] <- bead

  #----------------------------------------------------------------------------#
  # Seuart pre-processing steps
  #----------------------------------------------------------------------------#
  count <- NormalizeData(count)
  count <- FindVariableFeatures(count, nfeatures = nfeatures)
  count <- ScaleData(count)

  #----------------------------------------------------------------------------#
  # Embed UMAP projections into the RGB colour space and building image
  #----------------------------------------------------------------------------#
  umaps <- rgbUMAP(count, pcs =30, conserveSparse = FALSE)
  umaps <- buildImageArray(umaps,
                          filterThreshold = 0.995,
                          resolution = 40,
                          cores=5)
  plots[[1]] <- imagePlot(umaps, as.cimg =F, cex = 12) + theme_void()
  #----------------------------------------------------------------------------#
  # Embed PCA loadings into the RGB colour space and building image
  # Here we are going through multiple slices as only 3 pcs are represented
  # at a time. We will also increase contrast to help with visualisation
  #----------------------------------------------------------------------------#
  pcs <- rgbPCA(count, slices = 5 ,conserveSparse = FALSE)
  counter <- 2
  for(j in seq_along(unique(pcs$slice))){
    tmp <- buildImageArray(pcs,
                           sliceID = j,
                           resolution=40,
                           filterThreshold=0.995,
                           cores=5)
    tmp <- equalizeHistogram(tmp,sleft=3, sright=3)
    plots[[counter]] <- imagePlot(tmp,as.cimg =F, cex=12) + theme_void()
    counter <- counter + 1
  }
  #----------------------------------------------------------------------------#
  # Plot the outputs
  #----------------------------------------------------------------------------#
  fileOut <- paste0(slideTag[i],"_UMAP_vs_PCA_Vesalius.pdf")
  pdf(fileOut, width = 15, height = 10)

  grid.newpage()
  pushViewport(viewport(layout = grid.layout(2, 3)))
  vplayout <- function(x, y)
      viewport(layout.pos.row = x, layout.pos.col = y)


  for (j in seq_along(plots)) {
      curRow <- ceiling(j/3)
      curCol <- (j-1) %% 3 + 1
      print(plots[[j]], vp = vplayout(curRow, curCol ))

  }
  dev.off()

}

As seen in these figures, images generated in the embryo show a extremely diverse set of territories between PC slices compared to the brain. The approach provided by Vesalius enables a different way of gaining insight into ST assays. The use of PCA instead of UMAP projections decomposes the variance contained in the ST assay into smaller sections that can easily be visualised. We hope that this approach can contribute to capitalising on the information contained in high resolution ST data sets.

Conclusion

We hope that this has given you insights into how Vesalius works and what you can do with this tool. All figures produced here are the ones used in the manuscript. All ISH images are taken from the Allen Brain Atlas.

Patrick C.N. Martin

Preface

Set up

Libraries

Vesalius Images

RGB embeddings

Image Array Creation

Image processing

Territory Isolation

Vesalius Image and Territory visualisation

Isolating Vesalius Territories for in depth analysis

Slide-SeqV2 - Mouse Hippocampus

Slide-SeqV2 - Mouse Embryo (E12.5)

SeqScope - Liver

SeqScope - Colon

Visualising Gene expression patterns.

Expression Pattern in the isolated CA Field

Expression Patterns in the Isolated Eye

Territory specific expression patterns within cell types

Tissue border expression

Tissue Layering

UMAP vs PCA

Conclusion