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

set.seed(1234)

## ----loading, warning = FALSE, message = FALSE--------------------------------
# load the dataset as SpatialExperiment object
spe <- imcdatasets::Damond_2019_Pancreas("spe", full_dataset = TRUE)
# subset the example for computational reasons
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)

## ----Lfunction, warning = 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-----------------------------------------------
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")

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

# create meta info of the IDs
splitData <- dat$ID %>%
  str_replace("-","_") %>%
  str_split_fixed("x", 3) %>% 
  data.frame(stringsAsFactors = TRUE) %>%
  setNames(c("condition", "patient_id", "imageId")) %>%
  mutate(condition = relevel(condition,"Non_diabetic"))
dat <- cbind(dat, splitData)

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

## ----funcGamG, fig.height=10, warning = FALSE---------------------------------
library('refund')
# create a design matrix
mm <- model.matrix(~condition, data = dat)
colnames(mm)[1] <- "Intercept"
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 + pcre(id=patient_id, 
                                                efunctions=efunctions, 
                                                evalues=evalues, yind=r)),
    family = quasi(link = "identity", variance = "mu^2"),
    algorithm = "gam"
)
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()
resid(mdl) |> cov() |> filled.contour()
try(refund::pffr.check(mdl))

## ----intercept, warning = FALSE, eval = TRUE----------------------------------
# look at the smooth random intercepts per patient
#data <- coef(mdl)$smterms$`s(patient_id)`$coef
data <- coef(mdl)$smterms$`pcre(patient_id,efunctions,evalues,r)`$coef %>%
  dplyr::rename(patient_id = "patient_id.vec")

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

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

p

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