17 Interactivity with the R/Spatial ecosystem
Jeff Sheridan
August 7th 2024
17.1 Visium technology
Figure 17.1: Overview of Visium. Source: 10X Genomics.
Visium by 10x Genomics is a spatial gene expression platform that allows for the mapping of gene expression to high-resolution histology through RNA sequencing The process involves placing a tissue section on a specially prepared slide with an array of barcoded spots, which are 55 µm in diameter with a spot to spot distance of 100 µm. Each spot contains unique barcodes that capture the mRNA from the tissue section, preserving the spatial information. After the tissue is imaged and RNA is captured, the mRNA is sequenced, and the data is mapped back to the tissue”s spatial coordinates. This technology is particularly useful in understanding complex tissue environments, such as tumors, by providing insights into how gene expression varies across different regions.
17.2 Gene expression interpolation through kriging
Low resolution spatial data typically covers multiple cells making it difficult to delineate the cell contribution to gene expression. Using a process called kriging we can interpolate gene expression and map it to the single cell level from low resolution datasets. Kriging is a spatial interpolation technique that estimates unknown values at specific locations by weighing nearby known values based on distance and spatial trends. It uses a model to account for both the distance between points and the overall pattern in the data to make accurate predictions. By taking discrete measurement spots, such as those used for visium, we can interpolate gene expression to a finer scale using kriging.
Figure 17.2: 2D example of kriging. By Antro5 at English Wikipedia, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=52740780
By Antro5 at English Wikipedia, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=52740780
17.2.1 Dataset
For this tutorial we’ll be using the mouse brain dataset described in section 6. Visium datasets require a high resolution H&E or IF image to align spots to. Using these images we can identify individual nuclei and cells to be used for kriging. Identifying nuclei is outside the scope of the current tutorial but is required to perform kriging.
17.2.2 Generating a geojson file of nuclei location
For the following sections we will need to create a geojson that contains polygon information for the nuclei in the sample. We will be providing this in the following link, however when using for your own datasets this will need to be done outside of Giotto. A tutorial for this using qupath can be found here.
17.3 Downloading the dataset
We first need to import a dataset that we want to perform kriging on.
data_directory <- "data/03_session6"
dir.create(data_directory, showWarnings = F)
download.file(url = "https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Adult_Mouse_Brain/V1_Adult_Mouse_Brain_raw_feature_bc_matrix.tar.gz",
destfile = file.path(data_directory, "V1_Adult_Mouse_Brain_raw_feature_bc_matrix.tar.gz"))
download.file(url = "https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Adult_Mouse_Brain/V1_Adult_Mouse_Brain_spatial.tar.gz",
destfile = file.path(data_directory, "V1_Adult_Mouse_Brain_spatial.tar.gz"))17.5 Downloading giotto object and nuclei segmentation
We will need nuclei/cell segmentations to perform the kriging. Later in the tutorial we’ll also be using a pre-made giotto object. Download them using the following:
destfile <- file.path(data_directory, "subcellular_gobject.zip")
options(timeout = Inf) # Needed to download large files
download.file("https://zenodo.org/records/13144556/files/Day3_Session6.zip?download=1", destfile = destfile)
unzip(file.path(data_directory, "subcellular_gobject.zip"), exdir = data_directory)17.6 Importing visium data
We’re going to begin by creating a Giotto object for the visium mouse brain dataset. This tutorial won’t go into detail about each of these steps as these have been covered for this dataset in section 6. To get the best results when performing gene expression interpolation we need to identify spatially distinct genes. Therefore, we need to perform nearest neighbor to create a spatial network.
If you have a Giotto object from day 1 session 5, feel free to load that in and skip this first step.
library(Giotto)
save_directory <- "results/03_session6"
visium_save_directory <- file.path(save_directory, "visium_mouse_brain")
subcell_save_directory <- file.path(save_directory, "pseudo_subcellular/")
instrs <- createGiottoInstructions(show_plot = TRUE,
save_plot = TRUE,
save_dir = visium_save_directory)
v_brain <- createGiottoVisiumObject(data_directory,
gene_column_index = 2,
instructions = instrs)
# Subset to in tissue only
cm <- pDataDT(v_brain)
in_tissue_barcodes <- cm[in_tissue == 1]$cell_ID
v_brain <- subsetGiotto(v_brain,
cell_ids = in_tissue_barcodes)
# Filter
v_brain <- filterGiotto(gobject = v_brain,
expression_threshold = 1,
feat_det_in_min_cells = 50,
min_det_feats_per_cell = 1000,
expression_values = "raw")
# Normalize
v_brain <- normalizeGiotto(gobject = v_brain,
scalefactor = 6000,
verbose = TRUE)
# Add stats
v_brain <- addStatistics(gobject = v_brain)
# ID HVF
v_brain <- calculateHVF(gobject = v_brain,
method = "cov_loess")
fm <- fDataDT(v_brain)
hv_feats <- fm[hvf == "yes" & perc_cells > 3 & mean_expr_det > 0.4]$feat_ID
# Dimension Reductions
v_brain <- runPCA(gobject = v_brain,
feats_to_use = hv_feats)
v_brain <- runUMAP(v_brain,
dimensions_to_use = 1:10,
n_neighbors = 15,
set_seed = TRUE)
# NN Network
v_brain <- createNearestNetwork(gobject = v_brain,
dimensions_to_use = 1:10,
k = 15)
# Leiden Cluster
v_brain <- doLeidenCluster(gobject = v_brain,
resolution = 0.4,
n_iterations = 1000,
set_seed = TRUE)
# Spatial Network (kNN)
v_brain <- createSpatialNetwork(gobject = v_brain,
method = "kNN",
k = 5,
maximum_distance_knn = 400,
name = "spatial_network")
spatPlot2D(gobject = v_brain,
spat_unit = "cell",
cell_color = "leiden_clus",
show_image = TRUE,
point_size = 1.5,
point_shape = "no_border",
background_color = "black",
show_legend = TRUE,
save_plot = TRUE,
save_param = list(save_name = "03_ses6_1_vis_spat"))Here we can see the clustering of the regular visium spots is able to identify distinct regions of the mouse brain.
Figure 17.3: Mouse brain spatial plot showing leiden clustering
17.6.1 Identifying spatially organized features
We need to identify genes to be used for interpolation. This works best with genes that are spatially distinct. To identify these genes we’ll use binSpect(). For this tutorial we’ll only use the top 15 spatially distinct genes. The more genes used for interpolation the longer the analysis will take. When running this for your own datasets you should use more genes. We are only using 15 here to minimize analysis time.
# Spatially Variable Features
ranktest <- binSpect(v_brain,
bin_method = "rank",
calc_hub = TRUE,
hub_min_int = 5,
spatial_network_name = "spatial_network",
do_parallel = TRUE,
cores = 8) #not able to provide a seed number, so do not set one
# Getting the top 15 spatially organized genes
ext_spatial_features <- ranktest[1:15,]$feats17.7 Performing kriging
17.7.1 Interpolating features
Now we can perform gene expression interpolation. This involves creating a raster image for the gene expression of each of the selected genes. The steps from here can be time consuming and require large amounts of memory.
We will only be analyzing 15 genes to show the process of expression interpolation. For clustering and other analyses more genes are required.
future::plan(future::multisession()) # comment out for single threading
v_brain <- interpolateFeature(v_brain,
spat_unit = "cell",
feat_type = "rna",
ext = ext(v_brain),
feats = ext_spatial_features,
overwrite = TRUE)
print(v_brain)
Figure 17.4: Giotto object after to interpolating features. Addition of images for each interoplated feature (left) and an example of rasterized gene expression image (right).
For each gene that we interpolate a raster image is exported based on the gene expression. Shown below is an example of an output for the gene Pantr1.
Figure 17.5: Raster of gene expression interpolation for Pantr1
17.8 Adding cell polygons to Giotto object
17.8.1 Read in the poly information
First we need to read in the geojson file that contains the cell polygons that we’ll interpolate gene expression onto. These will then be added to the Giotto object as a new polygon object. This won’t affect the visium polygons. Both polygons will be stored within the same Giotto object.
# Read in the data
stardist_cell_poly_path <- file.path(data_directory, "segmentations/stardist_only_cell_bounds.geojson")
stardist_cell_gpoly <- createGiottoPolygonsFromGeoJSON(GeoJSON = stardist_cell_poly_path,
name = "stardist_cell",
calc_centroids = TRUE)
stardist_cell_gpoly <- flip(stardist_cell_gpoly)17.8.2 Vizualizing polygons
Below we can see a visualization of the polygons for the visium and the nuclei we identified from the H&E image. The visium dataset has 2698 spots compared to the 36694 nuclei we identified. Just using the visium spots we’re therefore losing a lot of the spatial data for individual cells. With the increased number of spots and them directly correlating with the tissue, through the spots alone we are able to better see the actual structure of the mouse brain.
Figure 17.6: Mouse brain cell polygons from the visium dataset
Figure 17.7: Mouse brain cell polygons with artifacts removed and flipped
17.8.3 Showing Giotto object prior to polygon addition
Before we add the polygons we can see the gobject contains “cell” as a spatial unit and a polygon.
Figure 17.8: Giotto object before adding subcellular polygons.
17.8.4 Adding polygons to giotto object
After we add the nuclei polygons we can see that a new polygon name, “stardist_cell” has been added to the gobject.
v_brain <- addGiottoPolygons(v_brain,
gpolygons = list("stardist_cell" = stardist_cell_gpoly))
print(v_brain)
Figure 17.9: Giotto object after to adding subcellular polygons.
17.8.5 Check polygon information
We can now see the addition of the new polygons under the name “stardist_cell”. Each of the new polyons is given a unique poly_ID as shown below. Each polygon is also added into same space as the original visium spots, therefore line up with the same image as the visium spots.
Figure 17.10: Polygon information for stardist_cell.
17.8.6 Expression overlap
The raster we created above gives the gene expression in a graphical form. We next need to determine how that relates to the nuclei location. To determine that we will calculate the overlap of the rasterized gene expression image to the polygons supplied earlier.
This step also takes more time the more genes that are provided. For large datasets please allow up to multiple hours for these steps to run.
v_brain <- calculateOverlapPolygonImages(gobject = v_brain,
name_overlap = "rna",
spatial_info = "stardist_cell",
image_names = ext_spatial_features)
v_brain <- Giotto::overlapToMatrix(x = v_brain,
poly_info = "stardist_cell",
feat_info = "rna",
aggr_function = "sum",
type="intensity")After performing the overlap we now have expression data for each gene provided. This can be seen below where we see the interpolated gene expression for genes in each of the nuclei we identified.
Figure 17.11: Gene expression for cells based on interpolation.
17.9 Reading in larger dataset
For better results more genes are required. The above data used only 15 genes. We will now read in a dataset that has 1500 interpolated genes an use this for the remained of the tutorial. If you haven’t downloaded this dataset please download it here.
17.10 Analyzing interpolated features
17.10.1 Filter and normalization
Now that we have a valid spat unit and gene expression data for each of the provided genes we can now perform the same analyses we used for the regular visium data. Please note that due to the differences in cell number that the values used for the current analysis aren’t identical to the visium analysis.
17.10.2 Visualizing gene expression from interpolated expression
Since we have the gene expression information for both the visium and the interpolated gene expression we can visualize gene expression for both from the same Giotto object. We will look at the expression for two genes “Sparc” and “Pantr1” for both the visium and interpolated data.
spatFeatPlot2D(v_brain,
spat_unit = "cell",
gradient_style = "sequential",
cell_color_gradient = "Geyser",
feats = "Sparc",
point_size = 2,
save_plot = TRUE,
show_image = TRUE,
save_param = list(save_name = "03_ses6_sparc_vis"))
spatFeatPlot2D(v_brain,
spat_unit = "stardist_cell",
gradient_style = "sequential",
cell_color_gradient = "Geyser",
feats = "Sparc",
point_size = 0.6,
save_plot = TRUE,
show_image = TRUE,
save_param = list(save_name = "03_ses6_sparc"))
spatFeatPlot2D(v_brain,
spat_unit = "cell",
gradient_style = "sequential",
feats = "Pantr1",
cell_color_gradient = "Geyser",
point_size = 2,
save_plot = TRUE,
show_image = TRUE,
save_param = list(save_name = "03_ses6_pantr1_vis"))
spatFeatPlot2D(v_brain,
spat_unit = "stardist_cell",
gradient_style = "sequential",
cell_color_gradient = "Geyser",
feats = "Pantr1",
point_size = 0.6,
save_plot = TRUE,
show_image = TRUE,
save_param = list(save_name = "03_ses6_pantr1"))Below we can see the gene expression for both datatypes. With the interpolated gene expression we’re able to get a better idea as to the cells that are expressing each of the genes. This is especially clear with Pantr1, which clearly localizes to the pyramidal layer.
Figure 17.12: Gene expression for visium (left) and interpolated (right) expression for Sparc (top) and Pantr1 (bottom).
17.10.4 Clustering
# UMAP
v_brain <- runUMAP(v_brain,
spat_unit = "stardist_cell",
dimensions_to_use = 1:15,
n_neighbors = 1000,
min_dist = 0.001,
spread = 1)
# NN Network
v_brain <- createNearestNetwork(gobject = v_brain,
spat_unit = "stardist_cell",
dimensions_to_use = 1:10,
feats_to_use = hv_feats,
expression_values = "normalized",
k = 70)
v_brain <- doLeidenCluster(gobject = v_brain,
spat_unit = "stardist_cell",
resolution = 0.15,
n_iterations = 100,
partition_type = "RBConfigurationVertexPartition")
plotUMAP(v_brain,
spat_unit = "stardist_cell",
cell_color = "leiden_clus")
Figure 17.13: UMAP for stardist_cell based on the 1500 interpolated gene expressions. Colored based on leiden clustering.
17.10.5 Visualizing clustering
Visualizing the clustering for both the visium dataset and the interpolated dataset we can get similar clusters. However, with the interpolated dataset we are able to see finer detail for each cluster.
spatPlot2D(gobject = v_brain,
spat_unit = "cell",
cell_color = "leiden_clus",
show_image = TRUE,
point_size = 0.5,
point_shape = "no_border",
background_color = "black",
save_plot = FALSE,
show_legend = TRUE)
spatPlot2D(gobject = v_brain,
spat_unit = "stardist_cell",
cell_color = "leiden_clus",
show_image = TRUE,
point_size = 0.1,
point_shape = "no_border",
background_color = "black",
show_legend = TRUE,
save_plot = TRUE,
save_param = list(save_name = "03_ses6_subcell_spat"))
Figure 17.14: Spatial plots showing leiden clustering mapped onto the base visium spots (left) and individual nuceli through interpolation (right)
17.10.6 Cropping objects
We are also able to crop both spat units simultaneously to zoom in on specific regions of the tissue such as seen below.
v_brain_crop <- subsetGiottoLocs(gobject = v_brain,
spat_unit = ":all:",
x_min = 4000,
x_max = 7000,
y_min = -6500,
y_max = -3500,
z_max = NULL,
z_min = NULL)
spatPlot2D(gobject = v_brain_crop,
spat_unit = "cell",
cell_color = "leiden_clus",
show_image = TRUE,
point_size = 2,
point_shape = "no_border",
background_color = "black",
show_legend = TRUE,
save_plot = TRUE,
save_param = list(save_name = "03_ses6_vis_spat_crop"))
spatPlot2D(gobject = v_brain_crop,
spat_unit = "stardist_cell",
cell_color = "leiden_clus",
show_image = TRUE,
point_size = 0.1,
point_shape = "no_border",
background_color = "black",
show_legend = TRUE,
save_plot = TRUE,
save_param = list(save_name = "03_ses6_subcell_spat_crop"))
Figure 17.15: Spatial plots showing leiden clustering mapped onto the base visium spots (left) and individual nuceli through interpolation (right)