## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  comment = "#>"
)

# load libraries without messages
library(SpaNorm)
library(ggplot2)
library(patchwork)
library(SpatialExperiment)
library(scater)

update_geom_defaults("point", aes(size = 0.5))

## ----eval=FALSE---------------------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# 
# # release version
# BiocManager::install("SpaNorm")
# # development version from GitHub
# BiocManager::install("bhuvad/SpaNorm")

## ----fig.width=4, fig.height=4.25---------------------------------------------
library(SpaNorm)
library(SpatialExperiment)
library(ggplot2)

# load sample data
data(HumanDLPFC)
# change gene IDs to gene names
rownames(HumanDLPFC) = rowData(HumanDLPFC)$gene_name
HumanDLPFC

# plot regions
p_region = plotSpatial(HumanDLPFC, colour = AnnotatedCluster, size = 0.5) +
  scale_colour_brewer(palette = "Paired", guide = guide_legend(override.aes = list(shape = 15, size = 5))) +
  ggtitle("Region")
p_region

## -----------------------------------------------------------------------------
# filter genes expressed in 20% of spots
keep = filterGenes(HumanDLPFC, 0.2)
table(keep)
# subset genes
HumanDLPFC = HumanDLPFC[keep, ]

## ----fig.width=7.5, fig.height=4.25-------------------------------------------
logcounts(HumanDLPFC) = log2(counts(HumanDLPFC) + 1)

p_counts = plotSpatial(
    HumanDLPFC,
    colour = MOBP,
    what = "expression",
    assay = "logcounts",
    size = 0.5
  ) +
  scale_colour_viridis_c(option = "F") +
  ggtitle("logCounts")
p_region + p_counts

## ----message=TRUE-------------------------------------------------------------
set.seed(36)
HumanDLPFC = SpaNorm(HumanDLPFC)
HumanDLPFC

## ----fig.width=7.5, fig.height=4.25-------------------------------------------
p_logpac = plotSpatial(
    HumanDLPFC,
    colour = MOBP,
    what = "expression",
    assay = "logcounts",
    size = 0.5
  ) +
  scale_colour_viridis_c(option = "F") +
  ggtitle("logPAC")
p_region + p_logpac

## -----------------------------------------------------------------------------
library(Seurat)

# create a Seurat object with the data
seurat_obj = HumanDLPFC
logcounts(seurat_obj) = counts(seurat_obj)
seurat_obj = suppressWarnings(Seurat::as.Seurat(seurat_obj))  
# add spatial coordinates to Seurat meta.data from SpatialExperiment
coords = SpatialExperiment::spatialCoords(HumanDLPFC)
seurat_obj@meta.data$x = coords[, 1]
seurat_obj@meta.data$y = coords[, 2]

# run SpaNorm on Seurat (CPU backend)
seurat_obj <- SpaNorm(seurat_obj, sample.p = 0.1, df.tps = 2, backend = "cpu", verbose = TRUE)
seurat_obj

## -----------------------------------------------------------------------------
# manually retrieve model
fit.spanorm = metadata(HumanDLPFC)$SpaNorm
fit.spanorm

## ----fig.width=11.5, fig.height=8.5-------------------------------------------
# Pearson residuals
HumanDLPFC = SpaNorm(HumanDLPFC, adj.method = "pearson")
p_pearson = plotSpatial(
    HumanDLPFC,
    colour = MOBP,
    what = "expression",
    assay = "logcounts",
    size = 0.5
  ) +
  scale_colour_viridis_c(option = "F") +
  ggtitle("Pearson")

# meanbio residuals
HumanDLPFC = SpaNorm(HumanDLPFC, adj.method = "meanbio")
p_meanbio = plotSpatial(
    HumanDLPFC,
    colour = MOBP,
    what = "expression",
    assay = "logcounts",
    size = 0.5
  ) +
  scale_colour_viridis_c(option = "F") +
  ggtitle("Mean biology")

# meanbio residuals
HumanDLPFC = SpaNorm(HumanDLPFC, adj.method = "medbio")
p_medbio = plotSpatial(
    HumanDLPFC,
    colour = MOBP,
    what = "expression",
    assay = "logcounts",
    size = 0.5
  ) +
  scale_colour_viridis_c(option = "F") +
  ggtitle("Median biology")

p_region + p_counts + p_logpac + p_pearson + p_meanbio + p_medbio + plot_layout(ncol = 3)

## ----fig.width=7.5, fig.height=4.25-------------------------------------------
# df.tps = 2
HumanDLPFC_df2 = SpaNorm(HumanDLPFC, df.tps = 2)
p_logpac_2 = plotSpatial(
    HumanDLPFC,
    colour = MOBP,
    what = "expression",
    assay = "logcounts",
    size = 0.5
  ) +
  scale_colour_viridis_c(option = "F") +
  ggtitle("logPAC (df.tps = 2)")

# df.tps = 6 (default)
p_logpac_6 = p_logpac +
  ggtitle("logPAC (df.tps = 6)")

p_logpac_2 + p_logpac_6

## ----fig.width=7.5, fig.height=4.25-------------------------------------------
# scale.factor = 1 (default)
HumanDLPFC = SpaNorm(HumanDLPFC, scale.factor = 1)
p_logpac_sf1 = plotSpatial(
    HumanDLPFC,
    colour = MOBP,
    what = "expression",
    assay = "logcounts",
    size = 0.5
  ) +
  scale_colour_viridis_c(option = "F") +
  ggtitle("logPAC (scale.factor = 1)")

# scale.factor = 4
HumanDLPFC = SpaNorm(HumanDLPFC, scale.factor = 4)
p_logpac_sf4 = plotSpatial(
    HumanDLPFC,
    colour = MOBP,
    what = "expression",
    assay = "logcounts",
    size = 0.5
  ) +
  scale_colour_viridis_c(option = "F") +
  ggtitle("logPAC (scale.factor = 4)")

p_logpac_sf1 + p_logpac_sf4 + plot_layout(ncol = 2)

## ----fig.width=7.5, fig.height=4.25-------------------------------------------
p1 = plotCovariate(HumanDLPFC, colour = MOBP, covariate = "biology") +
  scale_colour_viridis_c(option = "F") +
  ggtitle("Biology")
p2 = plotCovariate(HumanDLPFC, colour = MOBP, covariate = "ls") +
  scale_colour_viridis_c(option = "F") +
  ggtitle("Library size effect")
p1 + p2

## ----message=TRUE-------------------------------------------------------------
HumanDLPFC = SpaNormSVG(HumanDLPFC)
HumanDLPFC

## -----------------------------------------------------------------------------
svgs = topSVGs(HumanDLPFC, n = 10)
svgs

## ----fig.width=12, fig.height=12----------------------------------------------
# fix gene names
rownames(HumanDLPFC) = gsub("-", ".", rownames(HumanDLPFC))
rownames(svgs) = gsub("-", ".", rownames(svgs))

lapply(rownames(svgs)[1:9], function(g) {
  plotSpatial(HumanDLPFC, colour = !!sym(g), what = "expression", assay = "logcounts", size = 0.5) +
    scale_colour_viridis_c(option = "F") +
    ggtitle(g) +
    theme(legend.position = "bottom")
}) |> 
  wrap_plots(ncol = 3)

## ----fig.width=6, fig.height=5------------------------------------------------
library(scater)

HumanDLPFC = SpaNormPCA(HumanDLPFC, ncomponents = 50, svg.fdr = 0.2, nsvgs = Inf)
HumanDLPFC = runUMAP(HumanDLPFC, n_neighbors = 20, min_dist = 0.3)
plotUMAP(HumanDLPFC, colour_by = "AnnotatedCluster", size_by = "cell_count") +
  scale_colour_brewer(palette = "Paired", guide = guide_legend(override.aes = list(shape = 15, size = 5))) +
  labs(title = "UMAP derived from SpaNorm PCA", colour = "Cluster")

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