## ----smoppiximage, echo=FALSE-------------------------------------------------
htmltools::img(
    src = knitr::image_uri("../inst/Smoppix.jpg"),
    style = "position:absolute; top:0; right:0; padding:10px;",
    alt = "logo", width = "200px", heigth = "200px"
)

## ----setup, include=FALSE, echo=FALSE-----------------------------------------
knitr::opts_chunk$set(fig.heigh = 7, fig.width = 7)

## ----install, eval = FALSE----------------------------------------------------
# library(Biocmanager)
# install("smoppix")

## ----installAndLoadGitHub, eval = FALSE---------------------------------------
# library(devtools)
# install_github("sthawinke/smoppix")

## ----load---------------------------------------------------------------------
library(smoppix)

## ----biocparallel-------------------------------------------------------------
library(BiocParallel)

## ----nCores-------------------------------------------------------------------
nCores <- 2 # The number of CPUs

## ----nonwindows---------------------------------------------------------------
param <- if (.Platform$OS.type == "unix") {
    # On unix-based systems (linux and macOS), use MulticoreParam
    MulticoreParam(nCores)
} else {
    # On windows, build vignette without multithreading
    SerialParam()
}
register(param)

## ----multiThreadWindows, eval = FALSE-----------------------------------------
# register(SnowParam(workers = nCores, type = "SOCK")) # or doParallel

## ----serialparam, eval = FALSE------------------------------------------------
# register(SerialParam())

## ----loadYang-----------------------------------------------------------------
data(Yang)

## ----headYang-----------------------------------------------------------------
head(Yang)

## ----hypYangs-----------------------------------------------------------------
hypYang <- buildHyperFrame(Yang,
    coordVars = c("x", "y"),
    imageVars = c("day", "root", "section")
)

## ----headHypYang--------------------------------------------------------------
head(hypYang)

## ----zoomin-------------------------------------------------------------------
hypYang$ppp[[1]]

## ----plotExplore, fig.height = 6.5--------------------------------------------
plotExplore(hypYang, Cex = 2)

## ----numFeats-----------------------------------------------------------------
numFeats <- 10 # Limit number of genes in the analysis for computational reasons

## ----estpisyang---------------------------------------------------------------
nnObj <- estPis(hypYang,
    pis = c("nn", "nnPair"), null = "background", verbose = FALSE,
    features = c("SmVND2", "SmPINR", getFeatures(hypYang)[seq_len(numFeats)])
)

## ----adWf---------------------------------------------------------------------
nnObj <- addWeightFunction(nnObj, designVars = c("day", "root"))

## ----plotwfnn-----------------------------------------------------------------
plotWf(nnObj, pi = "nn")

## ----builddf------------------------------------------------------------------
dfUniNN <- buildDataFrame(nnObj, gene = "SmAHK4e", pi = "nn")

## ----lmerManual---------------------------------------------------------------
library(lmerTest, quietly = TRUE)
# Make sure to nest root within day, without including a random effect for day
dfUniNN$dayRoot <- paste0(dfUniNN$day, dfUniNN$root)
lmeMod <- lmerTest::lmer(pi - 0.5 ~ day + (1 | dayRoot),
    data = dfUniNN, na.action = na.omit,
    weights = weight, contrasts = list("day" = "contr.sum")
)
summary(lmeMod)

## ----fitLmms, include = FALSE-------------------------------------------------
allModsNN <- fitLMMs(nnObj, fixedVars = "day", randomVars = "root")

## ----headgetresults-----------------------------------------------------------
head(getResults(allModsNN, "nn", "Intercept"), 4)
head(getResults(allModsNN, "nn", "day"), 4)

## ----plotwf, fig.height = 5---------------------------------------------------
plotWf(nnObj, pi = "nnPair")

## ----builddfbi----------------------------------------------------------------
dfBiNN <- buildDataFrame(nnObj, gene = "SmVND2--SmPINR", pi = "nnPair")

## ----mixmodfit----------------------------------------------------------------
# Make sure to nest root within day, without including a random effect for day
dfBiNN$dayRoot <- paste0(dfBiNN$day, dfBiNN$root)
lmeModNN <- lmerTest::lmer(pi - 0.5 ~ day + (1 | dayRoot),
    data = dfBiNN, na.action = na.omit,
    weights = weight, contrasts = list("day" = "contr.sum")
)
summary(lmeModNN)

## ----visualInvest, fig.height = 6---------------------------------------------
plotExplore(hypYang, features = "SmVND2--SmPINR")

## ----fitlmmPair---------------------------------------------------------------
head(getResults(allModsNN, "nnPair", "Intercept"))

## ----writeXlsxYang, eval = FALSE----------------------------------------------
# writeToXlsx(allModsNN, file = "myfile.xlsx")

## ----readInEngData, warning=FALSE---------------------------------------------
data(Eng)
hypEng <- buildHyperFrame(Eng,
    coordVars = c("x", "y"),
    imageVars = c("fov", "experiment")
)

## ----engNoCells, fig.height = 6-----------------------------------------------
plotExplore(hypEng)

## ----CellEngData--------------------------------------------------------------
# Add cell identifiers
hypEng <- addCell(hypEng, EngRois, verbose = FALSE)

## ----plotExpCell, fig.height = 7----------------------------------------------
plotExplore(hypEng)

## ----engpis-------------------------------------------------------------------
engPis <- estPis(hypEng,
    pis = c("edge", "centroid", "nnCell", "nnPairCell"),
    null = "CSR", verbose = FALSE
)

## ----addWfCell----------------------------------------------------------------
engPis <- addWeightFunction(engPis)

## ----plotWfEng, fig.height = 5------------------------------------------------
plotWf(engPis, "nnPairCell")

## ----FitEngModels-------------------------------------------------------------
allModsNNcell <- fitLMMs(engPis, fixedVars = "experiment")

## ----plotEdgeEng, fig.height = 6.5--------------------------------------------
plotTopResults(hypEng,
    results = allModsNNcell, pi = "edge", what = "far",
    numFeats = 1
)

## ----plotnnCellEng, fig.height = 6.5------------------------------------------
plotTopResults(hypEng, results = allModsNNcell, pi = "nnCell", numFeats = 1)

## ----antiEng, fig.height = 6.5------------------------------------------------
plotTopResults(hypEng,
    results = allModsNNcell, pi = "nnPairCell", numFeats = 1,
    what = "anti"
)

## ----exploreSpatCell, fig.height = 7------------------------------------------
plotExplore(hypEng,
    piEsts = engPis, piColourCell = "edge", features =
        remEdge <- rownames(getResults(allModsNNcell, "edge", "Intercept"))[1]
)

## ----enggrads-----------------------------------------------------------------
engGrads <- estGradients(hypEng[seq_len(2), ], features = getFeatures(hypEng)[seq_len(3)])
pValsGrads <- getPvaluesGradient(engGrads, "cell")

## ----loadTNBC-----------------------------------------------------------------
if (reqFunky <- require(funkycells)) {
    data("TNBC_pheno")
    data("TNBC_meta")
    TNBC_pheno$Age <- TNBC_meta$Age[match(TNBC_pheno$Person, TNBC_meta$Person)]
    TNBC_pheno$Class <- ifelse(TNBC_pheno$Class == "0", "mixed", "compartentalised")
    hypTNBC <- buildHyperFrame(TNBC_pheno,
        coordVars = c("cellx", "celly"),
        imageVars = c("Person", "Class", "Age"), featureName = "Phenotype"
    )
    hypTNBC <- hypTNBC[order(hypTNBC$Class), ] # Order by class for plots
}

## ----plotTNBC, fig.cap = "Explorative plot of triple negative breast cancer dataset.\\parencite{Keren2018}", fig.height = 8----
if (reqFunky) {
    plotExplore(hypTNBC,
        Cex.main = 0.7, features = c("Tumor", "Macrophage", "CD8T"),
        CexLegend = 0.85
    )
}

## ----analyseTNBC--------------------------------------------------------------
if (reqFunky) {
    # Select feature through hypothesis tests
    pisTNBC <- estPis(hypTNBC, pis = c("nn", "nnPair"), null = "background", verbose = FALSE)
    pisTNBC <- addWeightFunction(pisTNBC)
    TNBClmms <- fitLMMs(pisTNBC, fixedVars = c("Age", "Class"), verbose = TRUE)
    # Extract PIs of significant features
    nnRes <- getResults(TNBClmms, "nn", "Class")
    sigLevel <- 0.05
    nnSig <- rownames(nnRes)[nnRes[, "pAdj"] < sigLevel]
    nnPairRes <- getResults(TNBClmms, "nnPair", "Class")
    nnPairSig <- rownames(nnPairRes)[nnPairRes[, "pAdj"] < sigLevel & !is.na(nnPairRes[, "pAdj"])]
    nnPis <- as.matrix(vapply(pisTNBC$hypFrame$pimRes, FUN.VALUE = double(nnn <- (length(nnSig) + length(nnPairSig))), function(x) {
        c(
            if (length(nnSig)) x$pointDists$nn[nnSig],
            if (length(nnPairSig)) {
                vapply(nnPairSig, FUN.VALUE = 1, function(y) {
                    getGp(gp = y, x$pointDists$nnPair, notFoundReturn = NA)
                })
            }
        )
    }))
    if (nnn != 1) {
        nnPis <- t(nnPis)
    }
    colnames(nnPis) <- c(nnSig, nnPairSig)
    nnPis[is.na(nnPis)] <- 0.5 # Replace missing values by 0.5
    if (require(glmnet, quietly = TRUE)) {
        set.seed(20062024)
        lassoModel <- cv.glmnet(factor(hypTNBC$Class), x = cbind("Age" = hypTNBC$Age, nnPis), family = "binomial", alpha = 1, nfolds = 10, type.measure = "class")
        lassoModel # In-sample misclassification error is quite low
    }
}

## ----schematic, echo=FALSE----------------------------------------------------
htmltools::img(
    src = knitr::image_uri("../inst/tikzPic.jpg"),
    alt = "scheme", width = "800px", heigth = "800px"
)

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

