## ----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)

## -----------------------------------------------------------------------------
p <- rMaxHeuristic(spe,
subsetby = "image_number", marks = "cell_type"
)
p

## ----crossSpatialInference, results='hide', message = FALSE, warning = FALSE----
#relevel to have non-diabetic as the reference category
spe$patient_stage <- relevel(factor(spe$patient_stage), "Non-diabetic")

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

names(resLs)

## ----bubblePlot, results='hide', message = FALSE, eval=TRUE-------------------
p <- plotCrossHeatmap(resLs, coefficientsToPlot = c("conditionLong_duration(x)", "conditionOnset(x)"), QCThreshold = 0, QCMetric = "medianMinIntensity")
p <- p + theme(legend.position = "bottom") + guides(shape = "none") + 
  facet_wrap(~factor(coefficient, levels = c("conditionOnset(x)", "conditionLong_duration(x)"), labels = c("conditionOnset(x)" = "Onset", "conditionLong_duration(x)" = "Long-Duration")))

p

## ----resSelection, eval=TRUE--------------------------------------------------
res <- resLs$alpha_Tc
names(res)

## ----plotFunction, message = FALSE, warning = FALSE---------------------------
#filtered and stabilised curves
metricRes <- res$metricRes
#raw spatial statistics curves
metricResRaw <- res$metricResRaw

# create a unique plotting ID
metricResRaw$ID <- paste0(
    metricResRaw$patient_stage, "|", metricResRaw$patient_id
)
# change levels for plotting
metricResRaw$ID <- factor(metricResRaw$ID, levels = c("Non-diabetic|6126",
                                                "Non-diabetic|6134",
                                                "Onset|6228","Onset|6414",
                                                "Long-duration|6089",
                                                "Long-duration|6180"))
# plot metrics
plotMetricPerFov(metricResRaw, 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-------------------------------------------------
residuals(mdl) |> cor() |> filled.contour(levels = seq(-1, 1, l = 40))
residuals(mdl) |> cov() |> filled.contour()

qqnorm(mdl$residuals, pch = 16)
qqline(mdl$residuals)

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