## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=6,
  fig.height=4,
  dpi = 72,
  dev = "png"
)

## ----eval=FALSE---------------------------------------------------------------
# if (!require("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# 
# BiocManager::install("SEraster")

## ----eval=FALSE---------------------------------------------------------------
# require(remotes)
# remotes::install_github('JEFworks-Lab/SEraster')

## ----eval=FALSE---------------------------------------------------------------
# require(remotes)
# remotes::install_github('satijalab/seurat-wrappers@SEraster')

## ----message=FALSE------------------------------------------------------------
library(SEraster)
library(SpatialExperiment)
library(nnSVG)
library(CooccurrenceAffinity)
library(ggplot2)

## -----------------------------------------------------------------------------
data("merfish_mousePOA")

# check the dimension of the genes-by-cells matrix at single-cell resolution
dim(merfish_mousePOA)

# check the number of cell-types
length(unique(colData(merfish_mousePOA)$celltype))

## -----------------------------------------------------------------------------
# plot at single-cell resolution
df <- data.frame(spatialCoords(merfish_mousePOA), celltype = colData(merfish_mousePOA)$celltype)

ggplot(df, aes(x = x, y = y, col = celltype)) +
  coord_fixed() +
  geom_point(size = 1.5, stroke = 0) +
  guides(col = guide_legend(override.aes = list(size = 3))) +
  labs(x = "x (μm)",
       y = "y (μm)",
       col = "Cell-types") +
  theme_bw() +
  theme(panel.grid = element_blank())

## -----------------------------------------------------------------------------
rastGexp <- SEraster::rasterizeGeneExpression(merfish_mousePOA, assay_name="volnorm", resolution = 50)

# check the dimension of the genes-by-cells matrix after rasterizing gene expression
dim(rastGexp)

## -----------------------------------------------------------------------------
# plot total rasterized gene expression
SEraster::plotRaster(rastGexp, name = "Total rasterized gene expression")

## -----------------------------------------------------------------------------
# plot a specific gene
SEraster::plotRaster(rastGexp, feature_name = "Esr1", name = "Esr1")

## -----------------------------------------------------------------------------
# rasterize cell-type specific gene expression by subsetting to cell-type of interest
ct_interest <- "Inhibitory"
spe_subset <- merfish_mousePOA[,merfish_mousePOA$celltype == ct_interest]

# rasterize gene expression
rastGexpSubset <- SEraster::rasterizeGeneExpression(spe_subset, assay_name="volnorm", resolution = 50)

## -----------------------------------------------------------------------------
# plot
SEraster::plotRaster(rastGexpSubset, name = paste0("Total rast gexp in ", ct_interest))

## -----------------------------------------------------------------------------
SEraster::plotRaster(rastGexpSubset, feature_name = "Esr1", name = paste0("Esr1 in ", ct_interest))

## -----------------------------------------------------------------------------
rastCt <- SEraster::rasterizeCellType(merfish_mousePOA, col_name = "celltype", resolution = 50)

# check the dimension of the cell-types-by-cells matrix after rasterizing cell-type labels
dim(rastGexp)

## -----------------------------------------------------------------------------
# plot total cell counts
SEraster::plotRaster(rastCt, name = "cell counts", option = "inferno")

## -----------------------------------------------------------------------------
# plot specific cell-type
SEraster::plotRaster(rastCt, feature_name = "Inhibitory", name = "Inhibitory neuron counts", option = "inferno")

## -----------------------------------------------------------------------------
resolutions <- c(50, 100, 200)
for (resolution in resolutions) {
  # rasterize at defined resolution
  temp <- SEraster::rasterizeGeneExpression(merfish_mousePOA, assay_name="volnorm", resolution = resolution)
  # plot a specific gene
  plt <- SEraster::plotRaster(temp, feature_name = "Esr1", name = "Esr1", plotTitle = paste0("resolution: ", resolution))
  show(plt)
}

## -----------------------------------------------------------------------------
for (resolution in resolutions) {
  # rasterize at defined resolution
  temp <- SEraster::rasterizeGeneExpression(merfish_mousePOA, assay_name="volnorm", resolution = resolution, square = FALSE)
  # plot a specific gene
  plt <- SEraster::plotRaster(temp, feature_name = "Esr1", name = "Esr1", plotTitle = paste0("resolution: ", resolution))
  show(plt)
}

## -----------------------------------------------------------------------------
# permutate
spe_list <- permutateByRotation(merfish_mousePOA, n_perm = 3)

# rasterize permutated datasets at once
out_list <- rasterizeGeneExpression(spe_list, assay_name = "volnorm", resolution = 50)

for (i in seq_along(out_list)) {
  # extract rotated angle
  angle <- gsub("rotated_", "", paste0("rotated ", names(out_list)[[i]], " degrees"))
  # plot a specific gene
  plt <- SEraster::plotRaster(out_list[[i]], feature_name = "Esr1", name = "Esr1", plotTitle = angle)
  show(plt)
}

## -----------------------------------------------------------------------------
# run nnSVG
set.seed(0)
rastGexp <- nnSVG(rastGexp, assay_name = "pixelval")

## -----------------------------------------------------------------------------
# number of significant SVGs based on the selected adjusted p-value threshold
table(rowData(rastGexp)$padj <= 0.05)

## -----------------------------------------------------------------------------
# plot rasterized gene expression of top-ranked SVG
top_svg <- which(rowData(rastGexp)$rank == 1)
top_svg_name <- rownames(rowData(rastGexp))[top_svg]
SEraster::plotRaster(rastGexp, feature_name = top_svg_name, name = top_svg_name)

## -----------------------------------------------------------------------------
# subset data
ct_interest <- "Excitatory"
spe_sub <- merfish_mousePOA[,merfish_mousePOA$celltype == ct_interest]

# run SEraster
rastGexp_sub <- SEraster::rasterizeGeneExpression(spe_sub, assay_name="volnorm", resolution = 50)

# run nnSVG
set.seed(0)
rastGexp_sub <- nnSVG(rastGexp_sub, assay_name = "pixelval")

## -----------------------------------------------------------------------------
# number of significant SVGs
table(rowData(rastGexp_sub)$padj <= 0.05)

## -----------------------------------------------------------------------------
# plot rasterized gene expression of top-ranked SVG
top_svg <- which(rowData(rastGexp_sub)$rank == 1)
top_svg_name <- rownames(rowData(rastGexp_sub))[top_svg]
SEraster::plotRaster(rastGexp_sub, feature_name = top_svg_name, name = top_svg_name)

## -----------------------------------------------------------------------------
# extract cell-type labels
ct_labels <- as.factor(colData(merfish_mousePOA)$celltype)

# compute relative enrichment (RE) metric
mat <- assay(rastCt, "pixelval")
mat_re <- do.call(rbind, lapply(rownames(rastCt), function(ct_label) {
    mat[ct_label,] / (sum(mat[ct_label,]) / sum(mat) * colSums(mat))
}))
rownames(mat_re) <- rownames(mat)

# binarize
mat_bin <- ifelse(mat_re >= 1, 1, 0)

# add RE and binarized layers to SpatialExperiment object
assays(rastCt) <- list(pixelval = assay(rastCt, "pixelval"), re = mat_re, bin = mat_bin)

## -----------------------------------------------------------------------------
ct_interest <- "Ependymal"

# plot pixel value for a cell-type of interest
plotRaster(rastCt, assay_name = "pixelval", feature_name = ct_interest, name = "cell-type counts", option = "inferno")

## -----------------------------------------------------------------------------
# plot RE value for a cell-type of interest
plotRaster(rastCt, assay_name = "re", feature_name = ct_interest, name = "RE", option = "inferno")

## -----------------------------------------------------------------------------
# plot binarized value for a cell-type of interest
plotRaster(rastCt, assay_name = "bin", feature_name = ct_interest, factor_levels = c(0,1), name = "binarized", option = "inferno")

## -----------------------------------------------------------------------------
# run CooccurrenceAffinity
ct_coocc <- CooccurrenceAffinity::affinity(data = mat_bin, row.or.col = "row", squarematrix = c("all"))

# plot maximum likelihood estimates of affinity metric (alpha MLE)
CooccurrenceAffinity::plotgg(data = ct_coocc, variable = "alpha_mle", legendlimit = "datarange")

## -----------------------------------------------------------------------------
sessionInfo()