params <-
list(test = FALSE)

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  cache = TRUE,
  message = FALSE,
  warning = FALSE
)
library(BiocStyle)

## ----warning = FALSE, message = FALSE-----------------------------------------
# Loading required packages
library(Statial)
library(spicyR)
library(ClassifyR)
library(lisaClust)
library(dplyr)
library(SingleCellExperiment)
library(ggplot2)
library(ggsurvfit)
library(survival)
library(tibble)

theme_set(theme_classic())
nCores <- 1

## ----eval = FALSE-------------------------------------------------------------
#  # Install the package from Bioconductor
#  if (!requireNamespace("BiocManager", quietly = TRUE)) {
#    install.packages("BiocManager")
#  }
#  BiocManager::install("Statial")

## -----------------------------------------------------------------------------
# Load head and neck data
data("kerenSCE")

kerenSCE

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

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

bcells <- c("B")
tcells <- c("CD3_Cell", "CD4_Cell", "CD8_Cell", "Tregs")
myeloid <- c("Dc/Mono", "DC", "Mono/Neu", "Macrophages", "Neutrophils")

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

tissue <- c(endothelial, mesenchymal)
immune <- c(bcells, tcells, myeloid, "NK", "other immune") # NK = Natural Killer cells

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

## -----------------------------------------------------------------------------
#Select image 6 from the kerenSCE dataset
kerenImage6 = kerenSCE[, kerenSCE$imageID =="6"]

#Select for all cells that express higher than baseline level of p53
p53Pos = assay(kerenImage6)["p53",]  |> 
  as.numeric() > -0.300460
kerenImage6$cellType[p53Pos & kerenImage6$cellType %in% c("Keratin_Tumour")] <- "p53+Tumour"

#Group all immune cells under the name "Immune"
kerenImage6$cellType[kerenImage6$cellType %in% immune] <- "Immune"

kerenImage6 |>
  colData() %>%
  as.data.frame() %>%
  filter(cellType %in% c("Keratin_Tumour", "Immune", "p53+Tumour")) %>%
  arrange(cellType) %>%
  ggplot(aes(x = x, y = y, color = cellType)) +
  geom_point(size = 1) +
  scale_colour_manual(values = c("#505050", "#D6D6D6", "#64BC46"))


## -----------------------------------------------------------------------------

#Select for all cells that express higher than baseline level of p53
kerenSCE$cellTypeNew <- kerenSCE$cellType
p53Pos = assay(kerenSCE)["p53",]  |> 
  as.numeric() > -0.300460
kerenSCE$cellTypeNew[p53Pos & kerenSCE$cellType %in% c("Keratin_Tumour")] <- "p53+Tumour"

#Group all immune cells under the name "Immune"
kerenSCE$cellTypeNew[kerenSCE$cellType %in% immune] <- "Immune"


curves <- kontextCurve(
  cells = kerenSCE,
  from = "p53+Tumour",
  to = "Immune",
  parent = c("p53+Tumour", "Keratin_Tumour"),
  rs = seq(10, 510, 100),
  image = "6",
  cellType = "cellTypeNew",
  cores = nCores
)

kontextPlot(curves)

## -----------------------------------------------------------------------------
p53_Kontextual <- Kontextual(
  cells = kerenSCE,
  r = 50,
  from = "p53+Tumour",
  to = "Immune",
  parent = c("p53", "Keratin_Tumour"),
  cellType = "cellTypeNew"
)

p53_Kontextual


## -----------------------------------------------------------------------------
# Get all relationships between cell types and their parents
parentDf <- parentCombinations(
  all = all,
  tumour,
  bcells,
  tcells,
  myeloid,
  endothelial,
  mesenchymal,
  tissue,
  immune
)


## ----eval=FALSE---------------------------------------------------------------
#  # Running Kontextual on all relationships across all images.
#  kerenKontextual <- Kontextual(
#    cells = kerenSCE,
#    parentDf = parentDf,
#    r = 50,
#    cores = nCores
#  )
#  

## -----------------------------------------------------------------------------
data("kerenKontextual")
head(kerenKontextual, 10)

## -----------------------------------------------------------------------------

kerenSCE <- getDistances(kerenSCE,
                    maxDist = 200,
                    nCores = 1)

kerenSCE <- getAbundances(kerenSCE,
                     r = 200,
                     nCores = 1)


## -----------------------------------------------------------------------------
stateChanges <- calcStateChanges(
  cells = kerenSCE,
  type = "distances",
  image = "6",
  from = "Keratin_Tumour",
  to = "Macrophages",
  marker = "p53",
  nCores = 1)

stateChanges

## -----------------------------------------------------------------------------
p <- plotStateChanges(
  cells = kerenSCE,
  type = "distances",
  image = "6",
  from = "Keratin_Tumour",
  to = "Macrophages",
  marker = "p53",
  size = 1,
  shape = 19,
  interactive = FALSE,
  plotModelFit = FALSE,
  method = "lm")

p


## -----------------------------------------------------------------------------
stateChanges <- calcStateChanges(
  cells = kerenSCE,
  type = "distances",
  nCores = 1,
  minCells = 100)

stateChanges |> head(n = 10)

stateChanges |> 
  filter(imageID == 6) |>
  head(n = 10)

## -----------------------------------------------------------------------------

p <- plotStateChanges(
  cells = kerenSCE,
  type = "distances",
  image = "6",
  from = "Keratin_Tumour",
  to = "Macrophages",
  marker = "HLA_Class_1",
  size = 1,
  shape = 19,
  interactive = FALSE,
  plotModelFit = FALSE,
  method = "lm")

p


## -----------------------------------------------------------------------------

p <- plotStateChanges(
  cells = kerenSCE,
  type = "distances",
  image = "35",
  from = "CD4_Cell",
  to = "B",
  marker = "CD20",
  size = 1,
  shape = 19,
  interactive = FALSE,
  plotModelFit = FALSE,
  method = "lm")

p


## -----------------------------------------------------------------------------
kerenSCE <- calcContamination(kerenSCE)

stateChangesCorrected <- calcStateChanges(
  cells = kerenSCE,
  type = "distances",
  nCores = 1,
  minCells = 100,
  contamination = TRUE)

stateChangesCorrected |> head(n = 20)


## -----------------------------------------------------------------------------
cellTypeMarkers <- c("CD3", "CD4", "CD8", "CD56", "CD11c", "CD68", "CD45", "CD20")

values = c("blue", "red")
names(values) <- c("None", "Corrected")

df <- rbind(data.frame(TP =cumsum(stateChanges$marker %in% cellTypeMarkers), FP = cumsum(!stateChanges$marker %in% cellTypeMarkers), type = "None"),
            data.frame(TP =cumsum(stateChangesCorrected$marker %in% cellTypeMarkers), FP = cumsum(!stateChangesCorrected$marker %in% cellTypeMarkers), type = "Corrected"))

ggplot(df, aes(x = TP, y = FP, colour = type)) + geom_line()+ labs(y = "Cell state marker", x = "Cell type marker") + scale_colour_manual(values = values)

ggplot(df, aes(x = TP, y = FP, colour = type)) + geom_line()+ xlim(0,100) + ylim(0,1000)+ labs(y = "Cell state marker", x = "Cell type marker") + scale_colour_manual(values = values)
            


## -----------------------------------------------------------------------------
# Helper function to run coxPh models

coxTests = function(measurementMat, Surv) {
  
  result = apply(measurementMat, 2, function(measurementCol) {
    fit = coxph(Surv ~ measurementCol)
    summary(fit)$coefficients[1, c(1, 3, 5)]
  })
  
  result = result |> 
    t() |> 
    data.frame() |> 
    rownames_to_column() 
  
  colnames(result) = c("imageID", "coef", "se.coef", "p.value")
  result$`adj.p.value` = p.adjust(result$`p.value`, "fdr")
  
  result = result |> arrange(p.value)
  
  result[,sapply(result,is.numeric)] <- signif(result[,sapply(result,is.numeric)],2)
  
  return(result)
}


## -----------------------------------------------------------------------------

# Extracting survival data to create survival object
survData = kerenSCE |>
    colData() |> 
    data.frame() |> 
    select(imageID, Survival_days_capped, Censored) |> 
    unique()

# Creating survival vector
kerenSurv = Surv(survData$Survival_days_capped, survData$Censored)
names(kerenSurv) = survData$imageID


# Converting Kontextual result into data matrix
kontextMat = prepMatrix(kerenKontextual)

# Ensuring rownames of kontextMat match up with rownames of the survival vector 
kontextMat = kontextMat[names(kerenSurv), ]


# Running survival analysis
survivalResults = coxTests(kontextMat, kerenSurv)

head(survivalResults)


## -----------------------------------------------------------------------------
# Selecting most significant relationship
survRelationship = kontextMat[["Mesenchymal__other immune__tissue"]]
survRelationship = ifelse(survRelationship > median(survRelationship), "Localised", "Dispersed")
    
# Plotting Kaplan-Meier curve
survfit2(kerenSurv ~ survRelationship) |>
    ggsurvfit() +
    add_pvalue() +
    ggtitle("Mesenchymal__other immune__tissue")

## -----------------------------------------------------------------------------
#Preparing features for Statial
stateMat <- prepMatrix(stateChanges)

# Ensuring rownames of stateMat match up with rownames of the survival vector 
stateMat = stateMat[names(kerenSurv), ]

survivalResults = coxTests(stateMat, kerenSurv)

head(survivalResults)


## -----------------------------------------------------------------------------
survRelationship = stateMat[["Keratin_Tumour__CD4_Cell__Keratin6"]]
survRelationship = ifelse(survRelationship > median(survRelationship), "Dispersed", "Localised")
    
# Plotting Kaplan-Meier curve
survfit2(kerenSurv ~ survRelationship) |>
    ggsurvfit() +
    add_pvalue() +
    ggtitle("Keratin_Tumour__CD4_Cell__Keratin6")


## -----------------------------------------------------------------------------
set.seed(51773)

# Preparing features for lisaClust
kerenSCE <- lisaClust::lisaClust(kerenSCE, k = 5)


## ----fig.height=5, fig.width=6.5----------------------------------------------
# Use hatching to visualise regions and cell types.
lisaClust::hatchingPlot(kerenSCE,
  useImages = "5",
  line.spacing = 41, # spacing of lines
  nbp = 100 # smoothness of lines
) 


## ----lisaClust----------------------------------------------------------------

cellTypeRegionMeans <- getMarkerMeans(kerenSCE,
                              imageID = "imageID",
                              cellType = "cellType",
                              region = "region")

survivalResults = coxTests(cellTypeRegionMeans[names(kerenSurv),], kerenSurv)

head(survivalResults)


## -----------------------------------------------------------------------------

# Calculate cell type and region proportions
cellTypeProp <- getProp(kerenSCE, 
                       feature = "cellType",
                       imageID = "imageID")
regionProp <- getProp(kerenSCE, 
                       feature = "region",
                       imageID = "imageID")

# Combine all the features into a list for classification 
featureList <- list(states = stateMat, 
                     kontextual = kontextMat,
                     regionMarkerMeans = cellTypeRegionMeans,
                     cellTypeProp = cellTypeProp,
                     regionProp = regionProp)

# Ensure the rownames of the features match the order of the survival vector
featureList <- lapply(featureList, function(x)x[names(kerenSurv),])


set.seed(51773)

kerenCV = crossValidate(
  measurements = featureList,
  outcome = kerenSurv,
  classifier = "CoxPH",
  selectionMethod  = "CoxPH",
  nFolds = 3,
  nFeatures = 10,
  nRepeats = 20,
  nCores = 1
  )


## -----------------------------------------------------------------------------
# Calculate AUC for each cross-validation repeat and plot.
performancePlot(kerenCV,
  characteristicsList = list(x = "Assay Name")
  ) + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

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