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

## ----Lfunction, warning = FALSE, message = FALSE------------------------------
metricRes <- calcMetricPerFov(spe = spe, selection = c("alpha", "Tc"),
                              subsetby = "image_number", fun = "Lcross", 
                              marks = "cell_type",
                              rSeq = seq(0, 50, length.out = 50), 
                              by = c("patient_stage", "patient_id",
                                     "image_number"),
                              ncores = 1)
metricRes %>% head(3)

## ----plotLfunction, warning = FALSE, fig.width=8, fig.height=8----------------
# 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 = "iso", x = "r",
                 imageId = "image_number", ID = "ID", ncol = 2)

## ----Gfunction, warning = FALSE, message = FALSE------------------------------
metricRes <- calcMetricPerFov(spe = spe, selection = c("alpha", "Tc"),
                              subsetby = "image_number", fun = "Gcross", 
                              marks = "cell_type",
                              rSeq = seq(0, 50, length.out = 50), 
                              by = c("patient_stage", "patient_id",
                                     "image_number"),
                              ncores = 1)
metricRes %>% head(3)

## ----plotGfunction, warning = FALSE, fig.width=8, fig.height=8----------------
# 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
)
#removing field of views that have as a curve only zeros - these are cases where
#there is no cells of one type
metricRes <- metricRes %>% dplyr::group_by(ID) %>% dplyr::filter(sum(rs) >= 1)

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

## ----variancetransform, warning = FALSE---------------------------------------
# can determine with a boxcox transformation what is the ideal parameter
# for the transformation use Fisher's variance-stabilising transformation
metricRes$rs <- asin(sqrt(metricRes$rs))

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

## ----fPCA, warning = FALSE----------------------------------------------------
# filter out all rows that have a constant zero part - all r=<10
metricRes <- metricRes %>% filter(r >= 10)

# prepare dataframe from calcMetricRes to be in the correct format for pffr
dat <- prepData(metricRes, "r", "rs", sample_id = "patient_id",
 image_id = "image_number", condition = "patient_stage")

dat$condition <- gsub("-","_", dat$patient_stage) %>% 
  as.factor() %>%
  relevel("Non_diabetic")

# drop rows with NA
dat <- dat |> drop_na()
# calculate the fPCA
pca <- functionalPCA(dat = dat, r = metricRes$r |> unique(), pve = 0.995)
evalues <- pca$evalues
efunctions <- pca$efunctions
# plot the mean curve and the two first eigenfunctions
p_mu <- ggplot(data.frame(r = unique(metricRes$r), mu = pca$mu), 
               aes(x = r, y = mu)) +
    geom_line() +
    theme_light() +
    xlab("r [µm]")

p_efunction1 <- ggplot(data.frame(r = unique(metricRes$r), 
                                  phi1 = pca$efunctions[,1]), 
                       aes(x = r, y = phi1)) +
    geom_line() +
    theme_light() +
    ylim(-0.3,0.3) +
    xlab("r [µm]")

p_efunction2 <- ggplot(data.frame(r = unique(metricRes$r),
                                  phi2 = pca$efunctions[,2]),
                       aes(x = r, y = phi2)) +
    geom_line() +
    theme_light() +
    ylim(-0.3,0.3) +
    xlab("r [µm]")

wrap_plots(list(p_mu, p_efunction1, p_efunction2), ncol = 3)
# plot the biplot of the first two PCs
plotFpca(dat = dat, res = pca, colourby = "condition")
# print the eigenvalues
evalues

## ----funcGamG, fig.height=10, warning = FALSE---------------------------------
library('refund')
# create a design matrix
mm <- model.matrix(~condition, data = dat)
mm %>% head()

r <- metricRes$r |> unique()
# fit the model
mdl <- functionalGam(
    data = dat, x = r,
    designmat = mm, weights = dat$npoints,
    formula = formula(Y ~ 1 + conditionLong_duration +
                          conditionOnset + s(patient_id, bs = "re")),
    family = mgcv::scat(link = "log"),
    algorithm = "bam"
)
summary(mdl)

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

## ----intercept, warning = FALSE, eval = TRUE----------------------------------
# look at the smooth random intercepts per patient
data <- coef(mdl)$smterms$`character(0)`$coef

data <- data %>% left_join(dat %>% 
                             select(patient_id, condition) %>% unique)

p <- ggplot(data, aes(x = yindex.vec, y = value, colour = condition, group = patient_id)) +
  geom_line() +
  theme_light() + 
  geom_smooth(aes(group = 1), col = 'black') +
  xlab("r [µm]")

p

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

