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

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

## ----setup, warning = FALSE, message = FALSE----------------------------------
library("spatialFDA")
library("dplyr")
library("ggplot2")
library("tidyr")
library("stringr")
library("dplyr")
library("patchwork")
library("SpatialExperiment")

set.seed(1234)

## ----loading, warning = FALSE, message = FALSE--------------------------------
# retrieve example data from Damond et al. (2019)
spe <- .loadExample(full = TRUE)

spe <- subset(spe, ,patient_id %in% c(6089,6180,6126,6134,6228,6414))
# set cell types as factors
colData(spe)$cell_type <- as.factor(colData(spe)$cell_type) 

## ----plotting fovs, warning = FALSE, fig.width=8, fig.height=15---------------
df <- data.frame(spatialCoords(spe), colData(spe))

dfSub <- df %>%
    subset(image_name %in% c("E02", "E03", "E04", "E05"))

p <- ggplot(dfSub, aes(x = cell_x, y = cell_y, color = cell_category)) +
    geom_point(size= 0.5) +
    facet_wrap(~image_name) +
    theme(legend.title.size = 20, legend.text.size = 20) +
    xlab("x") +
    ylab("y") +
    labs(color = "cell category")+
    coord_equal() +
    theme_light()

dfSub <- dfSub %>%
    subset(cell_type %in% c("alpha", "beta", "delta", "Th", "Tc"))

q <- ggplot(dfSub, aes(x = cell_x, y = cell_y, color = cell_type)) +
    geom_point(size= 0.5) +
    facet_wrap(~image_name) +
    theme(legend.title.size = 20, legend.text.size = 20) +
    xlab("x") +
    ylab("y") +
    labs(color = "cell type") +
    coord_equal() +
    theme_light()
wrap_plots(list(p,q), widths = c(1,1), heights = c(1,1), nrow = 2, ncol = 1)

## ----spatialInference---------------------------------------------------------
colData(spe)[["patient_stage"]] <- factor(colData(spe)[["patient_stage"]])
#relevel to have non-diabetic as the reference category
colData(spe)[["patient_stage"]] <- relevel(colData(spe)[["patient_stage"]],
"Non-diabetic")

#run the spatial statistics inference
res <- spatialInference(
    spe, 
    selection = c("alpha", "Tc"),
    subsetby = "image_number", 
    fun = "Gcross", 
    marks = "cell_type",
    rSeq = seq(0, 50, length.out = 50), 
    correction = "rs",
    sample_id = "patient_id",
    transformation = 0.5,
    eps = 0,
    delta = 10,
    family = mgcv::scat(link = "log"),
    image_id = "image_number", 
    condition = "patient_stage",
    ncores = 1
)

names(res)

## ----plotFunction-------------------------------------------------------------
metricRes <- res$metricRes

# create a unique plotting ID
metricRes$ID <- paste0(
    metricRes$patient_stage, "|", metricRes$patient_id
)
# change levels for plotting
metricRes$ID <- factor(metricRes$ID, levels = c("Non-diabetic|6126",
                                                "Non-diabetic|6134",
                                                "Onset|6228","Onset|6414",
                                                "Long-duration|6089",
                                                "Long-duration|6180"))
# plot metrics
plotMetricPerFov(metricRes, correction = "rs", x = "r",
                 imageId = "image_number", ID = "ID", ncol = 2)

## ----funcBoxPlot, warning = FALSE, results='hide'-----------------------------
# create a unique ID per row in the dataframe
metricRes$ID <- paste0(
    metricRes$patient_stage, "x", metricRes$patient_id,
    "x", metricRes$image_number
)

collector <- plotFbPlot(metricRes, "r", "rs", "patient_stage")

## ----funcGam------------------------------------------------------------------
mdl <- res$mdl
mm <- res$designmat

summary(mdl)

## ----plotFuncGam--------------------------------------------------------------
plotLs <- lapply(colnames(mm), plotMdl, mdl = mdl,
                 shift = mdl$coefficients[["(Intercept)"]])
wrap_plots(plotLs, nrow = 3, axes = 'collect')

## ----contour, warning = FALSE-------------------------------------------------
resid(mdl) |> cor() |> filled.contour(levels = seq(-1, 1, l = 40))
resid(mdl) |> cov() |> filled.contour()

qqnorm(resid(mdl), pch = 16)
qqline(resid(mdl))

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