## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)
set.seed(8048)

## ----installation, include = TRUE, eval=FALSE---------------------------------
# if (!requireNamespace("BiocManager")) {
#     install.packages("BiocManager")
# }
# BiocManager::install("sosta")

## ----setup, message=FALSE, warning=FALSE--------------------------------------
library("dplyr")
library("ExperimentHub")
library("ggplot2")
library("lme4")
library("lmerTest")
library("sf")
library("SpatialExperiment")
library("sosta")
library("tidyr")
library("ggfortify")

theme_set(theme_bw())

## ----loading, echo=FALSE, message=FALSE---------------------------------------
# load experiment hub
eh <- ExperimentHub()
oid <- names(eh[eh$title == "Damond_2019_Pancreas - sce - v1 - full"])

# Load single cell experiment object
spe <- eh[[oid]]
# Convert to spatial experiment object
spe <- toSpatialExperiment(spe,
    sample_id = "image_name",
    spatialCoordsNames = c("cell_x", "cell_y")
)
# We don't need expression, so we remove it to free up memory
assays(spe) <- list()
gc()

## ----plot-data----------------------------------------------------------------
df <- cbind(
    colData(spe[, spe$image_name %in% c("E04", "E03", "G01", "J02")]),
    spatialCoords(spe[, spe$image_name %in% c("E04", "E03", "G01", "J02")])
)
df |>
    as.data.frame() |>
    ggplot(aes(x = cell_x, y = cell_y, color = cell_category)) +
    geom_point(size = 0.25) +
    facet_wrap(~image_name, ncol = 2) +
    coord_equal() +
    scale_color_brewer(palette = "Dark2")

## -----------------------------------------------------------------------------
shapeIntensityImage(
    spe,
    marks = "cell_category",
    imageCol = "image_name",
    imageId = "G01",
    markSelect = "islet"
)

## -----------------------------------------------------------------------------
n <- estimateReconstructionParametersSPE(
    spe,
    marks = "cell_category",
    imageCol = "image_name",
    markSelect = "islet",
    nImages = 25,
    plotHist = TRUE
)

## -----------------------------------------------------------------------------
n |>
    ggplot(aes(x = bndw, y = thres)) +
    geom_point()

## -----------------------------------------------------------------------------
(thresSPE <- mean(n$thres))
(bndwSPE <- mean(n$bndw))

## ----eval=TRUE, include=TRUE--------------------------------------------------
# Sample 15 images per patient
sel <- colData(spe) |>
    as.data.frame() |>
    group_by(patient_id) |>
    select(image_name) |>
    sample_n(size = 20, replace = FALSE) |>
    pull(image_name)

# Select sampled images
speSel <- spe[, spe$image_name %in% sel]
speSel$image_name |>
    unique() |>
    length()

## ----eval=TRUE----------------------------------------------------------------
# Run on all images
allIslets <- reconstructShapeDensitySPE(
    speSel,
    marks = "cell_category",
    imageCol = "image_name",
    markSelect = "islet",
    bndw = bndwSPE,
    thres = thresSPE,
    nCores = 1
)

## -----------------------------------------------------------------------------
colsKeep <- c(
    "patient_stage", "tissue_slide", "tissue_region",
    "patient_id", "patient_disease_duration",
    "patient_age", "patient_gender", "sample_id"
)

patientData <- colData(speSel) |>
    as_tibble() |>
    group_by(image_name) |>
    select(all_of(colsKeep)) |>
    unique()

allIslets <- allIslets |>
    dplyr::left_join(patientData, by = "image_name")

## -----------------------------------------------------------------------------
allIslets |>
    st_drop_geometry() |> # we are only interested in metadata
    group_by(patient_id) |>
    summarise(n = n()) |>
    ungroup()

## -----------------------------------------------------------------------------
isletMetrics <- totalShapeMetrics(allIslets)

## -----------------------------------------------------------------------------
# specify factor levels
lv <- c("Non-diabetic", "Onset", "Long-duration")

allIslets <- allIslets |>
    cbind(t(isletMetrics)) |>
    mutate(patient_stage = factor(patient_stage, levels = lv))

## -----------------------------------------------------------------------------
autoplot(
    prcomp(t(isletMetrics), scale. = TRUE),
    x = 1,
    y = 2,
    data = allIslets,
    color = "patient_stage",
    size = 2,
    # shape = 'type',
    loadings = TRUE,
    loadings.colour = "steelblue3",
    loadings.label = TRUE,
    loadings.label.size = 3,
    loadings.label.repel = TRUE,
    loadings.label.colour = "black"
) +
    scale_color_brewer(palette = "Set2") +
    theme_bw() +
    coord_fixed()

## ----fig.width=8, fig.height=6------------------------------------------------
allIslets |>
    sf::st_drop_geometry() |>
    select(patient_stage, rownames(isletMetrics)) |>
    pivot_longer(-patient_stage) |>
    filter(name %in% c("Area", "Compactness", "Curl")) |>
    ggplot(aes(x = patient_stage, y = value, fill = patient_stage)) +
    geom_boxplot() +
    facet_wrap(~name, scales = "free") +
    scale_fill_brewer(palette = "Set2") +
    scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
    guides(fill = "none")

## ----fig.width=8, fig.height=6------------------------------------------------
allIslets |>
    sf::st_drop_geometry() |>
    select(patient_stage, patient_id, rownames(isletMetrics)) |>
    pivot_longer(-c(patient_stage, patient_id)) |>
    filter(name %in% c("Area")) |>
    ggplot(aes(
        x = as.factor(patient_id),
        y = log(value),
        fill = patient_stage
    )) +
    geom_boxplot() +
    facet_wrap(~patient_stage, scales = "free") +
    scale_fill_brewer(palette = "Set2") +
    scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
    guides(fill = "none")

## -----------------------------------------------------------------------------
mod <- lmer(
    log(Area) ~ patient_stage +
        (1 | patient_id) + (1 | image_name),
    data = allIslets
)

## -----------------------------------------------------------------------------
summary(mod)

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