Skip to content

Commit

Permalink
fixes for vignette
Browse files Browse the repository at this point in the history
  • Loading branch information
zskylarli committed May 29, 2024
1 parent e207a9c commit 7530040
Show file tree
Hide file tree
Showing 3 changed files with 198 additions and 11 deletions.
3 changes: 3 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,9 @@ export(RunQuantileNorm)
export(RunSNF)
export(RunVelocity)
export(Runtricycle)
export(rasterizeGeneExpression)
export(rasterizeCellType)
export(permutateByRotation)
export(StopCellbrowser)
export(as.cell_data_set)
export(scVIIntegration)
Expand Down
23 changes: 12 additions & 11 deletions R/SEraster.R
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ createRasterizedObject <- function(input, out, name) {
}
)

output_radius <- sqrt(nrow(input[[]]) / nrow(meta_rast)) * ifelse(!is.null(scale_factors), Radius(input_fov, scale = NULL), input_fov[['centroids']]@radius)
output_radius <- sqrt(nrow(input[[]]) / nrow(meta_rast)) * ifelse(!is.null(Radius(input_fov, scale = NULL)), Radius(input_fov, scale = NULL), Radius(input_fov[['centroids']], scale = NULL))

output_centroids <- CreateCentroids(
coords = output_coordinates,
Expand Down Expand Up @@ -108,7 +108,7 @@ rasterizeGeneExpression <- function(
image <- ifelse(is.null(image), Images(spe, assay=assay_name)[1], image)
coords <- GetTissueCoordinates(spe, image = image, scale = NULL)
if("cell" %in% colnames(coords)){
rownames(coords) <- coords$cell
rownames(coords) <- colnames(spe)
}

if (is.null(assay_name)) {
Expand Down Expand Up @@ -243,7 +243,7 @@ rasterizeCellType <- function(
image <- ifelse(is.null(image), Images(input, assay=assay_name)[1], image)
pos <- GetTissueCoordinates(input, image = image, scale = NULL)
if("cell" %in% colnames(pos)){
rownames(coords) <- pos$cell
rownames(pos) <- pos$cell
}
bbox <- sf::st_bbox(c(
xmin = floor(min(pos[,1])-resolution/2),
Expand Down Expand Up @@ -321,16 +321,15 @@ permutateByRotation <- function(input, n_perm = 1, verbose = FALSE) {
class <- class(spe[[image_name]])

if (class == 'VisiumV1') {
input <-updateVisiumV1(spe, pos_rotated, assay_name, angle, image_name)
input <- updateVisiumV1(spe, pos_rotated, assay_name, angle, image_name)
} else if (class == 'VisiumV2') {
input <-updateVisiumV2(spe, pos_rotated, assay_name, angle, image_name)
input <- updateVisiumV2(spe, pos_rotated, assay_name, angle, image_name)
} else if (class == 'FOV') {
input <-updateFOV(spe, pos_rotated, assay_name, angle, image_name)
input <- updateFOV(spe, pos_rotated, assay_name, angle, image_name)
}
})
}))

#names(output) <- ifelse(!is.null(names(input)), paste0(rep(names(input), each = length(angles)), "_rotated_", angles), NULL)
return(output)

} else {
Expand All @@ -339,6 +338,7 @@ permutateByRotation <- function(input, n_perm = 1, verbose = FALSE) {
colnames(pos_orig) <- c("x", "y")
stopifnot("Column 1 and 2 of the spatialCoords slot should be named x and y, respectively." = colnames(pos_orig)[1:2] == c("x", "y"))

output_list <- list()
for (angle in angles) {
image_name <- Images(input, assay = assay_name)[[1]]
image <- input[[image_name]]
Expand All @@ -350,14 +350,15 @@ permutateByRotation <- function(input, n_perm = 1, verbose = FALSE) {
class <- class(input[[image_name]])

if (class == 'VisiumV1') {
input <- updateVisiumV1(input, pos_rotated, assay_name, angle, image_name)
output <- updateVisiumV1(input, pos_rotated, assay_name, angle, image_name)
} else if (class == 'VisiumV2') {
input <-updateVisiumV2(input, pos_rotated, assay_name, angle, image_name)
output <- updateVisiumV2(input, pos_rotated, assay_name, angle, image_name)
} else if (class == 'FOV') {
input <-updateFOV(input, pos_rotated, assay_name, angle, image_name)
output <- updateFOV(input, pos_rotated, assay_name, angle, image_name)
}
output_list <- append(output_list, output)
}
return(input)
return(output_list)
}
}

Expand Down
183 changes: 183 additions & 0 deletions docs/SEraster.Rmd
Original file line number Diff line number Diff line change
@@ -0,0 +1,183 @@
---
title: "Getting Started With SEraster"
date: 'Compiled: `r format(Sys.Date(), "%B %d, %Y")`'
output:
github_document:
html_preview: true
toc: true
html_document:
df_print: kable
theme: simplex
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
message = FALSE,
warning = FALSE
)
```

# Getting Started with SEraster

This tutorial walks you through the basic functionalities of `SEraster` and two examples of downstream analysis that can be performed with the rasterized spatial omics data.

In the examples below, we assume the input data is provided as a `SpatialExperiment` Bioconductor object. Please refer to the following documentations to see how you would format your data into a `SpatialExperiment` object:

- [SpatialExperiment](https://bioconductor.org/packages/SpatialExperiment) package
- [Formatting a SpatialExperiment Object for SEraster](https://jef.works/SEraster/articles/formatting-SpatialExperiment-for-SEraster.html)
- [`merfish_mousePOA` dataset](https://jef.works/SEraster/reference/merfish_mousePOA.html)


## Load libraries

```{r}
library(Seurat)
#library(SeuratWrappers)
devtools::load_all()
library(SeuratData)
library(SEraster)
library(ggplot2)
```

## Load example dataset

This vignette uses the Tiny subset dataset provided in the Fresh Frozen Mouse Brain for Xenium Explorer Demo as well as the VisiumHD mouse brain dataset, both from from 10x Genomics. When using Seurat, we use different plotting functions for image-based and sequencing-based spatial datasets, and will therefore demonstrate SEraster behavior using both examples.

```{r}
path <- "/brahms/hartmana/vignette_data/xenium_tiny_subset"
# Load the Xenium data
xenium.obj <- LoadXenium(path, fov = "fov")
# remove cells with 0 counts
xenium.obj <- subset(xenium.obj, subset = nCount_Xenium > 0)
```
```{r}
# visualize
ImageDimPlot(xenium.obj, fov = "fov")
```
```{r}
localdir <- "/brahms/lis/visium_hd/mouse/new_mousebrain/"
# Load the VisiumHD data
visiumhd.obj <- Load10X_Spatial(data.dir = localdir, bin.size = 16)
```
```{r}
# visualize
SpatialDimPlot(visiumhd.obj, image = 'slice1.016um')
```

## SEraster basic functionalities

`SEraster` reduces the number of spatial points in spatial omics datasets for downstream analysis through a process of rasterization where single cells’ gene expression or cell-type labels are aggregated into equally sized square or hexagonal pixels (can be changed using the `square` argument) based on a user-defined resolution.

Here, we demonstrate the basic functionalities of `SEraster`.

### Rasterize gene expression

For continuous variables such as gene expression or other molecular information, `SEraster` aggregates the observed raw counts or normalized expression values for each molecule within each pixel using means by default (can be changed using the `fun` argument). When rasterizing Seurat objects, users can specify the `assay_name`, `slot`, and `image` for spatial objects with associated images.

Let's try rasterizing the gene expression of the Xenium mouse brain dataset we loaded.

```{r}
rastGexp <- SeuratWrappers::rasterizeGeneExpression(xenium.obj, assay_name = "Xenium", resolution = 100)
# check dimensions of spot by gene matrix after rasterization
dim(rastGexp[['Xenium']])
```

As you can see, `SEraster` aggregated 6,509 single cells into 1,301 pixels.

```{r}
# plot total rasterized gene expression
ImageDimPlot(rastGexp, size = 5)
```

```{r}
# plot specific genes on rasterized object
ImageFeaturePlot(rastGexp, size = 2.5, features = c("Cux2", "Rorb", "Bcl11b", "Foxp2"), max.cutoff = c(25, 35, 12, 10), cols = c("white", "red"))
# compare to original object
ImageFeaturePlot(xenium.obj, features = c("Cux2", "Rorb", "Bcl11b", "Foxp2"), max.cutoff = c(25, 35, 12, 10), size = 0.75, cols = c("white", "red"))
```

### Rasterize cell-type

For categorical variables such as cell-type or cluster labels, `SEraster` aggregates the number of cells for each label within each pixel using sums by default (can be changed using the `fun` argument).

Let's try rasterizing the cluster labels of the Xenium Tiny dataset.

```{r}
xenium.obj <- SCTransform(xenium.obj, assay = "Xenium")
xenium.obj <- RunPCA(xenium.obj, npcs = 30, features = rownames(xenium.obj))
xenium.obj <- RunUMAP(xenium.obj, dims = 1:30)
xenium.obj <- FindNeighbors(xenium.obj, reduction = "pca", dims = 1:30)
xenium.obj <- FindClusters(xenium.obj, resolution = 0.3)
```
```{r}
DimPlot(xenium.obj, group.by='seurat_clusters')
ImageDimPlot(xenium.obj, group.by='seurat_clusters')
```

```{r}
rastCt <- SeuratWrappers::rasterizeCellType(xenium.obj, col_name = "seurat_clusters", resolution = 200)
# check the dimension of the cell-types-by-cells matrix after rasterizing cell-type labels
dim(rastCt[['seurat_clusters']])
```

```{r}
# plot total cell counts per rasterized spot
ImageFeaturePlot(rastCt, size = 7.5, features = 'num_cell')
# plot number of seurat clusters for each rasterized spot
ImageFeaturePlot(rastCt, size = 7.5, features = 'nFeature_seurat_clusters')
```

### Setting rasterization resolution

Rasterization resolution can be controlled by the `resolution` argument of the `rasterizeGeneExpression` and `rasterizeCellType` functions. Here, we refer to a particular resolution of rasterization by the side length for square pixels and the distance between opposite edges for hexagonal pixels such that finer resolution indicates smaller pixel size and vice versa.

Let's see how the rasterized VisiumHD mouse brain dataset looks with various resolutions using square pixels, which is specified by `shape`.

```{r}
# rasterize at defined resolution
rast2000 <- SeuratWrappers::rasterizeGeneExpression(visiumhd.obj, assay_name="Spatial.016um", resolution = 2000)
rast2000 <- NormalizeData(rast2000)
p1 <- SpatialFeaturePlot(rast2000, shape = 22, pt.size.factor = 2, features = "Hpca") + ggtitle("Hpca expression (res = 2000)")
rast1000 <- SeuratWrappers::rasterizeGeneExpression(visiumhd.obj, assay_name="Spatial.016um", resolution = 1000)
rast1000 <- NormalizeData(rast1000)
p2 <- SpatialFeaturePlot(rast1000, shape = 22, pt.size.factor = 2, features = "Hpca") + ggtitle("Hpca expression (res = 1000)")
p1 | p2
```

### Creating and rasterizing permutations

Since rasterized values may be sensitive to edge effects such as the specific boundaries of grids upon rasterization, `SEraster` enables permutation by rotating the dataset at various angles before rasterization, which has been enabled in the SeuratWrappers implementation.

For example, we create 3 permutations of the Xenium Tiny mouse brain dataset, which would return the object with x,y coordinates rotated at 0, 120, and 240 degrees around the midrange point as new field of views.

In addition to a single `Seurat` spatial object, `rasterizeGeneExpression` and `rasterizeCellType` functions can both take a `list` of `Seurat` objects. This essentially allows users to streamline the preprocessing of permutations with `SEraster`; followed by a downstream analysis of choice. When`SEraster` rasterizes a `list` of objects, all objects in the inputted `list` are rasterized with the same pixel coordinate framework (same bounding box, resolution, centroid coordinates). This feature may not be particularly useful for permutations; however, it can potentially be applied to compare two or more datasets, such as structurally aligned tissues as well as healthy vs. disease tissues.

```{r}
# permutate
DefaultAssay(xenium.obj) <- "Xenium"
spe_list <- SeuratWrappers::permutateByRotation(xenium.obj, n_perm = 3)
# rasterize permutated datasets at once
out_list <- SeuratWrappers::rasterizeGeneExpression(spe_list, assay_name = "Xenium", resolution = 200)
for (i in seq_along(out_list)) {
# extract rotated angle
angle <- Images(out_list[[i]])
# plot a specific gene
plt <- ImageFeaturePlot(out_list[[i]], features = "Rorb", max.cutoff = 35, size = 7.5, cols = c("white", "red")) + ggtitle(angle)
show(plt)
}
```

<details>

<summary>**Session Info**</summary>

```{r}
sessionInfo()
```

</details>

0 comments on commit 7530040

Please sign in to comment.