## ----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---------------------------------------------------------------
if(.Platform$OS.type == "unix"){
    #On unix-based systems, use MulticoreParam
    register(MulticoreParam(nCores))
} else {
    #On windows, use makeCluster
    library(doParallel)
    Clus = makeCluster(nCores)
    registerDoParallel(Clus)
    register(DoparParam(), default = TRUE)
}

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

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

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

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

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

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

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

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

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

## -----------------------------------------------------------------------------
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 <- estPis(hypEng, pis = c("edge", "centroid", "nnCell", "nnPairCell"),
    null = "CSR", verbose = FALSE)

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

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

## ----a, 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
    }
}

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

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