## ----setup, include=FALSE-----------------------------------------------------
options(width=80)
knitr::opts_chunk$set(collapse=TRUE,
                      message=FALSE,
                      warning=FALSE,
                      comment="",
                      fig.align="center",
                      fig.wide=TRUE)

## ----message=FALSE, warning=FALSE---------------------------------------------
library(SpatialExperiment)
library(ExperimentHub)
library(spatialLIBD)

ehub <- ExperimentHub()
spe <- fetch_data(type="spe", eh=ehub)

## -----------------------------------------------------------------------------
table(spe$sample_id)

## -----------------------------------------------------------------------------
spe <- spe[, spe$sample_id == "151673"]
dim(spe)

## -----------------------------------------------------------------------------
stopifnot(all(colData(spe)$in_tissue))

## -----------------------------------------------------------------------------
table(colData(spe)$spatialLIBD, useNA="always")
spe <- spe[, !is.na(colData(spe)$spatialLIBD)]
dim(spe)
table(colData(spe)$spatialLIBD, useNA="always")

## ----message=FALSE, warning=FALSE---------------------------------------------
library(scrapper)

is_mito <- grepl("^MT-", rowData(spe)$gene_name)
table(is_mito)

## -----------------------------------------------------------------------------
spe <- quickRnaQc.se(spe, subsets=list(mito=is_mito))
spe$discard <- !spe$keep
table(spe$discard)

## ----dlpfcqcplot, echo=TRUE, message=FALSE, fig.height=5, fig.width=5, out.width="600px", fig.cap="Detected low-quality spots in the DLPFC sample 151673."----
library(ggspavis)

plotObsQC(spe, plot_type="spot", annotate="discard")

## -----------------------------------------------------------------------------
spe <- spe[, spe$keep]
dim(spe)

## -----------------------------------------------------------------------------
spotsxgene <- rowSums(counts(spe) > 0)
spe <- spe[spotsxgene > floor(ncol(spe)*0.01), ]
dim(spe)

## -----------------------------------------------------------------------------
spe <- normalizeRnaCounts.se(spe, size.factors=spe$sum)
assayNames(spe)

## ----message=FALSE, warning=FALSE---------------------------------------------
library(GSEABase)
library(GSVA)

fname <- file.path(system.file("extdata", package="GSVAdata"),
                   "human_brain_snRNAseq_cellType_markers.gmt.gz")
gsets <- readGMT(fname)
gsets

## -----------------------------------------------------------------------------
gsvapar <- gsvaParam(spe, gsets)
gsvapar

## -----------------------------------------------------------------------------
es <- gsva(gsvapar)
es

## -----------------------------------------------------------------------------
whmax <- apply(assay(es), 2, which.max)
es$cellType <- factor(rownames(es)[whmax])

## ----dlpfcgsva, echo=TRUE, message=FALSE, fig.height=10, fig.width=7, out.width="800px", dpi=100, fig.cap="(a) Annotated layers of the DLPFC sample 151673 and (b) spatial distribution of the most likely cell type in each spot, based on GSVA scores."----
library(ggspavis)
library(patchwork)
library(RColorBrewer)

ctpal <- colorRampPalette(brewer.pal(9, "Set1"))(nlevels(es$cellType))
plts <- list(plotVisium(spe, point_size=0.60, facets=NULL, annotate="spatialLIBD",
                        pal=brewer.pal(nlevels(spe$spatialLIBD), "Set1")) +
             labs(tag="a"),
             plotVisium(es, assay="es", point_size=0.60, facets=NULL,
                        annotate="cellType",
                        pal=ctpal) +
             labs(tag="b")
             )
wrap_plots(plts, nrow=1)

## -----------------------------------------------------------------------------
I <- spatCor(es)
I[order(I$observed, decreasing=TRUE), ]

## ----dlpfcgsvaoligo, echo=TRUE, message=FALSE, fig.height=10, fig.width=7, out.width="800px", dpi=100, fig.cap="Spatial distribution of the GSVA scores for the Oligodendrocytes gene set in the DLPFC sample 151673."----
plotVisium(es, assay="es", point_size=1, facets=NULL, annotate="Oligo")

## ----session_info, cache=FALSE------------------------------------------------
sessionInfo()