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

## ----setup, message=FALSE, warning=FALSE--------------------------------------
library("multipointR")
library("SpatialExperiment")
library("SpatialFeatureExperiment")
library("dplyr")
library("ggplot2")
library("spatstat.model")
library("patchwork")

## ----dataloading, warning=FALSE, message=FALSE--------------------------------
spe <- SpatialDatasets::spe_Keren_2018()
spe

## ----preprocessDataset, fig.height=50, fig.width=15---------------------------
# Code source: https://www.bioconductor.org/packages/release/bioc/vignettes/
# Statial/inst/doc/Statial.html
# #kontextual-identifying-discrete-changes-in-cell-state

# Examine all cell types in image
unique(spe$cellType)

# Set up cell populations
tumour <- c("Keratin_Tumour", "Tumour")

bcells <- c("B_cell")
tcells <- c("dn_T_CD3", "CD4_T_cell", "CD8_T_cell", "Tregs")
myeloid <- c("DC_or_Mono", "DC", "Mono_or_Neu", "Macrophages", "Neutrophils")

endothelial <- c("Endothelial")
mesenchymal <- c("Mesenchymal")

tissue <- c(endothelial, mesenchymal)
immune <- c(bcells, tcells, myeloid, "NK", "other immune")

all <- c(tumour, tissue, immune, "Unidentified")

# Lets define a new cell type vector
spe$cellTypeNew <- spe$cellType

# Select for all cells that express higher than baseline level of p53
p53Pos <- assay(spe)["p53", ] > -0.300460

# Find p53+ tumour cells
spe$cellTypeNew[spe$cellType %in% tumour] <- "Tumour"
spe$cellTypeNew[p53Pos & spe$cellType %in% tumour] <- "p53_Tumour"

# Group all immune cells under the name "Immune"

spe$cellTypeNew[spe$cellType %in% immune] <- "Immune"

spe$cellTypeNew <- as.factor(spe$cellTypeNew)
speSub <- subset(spe, , cellTypeNew %in% c("Immune", "Tumour", "p53_Tumour"))
speSub$cellTypeNew <- factor(speSub$cellTypeNew,
    levels = c("Immune", "Tumour", "p53_Tumour")
)

## ----singleImagePlot, echo=TRUE-----------------------------------------------
# Code source: https://www.bioconductor.org/packages/release/bioc/vignettes/
# Statial/inst/doc/Statial.html
# #kontextual-identifying-discrete-changes-in-cell-state

# Plot image 6

df <- spe |>
    colData() |>
    cbind(spatialCoords(spe)) |>
    as.data.frame() |>
    dplyr::filter(imageID == "6") |>
    dplyr::filter(cellTypeNew %in%
        c("Immune", "Tumour", "p53_Tumour"))

df$cellTypeNew <- factor(df$cellTypeNew,
    levels = c("Immune", "Tumour", "p53_Tumour")
)
p1 <- df |>
    arrange(cellTypeNew) |>
    ggplot(aes(x = x, y = y, color = cellTypeNew)) +
    geom_point(size = 1) +
    scale_colour_manual(
        values = c("#505050", "#D6D6D6", "#64BC46"),
        labels = c("Immune", "Tumour", "p53+ Tumour")
    ) +
    guides(colour = guide_legend(
        title = "Cell types",
        override.aes = list(size = 5)
    )) +
    coord_equal() +
    theme_light()

p1

## ----Ppm0---------------------------------------------------------------------
speSub <- subset(spe, , imageID == "6")

m0 <- fitModel(
    spe = speSub,
    marks = "cellTypeNew",
    interaction = "Hardcore",
    formula = as.formula("p53_Tumour ~ 1")
)

m0

## ----Ppm0Plot, fig.width=8, fig.height=8--------------------------------------
plot(m0)

## ----Ppm1---------------------------------------------------------------------
m1 <- fitModel(
    spe = speSub,
    marks = "cellTypeNew",
    formula = as.formula("p53_Tumour ~ s(x,y)"),
    interaction = "Fiksel",
    use.gam = TRUE
)
m1

## ----Ppm1Anova----------------------------------------------------------------
anova(m1, test = "LRT")

## ----Ppm1Plot, fig.width=8, fig.height=8--------------------------------------
plot(m1)

## ----Ppm2---------------------------------------------------------------------
m2 <- fitModel(
    spe = speSub,
    marks = "cellTypeNew",
    formula = as.formula("p53_Tumour ~ s(x,y) + distfun(Immune)"),
    interaction = "Fiksel",
    use.gam = TRUE
)

m2

## ----Ppm2Anova----------------------------------------------------------------
anova(m2, test = "LRT")

## ----Ppm2Plot, fig.width=8, fig.height=8--------------------------------------
plot(m2)

## ----Ppm2diagnose, fig.width=8, fig.height=8----------------------------------
diagnose.ppm(m2)

## ----Ppm2Qq, fig.width=8, fig.height=8, message = FALSE, results = "hide"-----
p <- qqplot.ppm(m2, nsim = 50)

## ----Ppm2QqPlot---------------------------------------------------------------
p

## ----sostaSegmentation, fig.width=8, fig.height=8-----------------------------
segmentedTumour <- sosta::reconstructShapeDensityImage(
    speSub,
    marks = "cellTypeNew",
    markSelect = c("Tumour"),
    thres = 1.1e-4
)

p1 + geom_sf(
    data = segmentedTumour, inherit.aes = FALSE,
    fill = NA, color = "red", linewidth = 1
)

## ----speToSfe-----------------------------------------------------------------
sfeSub <- toSpatialFeatureExperiment(speSub)
annotGeometry(sfeSub, "tumour_mask") <- segmentedTumour

## ----Ppm3---------------------------------------------------------------------
m3 <- fitModel(
    spe = sfeSub,
    marks = "cellTypeNew",
    formula = as.formula("p53_Tumour ~ s(x,y) + tumour_mask + distfun(Immune)"),
    interaction = "Fiksel",
    use.gam = TRUE
)

m3

## ----Ppm3Anova----------------------------------------------------------------
anova(m3, test = "LRT")

## ----Ppm3Plot, fig.width=8, fig.height=8--------------------------------------
plot(m3) + geom_sf(
    data = segmentedTumour, inherit.aes = FALSE,
    fill = NA, color = "darkred", linewidth = 0.6
)

## ----Ppm3Plotdiagnose, fig.width=8, fig.height=8------------------------------
diagnose.ppm(m3)

## ----mPpmLs-------------------------------------------------------------------
mdlLs <- fitModelAcrossImages(
    spe = spe,
    imageId = "imageID",
    marks = "cellTypeNew",
    sharedModel = FALSE,
    interaction = "Hardcore",
    formula = as.formula("p53_Tumour ~ s(x,y) + distfun(Immune)"),
    use.gam = TRUE,
    threshold = 10
)

## ----mPpmToDf, warning=FALSE, message = FALSE, error = FALSE, results = "hide"----
mdlDf <- mdlToDf(
    mdlLs = mdlLs,
    imageCovariates = c(
        "imageID",
        "tumour_type"
    )
)

## ----coefficientsmPpm,  message=FALSE, warning=FALSE, error=FALSE-------------
mdlDfSub <- mdlDf %>% filter(covariate %in% c("distfun.Immune."))

ggplot(mdlDfSub, aes(x = covariate, y = Estimate, label = imageID)) +
    geom_boxplot(outlier.shape = NA, alpha = 0.3) +
    geom_jitter(aes(color = log10(S.E.)),
        position = position_jitter(seed = 123)
    ) +
    geom_text(aes(color = log10(S.E.)),
        hjust = 0,
        vjust = 0, position = position_jitter(seed = 123)
    ) +
    theme_light() +
    facet_wrap(~tumour_type)

## ----model1213Plot, fig.width=8, fig.height=8---------------------------------
# model 12
plot(mdlLs[[12]])
diagnose.ppm(mdlLs[[12]])
# model 13
plot(mdlLs[[13]])
diagnose.ppm(mdlLs[[13]])

## ----SecondStageLm------------------------------------------------------------
mdl <- lm(Estimate ~ tumour_type,
    data = mdlDf,
    weights = 1 / (mdlDf$S.E.),
    subset = covariate == "distfun.Immune."
)

print(summary(mdl))

## ----mPpmShared---------------------------------------------------------------
mdl <- fitModelAcrossImages(
    spe = spe,
    imageId = "imageID",
    imageLs = seq(10, 15),
    marks = "cellTypeNew",
    interaction = "Fiksel",
    formula = as.formula("p53_Tumour ~  log(lambda) +
        tumour_type:distfun(Immune) + (1|DONOR_NO)"),
    threshold = 10
)
mdl |> summary()

## ----mPpmSharedResiduals------------------------------------------------------
res <- residuals(mdl, type = "raw")
resInt <- sapply(res, integral.msr)
print(resInt)

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