## ----echo=FALSE---------------------------------------------------------------
htmltools::img(src = knitr::image_uri("../inst/Smoppix.jpg"), 
               alt = 'logo', 
               style = 'position:absolute; top:0; right:0; padding:10px;',
               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)

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

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

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

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

## -----------------------------------------------------------------------------
library(lmerTest, quietly = TRUE)
lmeMod <- lmerTest::lmer(pi - 0.5 ~ day + (1 | root),
    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")

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

## -----------------------------------------------------------------------------
lmeModNN <- lmerTest::lmer(pi - 0.5 ~ day + (1 | root),
    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", addMoransI = TRUE)

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

## ----testMoransI--------------------------------------------------------------
getResults(allModsNNcell, "edge", "Intercept", moransI = TRUE)[remEdge, ]

## ----writeXlsxEng, eval = FALSE-----------------------------------------------
# writeToXlsx(allModsNNcell, file = "mysecondfile.xlsx")

## ----enggrads-----------------------------------------------------------------
engGrads <- estGradients(hypEng[seq_len(2), ], features = feat <- 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 = t(vapply(pisTNBC$hypFrame$pimRes, FUN.VALUE = double(length(nnSig)+ length(nnPairSig)), function(x){
        c(if(length(length(nnSig))) x$pointDists$nn[nnSig],
          if(length(length(nnPairSig))) vapply(nnPairSig, FUN.VALUE = 1, function(y) {
              tmp = getGp(gp = y, x$pointDists$nnPair)
              if(is.null(tmp)) NA else tmp
          }))
    }))
    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
    }
}

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