mxCODEX.Rmd

---
title: "Run modules"
output: html_document
date: '2022-08-16'
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Examples to run functions and modules

```{r}
Sys.setenv(KMP_DUPLICATE_LIB_OK = "TRUE")
# set test variables
# pID <- "pEdOoZ"
pID <- "Fu1zEC"
# versionID <- 2
versionID <- 1
# projectsDir <- "/Volumes/Analysis_SSD/Dominik/cecelia/projects/"
projectsDir <- "/Volumes/USER_data/Dominik/CECELIA_BACKUP/"
hpcDir <- "/data/scratch/projects/punim1124/cecelia/USERS/schienstockd/"
```

```{r}
library(flowWorkspace)

devtools::load_all("../")
cciaUse("~/cecelia/dev")

library(ggplot2)
library(tidyverse)
```

```{r}
# HPC config
HPC_CONF <-  list(
  hpc = list(
    conf = list(
      email = "schienstockd@student.unimelb.edu.au",
      emailOnBegin = FALSE,
      emailOnEnd = FALSE,
      emailOnFail = TRUE,
      numNodes = 1,
      numTasks = 1,
      numCPUperTask = 1,
      numGPUperTask = 1,
      memory = 50,
      # walltime = "00-01:00:00",
      walltime = "00-00:05:00",
      # projectPartitions = "cascade,sapphire",
      projectPartitions = "gpu-a100,gpu-a100-short",
      # projectID = "punim1124",
      projectID = "punim1031",
      # useGPU = FALSE
      useGPU = TRUE
    )
  ),
  utils = list(
    ssh = list(
      username = "schienstockd",
      address = "spartan.hpc.unimelb.edu.au",
      keyfile = "/Users/Dominik/.ssh/spartan_rsa"
    ),
    python = list(
      condaEnv = cciaConf()$python$conda$source$env
    )
  )
)
```

```{r}
# anaDir <- "/Volumes/USER_data/Dominik/Experiments/Phillips2021_CODEX/CECELIA"
# anaDir <- "/Volumes/USER_data/Dominik/Experiments/CODEX_THEO/ANALYSIS"
anaDir <- "/Volumes/MDHS-Research/5300/5260/Lab-Mueller/Dominik Schienstock/Notebook/Experiments/HuBMAP/ANALYSIS"
# anaDir <- "/Volumes/MDHS-Research/5300/5260/Lab-Mueller/Dominik Schienstock/Notebook/Experiments/CODEX_HCC/ANALYSIS"
```

```{r}
# baseChannels <- c(
#   # 'CD11b','CD16','CD164','EGFR','p53',
#   'MUC-1','Vimentin','Cytokeratin',
#   'PD-L1','CD15','CD30','CD2','GranzymeB',
#   'CD5','MMP-9','CD4','LAG-3','CD25','CD56',
#   'CD20','PD-1','CD11c','CD162',
#   'CD194','IDO-1','VISTA',
#   'HLA-DR','ICOS','BCL-2','CD3','CD69',
#   'CD8','CD7','CD45RA','CD45','CD1a',
#   'CD57','B-catenin','CD45RO','CD71','CD34',
#   'CD68','CD38','Collagen-IV','CD31',
#   'Podoplanin','CD138','CD163',
#   'Mast-cell-tryptase','MMP-12'
# )
# nucChannels <- c(
#   # 'Hoechst', 'DRAQ5', 'FOXP3', 'GATA3', 'T-bet', 'Ki-67'
#   'FOXP3', 'GATA3', 'T-bet', 'Ki-67'
# )

# get cluster channels
baseChannels <- list(
  init = c(
    'CD2','CD3','CD4','CD5','CD7','CD25','CD30','CD69','MMP-12',
    'MMP-9','CD45','CD45RA','CD45RO','GranzymeB','ICOS','LAG-3',
    'PD-1','PD-L1','CD11b','CD163','CD68','CD16','CD56','CD57',
    'CD138','CD20','CD38','CD15','Mast-cell-tryptase','CD1a',
    'CD11c','BCL-2','B-catenin','CD71','EGFR','HLA-DR','IDO-1',
    'Cytokeratin','MUC-1','CD31','CD34','Podoplanin','Collagen-IV',
    'Vimentin','VISTA','CD164','p53','CD162','CD194','CD8'
  )
  # CD4 = c(
  #   'GranzymeB','ICOS','MMP-9',
  #   'LAG-3','PD-1','PD-L1','VISTA',
  #   'CD2','CD5','CD7','CD25','CD30','CD69',
  #   'CD162','CD164','CD194','MMP-12','p53',
  #   'CD45RO', 'CD45RA'
  # )
)

nucChannels <- list(
  # init = c('Hoechst', 'DRAQ5','FOXP3','GATA3','T-bet','Ki-67')
  init = c('Hoechst', 'DRAQ5')
  # CD4 = c('FOXP3','GATA3','T-bet','Ki-67')
)

objectMeasures <- list(
  init = c('prolate','oblate','nuc_area','area','solidity','eccentricity','extent')
  # CD4 = c('prolate','oblate','nuc_area','area','solidity','eccentricity','extent')
  # 'orientation','asymmetry',
  # 'nc_ratio','aspect_ratio','fill'
  # 'perimeter', 'convex_area','equivalent_diameter',
  # 'feret_diameter_max','major_axis_length', 'minor_axis_length',
  # 'perimeter_to_area'
)

clusterRes <- list(
  init = 0.4
  # CD4 = 0.8
  )
```

```{r}
# check that the channel names are fine
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "hR7gDX", versionID = versionID, initReactivity = FALSE 
)

cciaObj$popPaths("flow")
```

```{r}
# check that the channel names are fine
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "akAcC1", versionID = versionID, initReactivity = FALSE # Set
)

uIDs <- names(cciaObj$cciaObjects())
uIDs <- uIDs[!uIDs %in% c("1I5oli")]

for (x in cciaObj$cciaObjects(uIDs = uIDs)) {
  channelNames <- sapply(
  xml2::xml_attrs(xml2::xml_find_all(
    x$omeXML(), x$omeXMLPath("//Image//Pixels//Channel"))),
  function(x) {x[[2]]})

  x$setImChannelNames(channelNames, valueName = "default")
  x$saveState()
}
```

```{r}
# check that the channel names are fine
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  # pID = pID, uID = "jGcZHp", versionID = versionID, initReactivity = FALSE # Set
  # pID = pID, uID = "akAcC1", versionID = versionID, initReactivity = FALSE # Set
  # pID = pID, uID = "pN5QJi", versionID = versionID, initReactivity = FALSE # Set
  pID = pID, uID = "DPKHD7", versionID = versionID, initReactivity = FALSE # Set
)

uIDs <- names(cciaObj$cciaObjects())
# uIDs <- uIDs[!uIDs %in% c("1I5oli")]

for (x in cciaObj$cciaObjects(uIDs = uIDs)) {
  print(paste(">>", x$getUID()))
  
  if (stringr::str_detect(x$oriFilepath(), "_core[0-9]+"))
    omexmlFile <- stringr::str_replace(x$oriFilepath(), "_core[0-9]+\\.tif", "-omexml.txt")
  else
    omexmlFile <- stringr::str_replace(x$oriFilepath(), ".qptiff", "-omexml.txt")
  # omexmlFile <- stringr::str_replace(x$oriFilepath(), ".ome.tif", ".txt")
  # omexmlFile <- stringr::str_replace(x$oriFilepath(), ".tif", ".txt")
  
  omeXML <- readLines(omexmlFile, warn = FALSE)
  channelNames <- stringr::str_extract(omeXML, "(?<=Name=\").+(?=\" )")
  channelNames <- channelNames[!is.na(channelNames)]
  
  # correct channel names
  omeXML[omeXML == "HLADR"] <- "HLA-DR"
  
  x$setImChannelNames(channelNames, valueName = "default")
  # x$setImChannelNames(omeXML, valueName = "default")
  x$saveState()
}
```

```{r}
# check that the channel names are fine
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "akAcC1", versionID = versionID, initReactivity = FALSE # Set
)

for (x in cciaObj$cciaObjects()) {
  x$setImFilepath("ccidAfCorrected.zarr", "afCorrected")
  x$saveState()
}
```

```{r}
# check that the channel names are fine
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "oUPuBX", versionID = versionID, initReactivity = FALSE # Human HCC
)

uIDs <- exp.info[Include == "Y"]$uID

for (x in cciaObj$cciaObjects(uIDs = uIDs)) {
  imChannelNames <- x$imChannelNames()
  imChannelNames[imChannelNames == "pancyto"] <- "pancytokeratin"
  imChannelNames[imChannelNames == "TCRValpha"] <- "TCRVa"
  imChannelNames[imChannelNames == "FOXP3"] <- "FoxP3"
  imChannelNames[imChannelNames == "CD4 good"] <- "CD4"
  imChannelNames[imChannelNames == "ILR18a"] <- "IL18Ra"
  x$setImChannelNames(imChannelNames)
  x$saveState()
}
```

```{r}
# upload objects
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE, projectsDir = projectsDir)

# init ccia object
cciaObj <- initCciaObject(
  # pID = pID, uID = "jGcZHp", versionID = versionID, initReactivity = FALSE 
  # pID = pID, uID = "oUPuBX", versionID = versionID, initReactivity = FALSE 
  pID = pID, uID = "DPKHD7", versionID = versionID, initReactivity = FALSE
)

# uIDs <- exp.info[Include == "Y"]$uID
# uIDs <- uIDs[!uIDs %in% c("AtIzds")]
uIDs <- names(cciaObj$cciaObjects())

# run task
funParams <- list(
)

# run task
task <- cciaObj$runTask(
  funName = "hpc.uploadCciaObj",
  funParams = funParams,
  envVars = HPC_CONF,
  hpcDir = paste(hpcDir, pID, "ANALYSIS", sep = "/"),
  env = "local",
  runInplace = FALSE,
  taskID = 1
)

task <- cciaObj$runTasks(
  funName = "hpc.uploadCciaObj",
  funParams = funParams,
  envVars = HPC_CONF,
  hpcDir = paste(hpcDir, pID, "ANALYSIS", sep = "/"),
  env = "local",
  runInplace = FALSE,
  taskID = 1,
  uIDs = uIDs,
  mc.cores = 3
)
```


```{r}
# check that the channel names are fine
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "akAcC1", versionID = versionID, initReactivity = FALSE # Set
)

uIDs <- names(cciaObj$cciaObjects())
# uIDs <- uIDs[uIDs %in% c("Rl72En", "lNNoqO", "J7z78G", "CcIXqP", "7iyMr9")]

for (x in cciaObj$cciaObjects(uIDs = uIDs)) {
  print(paste(">>", x$getUID()))
  
  imChannels <- unlist(x$imChannelNames())
  numChannels <- length(imChannels)
  seqChannels <- seq(numChannels) - 1
  names(seqChannels) <- seqChannels
  
  # prepare af combinations
  afCombinations <- lapply(seqChannels, function(y) {
    list(
      # divisionChannels = unname(seqChannels[seqChannels != x]),
      divisionChannels = unname(seqChannels),
      channelPercentile = 0,
      correctionPercentile = 0,
      correctionMode = 'divide',
      summaryMode = 'percentile',
      summaryPercentile = 75,
      correctionMin = 0.1,
      correctionMax = 99.9,
      generateInverse = FALSE,
      medianFilter = 0,
      topHatRadius = 0,
      rollingBallRadius = 0,
      rollingBallPadding = 0
    )
  })
  
  # run task
  funParams <- list(
    valueName = "default",
    ignoreVisibilities = TRUE,
    afCombinations = afCombinations,
    applyGaussian = FALSE,
    applyGaussianToOthers = FALSE
  )
  
  task <- x$runTask(
    funName = "cleanupImages.afCorrect",
    funParams = funParams,
    envVars = HPC_CONF,
    hpcDir = paste(hpcDir, pID, "ANALYSIS", sep = "/"),
    env = "hpc",
    # env = "local",
    runInplace = FALSE,
    taskID = 1
  )
}
```


```{r}
# set cluster property files
# TODO why is this not done.. ?
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  # pID = pID, uID = "jGcZHp", versionID = versionID, initReactivity = FALSE # Set
  pID = pID, uID = "pN5QJi", versionID = versionID, initReactivity = FALSE # Set
)

exp.info <- as.data.table(cciaObj$summary(withSelf = FALSE, fields = c("Attr")))

# uIDs <- exp.info[Include == "Y"]$uID
uIDs <- exp.info[Stain == "A"]$uID

for (x in cciaObj$cciaObjects(uIDs = uIDs[1])) {
  # x$setImAnndataFilepath("default.clust.h5ad", valueName = "default")
  # x$saveState()
}

x$imNeighboursFilepath()
```

```{r}
# run clustering
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "jGcZHp", versionID = versionID, initReactivity = FALSE # THEO MOUSE
  # pID = pID, uID = "akAcC1", versionID = versionID, initReactivity = FALSE # Set
  # pID = pID, uID = "pN5QJi", versionID = versionID, initReactivity = FALSE # Human LN
  # pID = pID, uID = "oUPuBX", versionID = versionID, initReactivity = FALSE # Human HCC
  # pID = pID, uID = "DPKHD7", versionID = versionID, initReactivity = FALSE # THEO Human
)

# get experimental info
exp.info <- as.data.table(cciaObj$summary(
  withSelf = FALSE, fields = c("Attr")
))

uIDs <- names(cciaObj$cciaObjects())
# uIDs <- exp.info[Include == "Y"]$uID
uIDs <- exp.info[Analyse == "Y"]$uID # THEO MOUSE
# uIDs <- exp.info[Stain == "A"]$uID
# uIDs <- exp.info[Include == "Y" & Stain == "B"]$uID
# uIDs <- exp.info[Include == "Y" & Stain == "A"]$uID

# similar LN
# uIDs <- c("1oBdEQ", "Eqiklk", "BnSXJp", "gDCQ7s", "sZBvMm", "TMy4g4", "YbP2O0")

# get common channels across images
channelNames <- lapply(cciaObj$cciaObjects(uIDs = uIDs), function(x) unname(x$imChannelNames()))
clusterChannels <- Reduce(intersect, channelNames)
# clusterChannels <- channelNames
# clusterChannels <- clusterChannels[clusterChannels != "DAPI"]
clusterChannels <- clusterChannels[!clusterChannels %in% c("DAPI", "DAPI-02", "HOECHST")]
# clusterChannels <- c("CD3", "CD20")

# exclude channels
clusterChannels <- clusterChannels[!clusterChannels %in% c(
  # THEO MOUSE
  # These markers did not seem to have worked
  "CD20", "CD11c", "LamB1", "CD8", "CD206/MRC1", "CD45",
  "Tryptase", "NF", "TMEM119", "LamB1", "SOX2"
  # NEXT
  # "Lyve1", "SMActin", "Vimentin", "ECAD", "CollagenIV", "Podoplanin"
  # "ECAD", "CollagenIV", "IgM"
  )]

# run task
funParams <- list(
  # popType = "clust",
  popType = "flow",
  popsToCluster = c(),
  # popsToCluster = c("Zero"),
  savePops = TRUE,
  # keepPops = TRUE,
  keepPops = FALSE,
  # resolution = 0.2, # THEO MOUSE BRAIN, initial
  # resolution = 0.4, # Human LN
  # resolution = 0.6, # Human HCC
  resolution = 0.8, # THEO MOUSE BRAIN, 2nd
  mergeUmap = TRUE,
  normalisePercentile = 99.8,
  normalisePercentileBottom = 0.0,
  normaliseAxis = "channels",
  normaliseToMedian = FALSE, # Human LN and HCC
  # normaliseToMedian = TRUE, # MOUSE THEO
  normaliseIndividually = FALSE,
  # normaliseIndividually = TRUE,
  transformation = "NONE",
  # transformation = "log",
  logBase = 0,
  correctBatch = "NONE",
  # correctBatch = "combat",
  # correctBatch = "harmony",
  valueName = "default",
  # valueName = "detailed",
  uIDs = uIDs,
  clusterChannels = list(
    base = list(channels = clusterChannels)
  ),
  objectMeasures = list(),
  # usePaga = TRUE,
  usePaga = FALSE,
  pagaThreshold = 0.1
)

# go through and cluster individually
i <- 1
for (x in cciaObj$cciaObjects(uIDs = uIDs)) {
  funParams$uIDs <- c(x$getUID())
  
  # get populations to do clustering on
  # THEO MOUSE
  funParams$popsToCluster <- unlist(unname(x$popPaths("flow")))
  
  task <- cciaObj$runTask(
    funName = "clustPopulations.leidenClustering",
    funParams = funParams,
    # envVars = HPC_CONF,
    # hpcDir = paste(hpcDir, pID, "ANALYSIS", sep = "/"),
    # env = "hpc",
    env = "local",
    # runInplace = FALSE,
    runInplace = TRUE,
    taskID = i
  )

  i <- i + 1
}

# task <- cciaObj$runTask(
#   funName = "clustPopulations.leidenClustering",
#   funParams = funParams,
#   envVars = HPC_CONF,
#   hpcDir = paste(hpcDir, pID, "ANALYSIS", sep = "/"),
#   # env = "hpc",
#   env = "local",
#   # runInplace = FALSE,
#   runInplace = TRUE,
#   taskID = 1
# )
```

```{r}
# run clustering
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  # pID = pID, uID = "jGcZHp", versionID = versionID, initReactivity = FALSE # THEO Mouse
  # pID = pID, uID = "akAcC1", versionID = versionID, initReactivity = FALSE # Set
  # pID = pID, uID = "pN5QJi", versionID = versionID, initReactivity = FALSE # Set
  # pID = pID, uID = "oUPuBX", versionID = versionID, initReactivity = FALSE # Set
  pID = pID, uID = "DPKHD7", versionID = versionID, initReactivity = FALSE # THEO Human
)

# get experimental info
exp.info <- as.data.table(cciaObj$summary(
  withSelf = FALSE, fields = c("Attr")
))

uIDs <- names(cciaObj$cciaObjects())
# uIDs <- exp.info[Include == "Y"]$uID
# uIDs <- exp.info[Stain == "A"]$uID

for (x in cciaObj$cciaObjects(uIDs = uIDs)) {
  x$setImAnndataFilepath("default.clust.h5ad")
  # x$valuePartOf("imAnndataFilepath")
  # x$setImNeighboursFilepath("default.sq.h5ad")
  # x$setImRegionsFilepath("default.region.h5ad")
  x$saveState()
}
```

```{r fig_umap_clustered, fig.height=6, fig.width=6}
# run clustering
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "jGcZHp", versionID = versionID, initReactivity = FALSE # THEO MOUSE
  # pID = pID, uID = "oUPuBX", versionID = versionID, initReactivity = FALSE # Human HCC
)

# get experimental info
exp.info <- as.data.table(cciaObj$summary(
  withSelf = FALSE, fields = c("Attr")
))

uIDs <- names(cciaObj$cciaObjects())
uIDs <- exp.info[Include == "Y"]$uID

# run task
funParams <- list(
  valueName = "default",
  popType = "clust",
  uIDs = uIDs,
  prepClusterNum = "min",
  numClusters = 0
)

task <- cciaObj$runTask(
  funName = "clustPopulations.integrateClustering",
  funParams = funParams,
  env = "local",
  runInplace = TRUE,
  taskID = 1
)

# check UMAP
umapDF <- readRDS(file.path(cciaObj$persistentObjectDirectory(), "data", "default.cluster.groups.rds"))
colPal <- randomcoloR::distinctColorPalette(length(unique(umapDF$clusters.group)))

# umapDF <- .mapClustNames(umapDF, clusterMapping, "clusters.group", removeNone = FALSE)

p1 <- ggplot(umapDF, aes(V1, V2)) +
  # geom_point(aes(color = uID), size = 2) +
  geom_point(aes(color = clusters.group), size = 2) +
  # geom_point(aes(color = clusters.name), size = 3) +
  theme_classic() +
  # scale_color_brewer(name = NULL, palette = "Set1") +
  scale_color_manual(values = colPal) +
  # scale_color_manual(name = NULL, values = clusterColours) +
  # viridis::scale_colour_viridis() +
  theme(
    axis.text = element_text(size = 15),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    legend.justification = "right",
    legend.position = "bottom",
    # legend.position = "none",
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
    ) +
  guides(colour = guide_legend(override.aes = list(size=5)))

ggsave(file.path(anaDir, "group_umap.png"), width = 4, height = 5, plot = p1)
```

```{r}
# can you cluster and check if there is a difference in R and N?
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "oUPuBX", versionID = versionID, initReactivity = FALSE # Set
)

# get experimental info
exp.info <- as.data.table(cciaObj$summary(
  withSelf = FALSE, fields = c("Attr")
))

uIDs <- names(cciaObj$cciaObjects())
uIDs <- exp.info[Include == "Y"]$uID

# run task
funParams <- list(
  uIDs = names(cciaObj$cciaObjects()),
  popType = c("clust"),
  pops = c("root"),
  objectMeasures = c(),
  includeCounts = FALSE,
  includeTotalCount = FALSE,
  colName = "default",
  valueName = "default",
  excludeSelfConnections = FALSE,
  useClusters = FALSE,
  useClustersGroup = TRUE,
  expandPops = FALSE,
  numClusters = 6,
  uIDs = uIDs
)

task <- cciaObj$runTask(
  funName = "clustRegions.kmeansClust",
  funParams = funParams,
  env = "local",
  runInplace = TRUE,
  taskID = 1
)
```


```{r}
x <- cciaObj$cciaObjectsByUIDs(c("VRZml0"))[[1]]

labels <- cciaEnv()$LabelPropsUtils(
  task_dir = x$persistentObjectDirectory(),
  labels_file = file.path("labelProps", "default.clust.h5ad")
)$label_props_view()

labelsDT <- labels$as_df()
labels$close()
```


```{r}
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  # pID = pID, uID = "jGcZHp", versionID = versionID, initReactivity = FALSE # TEHO Mouse
  # pID = pID, uID = "akAcC1", versionID = versionID, initReactivity = FALSE # Set
  pID = pID, uID = "pN5QJi", versionID = versionID, initReactivity = FALSE # Human LN
  # pID = pID, uID = "oUPuBX", versionID = versionID, initReactivity = FALSE # Human HCC
  # pID = pID, uID = "DPKHD7", versionID = versionID, initReactivity = FALSE # THEO Human
)

# get experimental info
exp.info <- as.data.table(cciaObj$summary(
  withSelf = FALSE, fields = c("Attr")
))

uIDs <- names(cciaObj$cciaObjects())
# uIDs <- exp.info[Include == "Y"]$uID
# uIDs <- exp.info[Analyse == "Y"]$uID
uIDs <- exp.info[Stain == "A"]$uID
# uIDs <- exp.info[Include == "Y" & Stain == "B"]$uID
# uIDs <- exp.info[Include == "Y" & Stain == "A"]$uID
# uIDs <- c("XFHvVd")

pops <- cciaObj$popPaths(uIDs = uIDs[[1]], popType = "clust", includeFiltered = TRUE)
pops <- pops[is.na(stringr::str_match(pops, "/"))]
# pops <- c("root")

# regionPops <- cciaObj$popPaths(uIDs = uIDs[[1]], popType = "region", includeFiltered = TRUE)

# get objects from selected set
regionsDT <- cciaObj$popDT(
  popType = "region",
  uIDs = uIDs,
  # pops = regionPops,
  includeFiltered = TRUE,
  completeDT = TRUE, replaceNA = TRUE,
  filterMeasures = c("regions")
)
unique(regionsDT$regions)

# get objects from selected set
# popDT <- cciaObj$popDT(
#   popType = "clust",
#   uIDs = uIDs,
#   pops = pops,
#   includeFiltered = TRUE,
#   completeDT = TRUE, replaceNA = TRUE
# )
```

```{r}
# save out populations
summaryDF <- regionsDT %>%
  group_by(uID, pop) %>%
  summarise(n = n()) %>%
  mutate(freq = n/sum(n) * 100) %>%
  ungroup() %>%
  complete(uID, pop, fill = list(freq = 0)) %>%
  left_join(exp.info)

write.csv(summaryDF, file = file.path(anaDir, "freq_regions.csv"))
```

```{r}
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "pN5QJi", versionID = versionID, initReactivity = FALSE # Human LN
)

# get experimental info
exp.info <- as.data.table(cciaObj$summary(withSelf = FALSE, fields = c("Attr")))
uIDs <- exp.info[Stain == "A"]$uID

pops <- cciaObj$popPaths(uIDs = uIDs[[1]], popType = "clust", includeFiltered = TRUE)
pops <- pops[is.na(stringr::str_match(pops, "/"))]

spe <- cciaObj$spe(
  popType = "clust", pops = pops, uIDs = uIDs, includeFiltered = TRUE, completeDT = TRUE)
```

```{r fig_spiat_pops, fig.height=6, fig.width=10}
x <- cciaObj$cciaObjects(uIDs = uIDs)[[1]]
imPopMap <- x$imPopMap(popType = "clust", includeFiltered = TRUE)
popColours <- sapply(imPopMap, function(x) {x$colour})
names(popColours) <- sapply(imPopMap, function(x) {x$path})

phenotypes <- unique(spe[[1]]$Phenotype)
popColours <- popColours[names(popColours) %in% phenotypes]

SPIAT::plot_cell_categories(
  spe[[1]], phenotypes, unname(popColours), "Phenotype", cex = 0.1) +
  coord_fixed()

ggsave(file.path(anaDir, "SPIAT", paste0(names(spe)[[1]], "_map.png")), height = 4, width = 6)
```

```{r}
gradient_pos <- seq(50, 500, 50) # radii

spiat.entropies <- lapply(spe, function(x) {
  as.data.table(SPIAT::entropy_gradient_aggregated(
    x, cell_types_of_interest = pops,
    feature_colname = "Phenotype", radii = gradient_pos)$gradient_df)
})

entropiesDT <- rbindlist(spiat.entropies, idcol = "uID")
```

```{r fig_spe_entity, fig.height=10, fig.width=12}
for (i in unique(entropiesDT$uID)) {
  datToPlot <- entropiesDT[uID == i] %>%
    pivot_longer(
      cols = starts_with("Pos_"),
      names_to = "radius",
      names_pattern = ".*_(.*)",
      values_to = "value") %>%
    mutate(radius = as.numeric(radius))
  
  p1 <- ggplot(datToPlot, aes(radius, value, color = Celltype2, fill = Celltype2, group = Celltype2)) +
    theme_classic() +
    # geom_line(size = 2) +
    geom_smooth() +
    facet_wrap(.~Celltype1, ncol = 3) +
    scale_color_brewer(palette = "Set1") +
    scale_fill_brewer(palette = "Set1") 
  
  print(p1)
}
```

```{r fig_regions, fig.height=20, fig.width=30}
# init ccia object
# cciaObj <- initCciaObject(
#   pID = "Fu1zEC", uID = "pN5QJi", versionID = 1, initReactivity = FALSE # LN Set
# )
# 
# # get experimental info
# exp.info <- as.data.table(cciaObj$summary(withSelf = FALSE, fields = c("Attr")))
# uIDs <- exp.info[Stain == "A"]$uID
# 
# popDT <- cciaObj$popDT(
#   popType = "clust", uIDs = uIDs, includeFiltered = TRUE, completeDT = TRUE)

# get colours
colPal2 <- randomcoloR::distinctColorPalette(length(unique(popDT$clusters)))

# for (x in cciaObj$cciaObjects(uIDs = uIDs)) {
#   # plot UMAP
#   p1 <- ggplot(popDT[uID == x$getUID()], aes(UMAP_1, UMAP_2)) +
#     geom_point(aes(color = clusters)) +
#     theme_classic() +
#     scale_color_manual(name = "Clusters", values = colPal) +
#     guides(colour = guide_legend(override.aes = list(size=5))) +
#     ggtitle(x$getUID())
# 
#   print(p1)
# }

p1 <- ggplot(popDT, aes(UMAP_1, UMAP_2)) +
  geom_point(aes(color = clusters), size = 0.5) +
  theme_classic() +
  scale_color_manual(name = "Clusters", values = colPal2) +
  guides(colour = guide_legend(override.aes = list(size=5))) +
  facet_wrap(.~uID, ncol = 4) +
  theme(
    strip.background = element_blank(),
    strip.text.x = element_blank(),
    legend.position = "none",
    axis.text = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )

ggsave(file.path(anaDir, "umaps.tiff"), width = 6, height = 4, plot = p1)
```

```{r}
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# could you cluster clusters across samples?
# because I have no idea what most of these markers are ...
# get matrix
popMat <- t(adataMatFromPopDT(popDT[
  , -c("clusters.group", "label", "pop", paste0("centroid_", c("x", "y", "z", "t")))],
  popKeys = c("uID", "clusters")))
```

```{r}
# now plot on UMAP?
popUMAP <- umap::umap(popMat)
```


```{r fig_umap_clustered, fig.height=6, fig.width=6}
umapDF <- as.data.frame(popUMAP$layout) %>%
  rownames_to_column() %>%
  separate_wider_delim(cols = rowname, delim = ".", names = c("uID", "clusters"))

# run KMeans to get the same number of clusters?
# clustToFind <- max(as.numeric(umapDF$clusters))
# clustToFind <- median((umapDF %>% group_by(uID) %>% summarise(max.group = max(as.numeric(as.character(clusters)))))$max.group)
clustToFind <- median((umapDF %>% group_by(uID) %>% summarise(max.group = max(as.numeric(as.character(clusters)))))$max.group)
# clustToFind <- 9
# clusters <- kmeans(popMat, clustToFind)
clusters <- dbscan::dbscan(
  umapDF[, c("V1", "V2")], eps = 2, minPts = 4, borderPoints = TRUE)

umapDF$clusters.group <- as.factor(clusters$cluster)
# umapDF <- .mapClustNames(umapDF, clusterMapping, "clusters.group", removeNone = FALSE)

colNames <- unique(umapDF$clusters.group)
# colPal <- randomcoloR::distinctColorPalette(length(unique(umapDF$uID)))
colPal <- randomcoloR::distinctColorPalette(length(unique(umapDF$clusters.group)))
# colPal <- randomcoloR::distinctColorPalette(length(unique(colNames)))
# colPal <- randomcoloR::distinctColorPalette(length(colNames))
names(colPal) <- colNames

# save to set
# saveRDS(umapDF, file = file.path(cciaObj$persistentObjectDirectory(), "data", "cluster.groups.rds"))
# umapDF <- readRDS(file.path(cciaObj$persistentObjectDirectory(), "data", "cluster.groups.rds"))

ggplot(umapDF, aes(V1, V2)) +
  # geom_point(aes(color = uID), size = 2) +
  geom_point(aes(color = clusters.group), size = 2) +
  # geom_point(aes(color = clusters.name), size = 3) +
  theme_classic() +
  # scale_color_brewer(name = NULL, palette = "Set1") +
  scale_color_manual(values = colPal) +
  # scale_color_manual(values = clusterColours) +
  # viridis::scale_colour_viridis() +
  theme(
    axis.text = element_text(size = 15),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    legend.justification = "right",
    legend.position = "bottom",
    # legend.position = "none",
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
    )
```

```{r}
# push cluster groups back to DF
# popDT[, clusters.group := NULL]
# for (i in rownames(umapDF)) {
#   x <- umapDF[i, ]
#   
#   popDT[uID == x$uID & clusters == x$clusters, clusters.group := as.numeric(x$clusters.group)]
# }
# 
# popDT[, clusters.group := as.factor(clusters.group)]

clusterMapping <- list(
  CD3 = c(8),
  CD68.CD11c = c(3),
  CD11b = c(5),
  EPCAM = c(6),
  PD1 = c(10),
  HNFa = c(2),
  Others = c(12,7,11,4,9,1)
)

clusterColours <- list(
  CD3 = "#ff1493",
  CD68.CD11c = "#00bfff",
  CD11b = "#ffd700",
  EPCAM = "#e9967a",
  PD1 = "#00fa9a",
  # HNFa = "white",
  HNFa = "black",
  Others = "#d8bfd8"
)

popDT <- .mapClustNames(popDT, clusterMapping, "clusters.group", removeNone = FALSE)
```

```{r}
      - "#ff0000"
      - "#ff8c00"
      - "#ffd700"
      - "#0000cd"
      - "#00ff00"
      - "#00fa9a"
      - "#4169e1"
      - "#e9967a"
      - "#00ffff"
      - "#00bfff"
      - "#d8bfd8"
      - "#ff00ff"
      - "#f0e68c"
      - "#ff1493"
      - "#ee82ee"
```

```{r}
# show pops on image
Sys.setenv(KMP_DUPLICATE_LIB_OK = "TRUE")
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = TRUE, initJupyter = TRUE)
# cciaUse("~/cecelia/dev", initConda = FALSE, initJupyter = TRUE)

# start ipython kernel
viewer <- NapariUtils$new()

viewer$initNapari()
```

```{r}
generateScreenshots <- function(
    anaDir, cciaObj, uIDs, channelsToShow, popsToShow, fileAttrs,
    valueName = "default", matchChannels = FALSE, save = TRUE, ...) {
  # go through images and save movies
  for (x in cciaObj$cciaObjects(uIDs = uIDs)) {
    i <- x$getUID()
    viewer$setTaskDir(x$persistentObjectDirectory())
    
    # get channels
    imChannelNames <- x$imChannelNames()
    
    # match channels
    curChannelsToShow <- channelsToShow
    if (matchChannels == TRUE) {
      # TODO this is very crude    
      names(curChannelsToShow) <- sapply(names(channelsToShow), function(y) {
        imChannelNames[!is.na(stringr::str_extract(imChannelNames, y))][[1]]
      })
    }

    # get cmaps
    layersVisible <- imChannelNames %in% names(curChannelsToShow)
    channelColormaps <- rep("gray", length(layersVisible))
    
    for (j in names(curChannelsToShow))
      channelColormaps[[which(imChannelNames == j)]] <- curChannelsToShow[[j]]
    
    names(layersVisible) <- seq(length(layersVisible))
    
    # open image
    viewer$openImage(
      x$imFilepath(),
      imChannelNames = imChannelNames,
      show3D = FALSE,
      layersVisible = layersVisible,
      channelColormaps = channelColormaps
      # multiscales = as.integer(1)
    )
    
    # add populations
    popsAdded <- list()
    
    # show on image
    removePrevious <- TRUE
    pointsSize <- max(x$omeXMLPixels()$SizeX, x$omeXMLPixels()$SizeY) * 0.005
    # pointsSize <- max(x$omeXMLPixels()$SizeX, x$omeXMLPixels()$SizeY) * 0.003
    
    # add populations to popMap
    for (popType in names(popsToShow)) {
      popMap <- x$imPopMap(popType, includeFiltered = TRUE)
      curPops <- popsToShow[[popType]]
    
      # get popIDs for paths
      popIDs <- lapply(names(curPops), function(y) x$popIDsByAttr(popType, "path", y, includeFiltered = TRUE))
      names(popIDs) <- names(curPops)
      
      for (pop in names(curPops)) {
        x$editPopColour(popType, popIDs[[pop]], curPops[[pop]], includeFiltered = TRUE)
      }
    
      # save to disk
      x$savePopMap(popType, includeFiltered = TRUE)
      x$savePops(popType, pops = names(popIDs),
                 purge = TRUE, includeFiltered = TRUE)
      
      for (j in names(curPops)) {
        viewer$showPopMapping(popType, valueName = popMap[[popIDs[[j]]]]$valueName,
                              pops = list(j), pointsSize = pointsSize,
                              execInteractive = FALSE, removePrevious = removePrevious)
        removePrevious <- FALSE
      }
    }
    
    # adjust channel intensities
    browser()
  
    # save screenshot
    if (save == TRUE)
      viewer$saveScreenshot(
        file.path(anaDir, paste0(
          paste(sapply(fileAttrs, function(i) x$getCciaAttr(i)), collapse = "_"), "_", i, ".png")),
        ...)
  }
}

# init ccia object
cciaObj <- initCciaObject(
  # pID = pID, uID = "jGcZHp", versionID = versionID, initReactivity = FALSE
  pID = pID, uID = "DPKHD7", versionID = versionID, initReactivity = FALSE
)

exp.info <- as.data.table(cciaObj$summary(withSelf = FALSE, fields = c("Attr")))
uIDs <- names(cciaObj$cciaObjects())
# uIDs <- exp.info[Analyse == "Y"]$uID

channelsToShow <- list(
  "DAPI" = "bop blue",
  "GFAP" = "gray",
  # "CD3" = "magenta", # THEO Mouse
  "CD3e" = "magenta", # THEO Human
  "Iba-1" = "green" # THEO Mouse
  )
popsToShow <- list(
  # flow = list(
  #   "/cells/O1/O2/B220" = "#327EBA",
  #   "/cells/O1/CD11c" = "yellow",
  #   "/cells/O1/XCR1" = "green",
  #   "/cells/P14" = "magenta",
  #   "/cells/P14/clustered" = "white"
  # ),
  # branch = list(
  #   "LCMV.branch/two" = "red"
  # )
)
fileAttrs <- c("Treatment")

generateScreenshots(
  file.path(anaDir, "images"), cciaObj, uIDs, channelsToShow,
  # file.path(anaDir, "pops"), cciaObj, uIDs, list(),
  # list(), fileAttrs, size = list(800, 800), save = TRUE)
  list(), fileAttrs, size = NULL, save = TRUE)
  # popsToShow, fileAttrs, size = list(1200, 1200), save = TRUE)
```

```{r fig_umap, fig.height=4, fig.width=4}
dfToPlot <- copy(popDT) %>% drop_na(clusters) %>% left_join(exp.info)
# dfToPlot <- .mapClustNames(dfToPlot, clusterMapping, "clusters")

# get mean positions of clusters
meanClusterPos <- dfToPlot %>%
  group_by(clusters) %>%
  # group_by(clusters.name) %>%
  summarise(
    UMAP_1 = mean(UMAP_1),
    UMAP_2 = mean(UMAP_2)
  )
dfToPlot$clusters <- factor(dfToPlot$clusters)

# colPal <- randomcoloR::distinctColorPalette(length(unique(dfToPlot$clusters)))

markers <- colnames(dfToPlot)
markers <- markers[!markers %in% c(
  "label", "clusters", "UMAP_1", "UMAP_2", "pop", "centroid_y", "centroid_x", colnames(exp.info))]

for (i in markers) {
  print(i)
  
    # plot UMAP
  p1 <- ggplot(
    dfToPlot,
    aes(UMAP_1, UMAP_2)
    ) +
    # geom_point(aes(color = clusters)) +
    # geom_point(aes(color = Type)) +
    geom_point(aes(color = get(i))) +
    theme_classic() +
    # scale_color_brewer(name = NULL, palette = "Set1") +
    # scale_color_manual(values = colPal) +
    viridis::scale_colour_viridis() +
    theme(
      axis.text = element_text(size = 15),
      axis.text.x = element_blank(),
      axis.ticks.x = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks.y = element_blank(),
      legend.justification = "right",
      # legend.position = "bottom",
      legend.position = "none",
      axis.line = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank()
      ) +
    geom_label(
      data = meanClusterPos,
      aes(label = clusters),
      # aes(label = clusters.name),
      label.size = 0.25,
      color = "black"
      ) +
    ggtitle(i)
    # facet_grid(.~Type)
  
  print(p1)
}

ggsave(file.path(anaDir, "umap.tiff"), width = 2, height = 2)
# ggsave(file.path(anaDir, "umap_coating.tiff"), width = 2, height = 2)
# ggsave(file.path(anaDir, "umap_coating_treatment.tiff"), width = 4, height = 2.5)
```

```{r fig_umap, fig.height=5, fig.width=5}
# clusterMapping <- list(
#   `0` = c(0, 7),
#   `1` = c(1),
#   `2` = c(2),
#   `3` = c(3, 8),
#   `4` = c(4),
#   `5` = c(5),
#   `6` = c(6)
# )

colNames <- unique(popDT$clusters)
# colNames <- unique(umapDF$clusters.group)
colPal <- randomcoloR::distinctColorPalette(length(colNames))
names(colPal) <- colNames

for (x in cciaObj$cciaObjects(uIDs = uIDs)) {
  dfToPlot <- popDT[uID == x$getUID()] %>% drop_na(clusters) %>% left_join(exp.info)
  # dfToPlot <- popDT %>% drop_na(clusters) %>% left_join(exp.info)
  # dfToPlot <- popDT[!clusters %in% seq(0, 7)] %>% drop_na(clusters) %>% left_join(exp.info)
  # dfToPlot <- .mapClustNames(dfToPlot, clusterMapping, "clusters")
  
  # get mean positions of clusters
  meanClusterPos <- dfToPlot %>%
    group_by(clusters) %>%
    # group_by(clusters.group) %>%
    # group_by(clusters.name) %>%
    summarise(
      UMAP_1 = median(UMAP_1),
      UMAP_2 = median(UMAP_2)
    )
  dfToPlot$clusters <- factor(dfToPlot$clusters)
  # dfToPlot$clusters.group <- factor(dfToPlot$clusters.group)
  
  # colPal <- randomcoloR::distinctColorPalette(length(unique(dfToPlot$clusters)))
  # colPal <- randomcoloR::distinctColorPalette(length(unique(dfToPlot$clusters.name)))
  
  # plot UMAP
  p1 <- ggplot(
    dfToPlot,
    aes(UMAP_1, UMAP_2)
    ) +
    geom_point(aes(color = clusters)) +
    # geom_point(aes(color = clusters.group)) +
    # geom_point(aes(color = uID), alpha = 0.2) +
    # geom_point(aes(color = clusters.name)) +
    # geom_point(aes(color = Treatment)) +
    theme_classic() +
    # scale_color_brewer(name = "Clusters", palette = "Accent") +
    # scale_color_brewer(name = "Clusters", palette = "Set1") +
    scale_color_manual(name = "Clusters", values = colPal) +
    # scale_color_brewer(palette = "Set3", name = "Clusters") +
    # viridis::scale_colour_viridis() +
    theme(
      axis.text = element_text(size = 15),
      axis.text.x = element_blank(),
      axis.ticks.x = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks.y = element_blank(),
      legend.justification = "left",
      legend.position = "bottom",
      # legend.position = "none",
      axis.line = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank()
      ) +
    # geom_label(
    #   data = meanClusterPos,
    #   aes(label = clusters),
    #   # aes(label = clusters.group),
    #   # aes(label = clusters.name),
    #   label.size = 0.25,
    #   color = "black"
    #   ) 
    guides(colour = guide_legend(override.aes = list(size=5)))
    # facet_wrap(.~uID, nrow = 1)
    # facet_wrap(.~uID, ncol = 4)
  
  # ggsave(file.path(anaDir, "UMAPs", paste0("umap-", cciaObj$getUID(), ".png")), width = 4, height = 4, plot = p1)
  ggsave(file.path(anaDir, "UMAPs", paste0("umap-", exp.info[uID == x$getUID()]$Treatment, "-", x$getUID(), ".png")), width = 4, height = 5, plot = p1)
  # ggsave(file.path(anaDir, "UMAPs", paste0("umap-", x$getUID(), ".png")), width = 4, height = 4, plot = p1)
}
# ggsave(file.path(anaDir, "umap-uids.png"), width = 15, height = 15, plot = p1)
# ggsave(file.path(anaDir, "umap_coating.png"), width = 2, height = 2)
# ggsave(file.path(anaDir, "umap_coating_treatment.png"), width = 4, height = 2.5)
```

```{r fig_pops_im, fig.height=40, fig.width=20}
# popClusters <- unique(popDT$clusters)
# colPal <- randomcoloR::distinctColorPalette(length(popClusters))
# names(colPal) <- popClusters

# uidsToFlip <- c("XFHvVd", "Lm5FOH")

for (i in unique(popDT$uID)) {
  x <- popDT[uID == i]
  print(paste(">>", i))
  
  # go through clusters
  # for (j in unique(x$clusters)) {
  for (j in unique(x$pop)) {
    print(paste(">", j))
    
    # plot out
    p1 <- ggplot(x,
                 # aes(centroid_x, centroid_y, color = as.factor(clusters))) +
                 aes(centroid_x, -centroid_y)) +
                 # aes(centroid_x, centroid_y, color = as.factor(clusters.name))) +
      theme_classic() +
      plotThemeDark(angle = 0) +
      # geom_point(size = 0.2) +
      geom_point(data = x, size = 0.2, color = "darkgrey") +
      # geom_point(data = x[clusters == j], size = 0.2, color = "red") +
      geom_point(data = x[pop == j], size = 0.2, color = "red") +
      # scale_color_brewer(palette = "Set3", name = "Clusters") +
      # scale_color_manual(values = colPal, name = "Clusters") +
      theme(
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        # legend.position = "none"
        legend.position = "bottom",
        strip.text = element_text(colour = "black")
      ) + coord_fixed()
      # guides(colour = guide_legend(override.aes = list(size=10))) +
      # facet_wrap(.~uID, ncol = 4)
      # facet_grid(clusters~uID)
      # facet_grid(clusters.name~uID)
      # ggtitle(i)
    
    # if (i %in% uidsToFlip) {
    #   p1 <- p1 + coord_flip()
    # }
    
    # print(p1)
    
    # ggsave(file.path(anaDir, "clusters_im.png"), width = 40, height = 20, plot = p1)
    # ggsave(file.path(anaDir, "maps", paste0("clusters-", exp.info[uID == i]$Treatment, "-", i, "-", j, ".png")), width = 5, height = 5, plot = p1)
    ggsave(file.path(anaDir, "maps", paste0("pop-", exp.info[uID == i]$Treatment, "-", i, "-", j, ".png")), width = 5, height = 5, plot = p1)
  }
}
```

```{r}
popDefs <- list()

# print out population definitions
for (x in cciaObj$cciaObjects(uIDs = uIDs)) {
  popDef <- sapply(
    x$imPopMap("clust", includeFiltered = TRUE),
    function(y) if (length(y$filterValues) > 0)paste(y$path, paste(y$filterValues, collapse = ","), sep = ": "))
  
  popDefs[[x$getUID()]] <- as.data.table(t(as.data.frame(popDef[lengths(popDef) > 0])))
}

write.csv(rbindlist(popDefs, idcol = "uID"), file.path(anaDir, "population_definitions.csv"))
```


```{r fig_clusters_freq, fig.height=5, fig.width=15}
# show frequency of populations for treatments
summaryDF <- popDT %>%
  # group_by(uID, clusters) %>%
  group_by(uID, pop) %>%
  summarise(n = n()) %>%
  mutate(freq = n/sum(n) * 100) %>%
  ungroup() %>%
  left_join(exp.info)

levelOrder <- unique(exp.info[uID %in% uIDs]$Treatment)
summaryDF$Treatment <- factor(summaryDF$Treatment, levels = levelOrder)

colPal <- as.list(randomcoloR::distinctColorPalette(length(unique(summaryDF$pop))))

# plot
ggplot(summaryDF,
       # aes(paste(Treatment, uID), freq, fill = clusters)) +
       aes(paste(Treatment, uID), freq, fill = pop)) +
  theme_classic() +
  geom_col() +
  coord_flip() +
  # scale_fill_manual(values = colPal, name = "Clusters") +
  # xlab("") + ylab("Clusters (%)")
  scale_fill_manual(values = colPal, name = "Pops") +
  xlab("") + ylab("Population (%)")

# ggsave(file.path(anaDir, "clusters-freq.png"), width = 8, height = 4)
ggsave(file.path(anaDir, "pops-freq.png"), width = 8, height = 5)
```


```{r fig_pops_im, fig.height=10, fig.width=10}
# uidsToFlip <- c("XFHvVd", "Lm5FOH")

x <- cciaObj$cciaObjects(uIDs = uIDs[[1]])[[1]]
xres <- x$omeXMLPixelRes()$x

# # imPopMap <- x$imPopMap(popType = "clust", includeFiltered = TRUE)
# imPopMap <- x$imPopMap(popType = "region", includeFiltered = TRUE)
# popColours <- sapply(imPopMap, function(x) {x$colour})
# names(popColours) <- sapply(imPopMap, function(x) {x$path})

# channelCols <- unname(x$imChannelNames(correctChannelNames = TRUE))
# channelCols <- channelCols[channelCols %in% colnames(popDT)]
# channelCols <- channelCols[!channelCols %in% c("HOECHST")]
# dfToPlot <- as.data.table(popDT %>%
#   pivot_longer(cols = channelCols, names_to = "channels", values_to = "intensity"))

# clusters.groups <- popDT %>% group_by(uID, clusters.group) %>% summarise() %>% group_map(~.x)
# clusters.groups <- c(unname(unlist(Reduce(intersect, clusters.groups))), 10) # HCC
# clusters.groups <- c(unname(unlist(Reduce(intersect, clusters.groups))), 4, 11, 5, 9) # THEO MOUSE
# clusters.groups <- c(unname(unlist(Reduce(intersect, clusters.groups))))

popClusters <- unique(popDT$clusters)
# popClusters <- unique(popDT$clusters.group)
# popClusters <- unique(popDT$clusters.name)
# popClusters <- unique(regionsDT$regions)
# popClusters <- clusters.groups
colPal <- as.list(randomcoloR::distinctColorPalette(length(popClusters)))
# colPal <- list("#4169e1", "#ff1493", "#ffd700", "#00fa9a", "#ff8c00", "#d8bfd8", "#e9967a", "#ff0000")
# colPal <- list("#ff0000", "#ff8c00", "#ffd700", "#0000cd", "#00ff00", "#00fa9a",
#                "#4169e1", "#e9967a", "#00ffff", "#ff1493", "#d8bfd8", "#ff00ff",
#                "#f0e68c", "#ee82ee", "#00bfff")
names(colPal) <- popClusters

# for (i in unique(dfToPlot$uID)) {
# for (i in unique(regionsDT$uID)) {
for (i in unique(popDT$uID)) {
# for (i in c("XFHvVd")) {
# for (i in c("SiaqvM")) { 
  # x <- popDT[uID == i & clusters.name == "CD68.CD11c"] %>%
  #                left_join(spatialDT[pop.from == "CD68.CD11c" & pop.to == "CD3"], by = c("uID", "label" = "from"))
  # z <- x[is.na(pop.to)]
  # y <- x[!is.na(pop.to)]
  # 
  # if (nrow(x) == 0) next
  
  # plot out
  # p1 <- ggplot(dfToPlot[uID == i],
  # p1 <- ggplot(spatialDT[uID == i & pop.from == "CD3" & pop.to == "CD11c"] %>%
  # p1 <- ggplot(x,
  p1 <- ggplot(popDT[uID == i],
  # p1 <- ggplot(popDT[uID == i & clusters.group %in% clusters.groups],
  # p1 <- ggplot(popDT[uID == i &
  # p1 <- ggplot(popDT[uID == i & !clusters %in% seq(0,7) &
                       # centroid_y >= 1849/xres & centroid_y <= 7196/xres &
                       # centroid_x >= 3535/xres & centroid_x <= 7920/xres],
                       # centroid_y >= 2096/xres & centroid_y <= 4406/xres &
                       # centroid_x >= 1200/xres & centroid_x <= 4597/xres],
  # p1 <- ggplot(popDT[uID == i & !clusters.group %in% c(3, 12, 5, 13)],
  # p1 <- ggplot(regionsDT[uID == i],
               aes(centroid_x, -centroid_y, color = as.factor(clusters))) +
               # aes(centroid_x, -centroid_y, color = as.factor(regions))) +
               # aes(centroid_x, -centroid_y, color = GFAP)) +
               # aes(centroid_x, -centroid_y, color = as.factor(pop))) +
               # aes(centroid_x, -centroid_y, color = as.factor(clusters.group))) +
               # aes(centroid_x, -centroid_y, color = as.factor(clusters.name))) +
               # aes(centroid_x, -centroid_y)) +
    theme_classic() +
    plotThemeDark(angle = 0) +
    # geom_point(size = 0.1, alpha = 1) +
    geom_point(size = 0.5, alpha = 1) +
    # geom_point(aes(color = intensity), size = 0.1, alpha = 1) +
    # geom_point(data = z, color = .flowColours(z$centroid_x, z$centroid_y, colramp = viridisLite::inferno)) +
    # geom_point(data = y, color = .flowColours(y$centroid_x, y$centroid_y, colramp = viridisLite::inferno)) +
    # scale_color_brewer(palette = "Accent", name = "Clusters") +
    # scale_color_brewer(palette = "Set1", name = "Clusters") +
    scale_color_manual(values = colPal, name = "Clusters") +
    # scale_color_manual(values = clusterColours, name = "pop") +
    # scale_color_manual(values = popColours, name = "Pop") +
    # scale_color_manual(values = c("NA" = "grey", "CD3" = "red"), name = "Pop") +
    # viridis::scale_color_viridis() +
    # scale_color_gradient2(low = "blue", high = "red", mid = "white", midpoint = 25) +
    theme(
      axis.text.x = element_blank(),
      axis.ticks.x = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks.y = element_blank(),
      axis.line = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      # legend.position = "none",
      legend.position = "bottom",
      strip.text = element_text(colour = "black")
    ) + coord_fixed() +
    guides(colour = guide_legend(override.aes = list(size=10))) 
    # facet_wrap(.~clusters, ncol = 4)
    # facet_wrap(.~clusters.group, ncol = 6)
    # facet_wrap(.~clusters.name, ncol = 6)
    # facet_wrap(.~regions, ncol = 6)
    # facet_grid(.~pop.to)
    # facet_wrap(.~channels, ncol = 8)
    # facet_grid(clusters~uID)
    # facet_grid(clusters.name~uID)
    # ggtitle(i)
  
  ggsave(file.path(anaDir, "maps", paste0("clusters_", i, ".png")), width = 6, height = 8, plot = p1)
  # ggsave(file.path(anaDir, "maps", paste0("clusters_group_", i, ".png")), width = 6, height = 8, plot = p1)
  # ggsave(file.path(anaDir, "maps", paste0("clusters_name_", i, ".png")), width = 12, height = 16, plot = p1)
  # ggsave(file.path(anaDir, "maps", paste0("labels_", i, ".png")), width = 30, height = 20, plot = p1)
  # ggsave(file.path(anaDir, "maps", paste0("clusters_", i, "-dist.png")), width = 8, height = 5, plot = p1)
  # ggsave(file.path(anaDir, "maps", paste0("clusters_", i, "-zoom.png")), width = 12, height = 16, plot = p1)
  # ggsave(file.path(anaDir, "maps", paste0("pops_", i, ".png")), width = 12, height = 12, plot = p1)
  # ggsave(file.path(anaDir, "maps", paste0("regions_", i, ".png")), width = 6, height = 8, plot = p1)
  # ggsave(file.path(anaDir, "maps", paste0("clusters_facet_", i, ".png")), width = 20, height = 20, plot = p1)
}
```

```{r}
# plot all images out with the same populations
x <- cciaObj$cciaObjects(uIDs = uIDs[[1]])[[1]]
xres <- x$omeXMLPixelRes()$x

clusters.groups <- popDT %>% group_by(uID, clusters.group) %>% summarise() %>% group_map(~.x)
# clusters.groups <- c(unname(unlist(Reduce(intersect, clusters.groups))), 10) # HCC
clusters.groups <- c(unname(unlist(Reduce(intersect, clusters.groups))), 1, 4, 11, 5, 9, 14) # THEO MOUSE
# clusters.groups <- c(unname(unlist(Reduce(intersect, clusters.groups))))
# clusters.groups <- clusters.groups[clusters.groups != 1]

colPal <- list("darkgrey", "#ff0000", "#ff8c00", "#ffd700", "#0000cd", "#00ff00", "#00fa9a",
               "#4169e1", "#e9967a", "#00ffff", "#ff1493", "#d8bfd8", "#ff00ff",
               "#f0e68c", "#ee82ee", "#00bfff")

p1 <- ggplot(popDT[clusters.group %in% clusters.groups],
# p1 <- ggplot(popDT,
             aes(centroid_x, -centroid_y, color = as.factor(clusters.group))) +
             # aes(centroid_x, -centroid_y, color = as.factor(clusters.name))) +
  theme_classic() +
  plotThemeDark(angle = 0) +
  geom_point(size = 0.1, alpha = 1) +
  # geom_point(data = popDT[clusters.name == "Others"], size = 0.1, alpha = 1) +
  # geom_point(data = popDT[clusters.name != "Others"], size = 0.1, alpha = 1) +
  geom_point(data = popDT[clusters.group == 1 & clusters.group %in% clusters.groups], size = 0.1, alpha = 1) +
  geom_point(data = popDT[clusters.group != 1 & clusters.group %in% clusters.groups], size = 0.1, alpha = 1) +
  # scale_color_brewer(palette = "Set3", name = "Clusters") +
  # scale_color_manual(values = clusterColours, name = "pop") +
  scale_color_manual(values = colPal, name = "Clusters") +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    # legend.position = "none",
    legend.position = "bottom",
    strip.text = element_text(colour = "black")
  ) + coord_fixed() +
  guides(colour = guide_legend(override.aes = list(size=10))) +
  facet_wrap(.~uID, ncol = 4)
  
ggsave(file.path(anaDir, "maps", paste0("clusters_group_", cciaObj$getUID(), ".png")), width = 32, height = 18, plot = p1)
# ggsave(file.path(anaDir, "maps", paste0("clusters_name_", cciaObj$getUID(), ".png")), width = 24, height = 18, plot = p1)
```


```{r}
# is there any spatial interaction?
spatialDT <- cciaObj$spatialDT(valueName = "default", uIDs = uIDs)

# join pops
# spatialDT[popDT[, c("uID", "label", "clusters.group")],
spatialDT[popDT[, c("uID", "label", "clusters.name")],
          on = c("uID", "to" = "label"),
          # pop.to := clusters.group]
          pop.to := clusters.name]
# spatialDT[popDT[, c("uID", "label", "clusters.group")],
spatialDT[popDT[, c("uID", "label", "clusters.name")],
          on = c("uID", "from" = "label"),
          # pop.from := clusters.group]
          pop.from := clusters.name]

# filter same type associations
spatialDT <- spatialDT[pop.to != pop.from]
```

```{r fig_interactions, fig.height=2, fig.width=6}
# get "i-niches" and determine a given number of defined clusters
freqInt <- spatialDT %>%
  group_by(uID, pop.from, pop.to) %>%
  summarise(n = n()) %>%
  mutate(freq = n/sum(n) * 100) %>%
  drop_na() %>%
  ungroup() %>%
  complete(uID, pop.from, pop.to, fill = list(freq = 0)) %>%
  left_join(exp.info)

# show interactions of individuals
for (i in unique(freqInt$pop.from)) {
  x <- freqInt %>% dplyr::filter(pop.from == i) 
  
  p1 <- ggplot(x, aes(as.factor(pop.to), freq)) +
    theme_classic() +
    geom_boxplot(outlier.alpha = 0) +
    geom_jitter(
      # position = position_jitterdodge(jitter.width = 0.10), alpha = 1.0) +
      width = 0.3, alpha = 0.6) +
    # stat_summary(fun=mean, geom="point", size=3, shape=18, color="Magenta") +
    ylab("Population (%)") + xlab("") +
    # ylim(0, 1) +
    scale_fill_brewer(name = NULL, palette = "Set3") +
    ggtitle(i)

  print(p1)
}
```


```{r fig_interactions_heat, fig.height=7, fig.width=7}
freqInt <- spatialDT %>%
  group_by(pop.from, pop.to) %>%
  summarise(n = n()) %>%
  mutate(freq = n/sum(n) * 100) %>%
  drop_na() %>%
  ungroup() %>%
  complete(pop.from, pop.to, fill = list(freq = 0))

ggplot(freqInt, aes(as.factor(pop.from), as.factor(pop.to))) +
  theme_classic() +
  geom_tile(aes(fill = freq), colour = "white", size = 0.5) +
  viridis::scale_fill_viridis(
    breaks = c(0, 1),
    labels = c(0, 1)
  ) +
  theme(
    legend.position = "none",
    legend.title = element_blank(),
    legend.text = element_text(size = 20),
    legend.key.size = unit(10, "mm"),
    axis.text.y = element_text(size = 20),
    axis.text.x = element_text(size = 20, angle = 45, hjust = 1, vjust = 1),
    # axis.text.x = element_text(size = 20),
    strip.background = element_blank(),
    strip.text.x = element_blank()
    ) + xlab("") + ylab("")
```

```{r}
# clustMarkers <- colnames(dfToPlot)
clustMarkers <- colnames(popDT)
# clustMarkers <- colnames(regionsDT)
clustMarkers <- clustMarkers[!clustMarkers %in% c(
  "label", "clusters", "clusters.id", "clusters.name", "clusters.group", "regions",
  "UMAP_1", "UMAP_2", "pop", "centroid_y", "centroid_x", colnames(exp.info))]

# define populations
clust <- list(
  init = c("root")
  # regions = c("root")
)

markers <- list(
  init = clustMarkers
  # regions = clustMarkers
)

groupVars <- list(
  # init = "clusters.name"
  init = "clusters"
  # init = "clusters.group"
  # regions = "regions"
)

combinePops <- list(
  init = list()
  # regions = list()
)

plotSizes <- list(
  # init = c(6, 10)
  # init = c(6, 10)
  # init = c(8, 12) # THEO Mouse
  init = c(8, 16) # THEO Human
  # regions = c(4, 10)
)
```

```{r fig_clusters_heat, fig.height=4, fig.width=4}
for (i in names(clust)) {
  popsToShow <- clust[[i]]
  markersToShow <- .flowCorrectChannelNames(markers[[i]])
  markersToShow <- markersToShow[!startsWith(markersToShow, "flow.")]
  groupVar <- groupVars[[i]]
  popsToCombine <- combinePops[[i]]
  plotSize <- plotSizes[[i]]
  
  markersToShowNoDye <- stringr::str_extract(markersToShow, ".*(?=_)")
  
  for (j in unique(popDT$uID)) {
    y <- popDT[uID == j]
    
    colsToShow <- c("label", "clusters", "clusters.name", "clusters.group", "pop", "regions", markersToShow)
    colsToShow <- colsToShow[colsToShow %in% colnames(y)]
    # colsToShow <- colsToShow[colsToShow %in% colnames(regionsDT)]
    # colsToShow <- colsToShow[colsToShow %in% colnames(dfToPlot)]
    
    # show heatmap
    popLymphs <- y[pop %in% popsToShow, ..colsToShow]
    # popLymphs <- regionsDT[pop %in% popsToShow, ..colsToShow]
    # popLymphs <- dfToPlot[pop %in% popsToShow, ..colsToShow]
    # popLymphs <- popDT[, ..colsToShow]
    
    # combine pops
    if (length(popsToCombine) > 0) {
      for (j in names(popsToCombine))
        popLymphs[pop %in% popsToCombine[[j]], pop := j]
    }
    
    markersToPlot <- popLymphs %>%
      pivot_longer(
        cols = all_of(markersToShow), names_to = "marker", values_to = "value"
        )
    markersToPlot$marker <- as.factor(markersToPlot$marker)
    
    # show heatmap for clusters
    markersSummary <- markersToPlot %>%
      drop_na() %>%
      group_by(marker, get(groupVar)) %>%
      summarise(mean = mean(value, rm.na = TRUE)) %>%
      mutate(norm = (mean - min(mean))/(max(mean) - min(mean))) %>%
      replace_na(list(norm = 0)) %>%
      rename_with(~ c(groupVar), all_of(c("get(groupVar)")))
    
    # find order for clusters and markers
    markersMatrix <- markersSummary %>%
      pivot_wider(names_from = marker, values_from = norm, id_cols = groupVar)
    
    hcClusters <- hclust(dist(as.matrix(markersMatrix %>% dplyr::select(-all_of(groupVar)))), "ave")
    hcMarkers <- hclust(dist(t(as.matrix(markersMatrix %>% dplyr::select(-all_of(groupVar))))), "ave")
    
    clustersOrder <- markersMatrix[[groupVar]][hcClusters$order]
    markerOrder <- colnames(markersMatrix)[2:length(colnames(markersMatrix))][hcMarkers$order]
    
    markersSummary[[groupVar]] <- factor(markersSummary[[groupVar]], levels = clustersOrder)
    
    p1 <- ggplot(markersSummary %>%
             mutate(
               # clusters = factor(clusters, levels = clustersOrder),
               # clusters = factor(clusters, levels = sort(as.numeric(levels(clustersOrder)))),
               marker = factor(marker, levels = markerOrder)
             ), aes(get(groupVar), marker)) +
      theme_classic() +
      # geom_tile(aes(fill = mean), colour = "white", size = 0.5) +
      geom_tile(aes(fill = norm), colour = "white", size = 0.5) +
      viridis::scale_fill_viridis(
        breaks = c(0, 1),
        labels = c(0, 1)
      ) +
      theme(
        legend.position = "none",
        legend.title = element_blank(),
        legend.text = element_text(size = 20),
        legend.key.size = unit(10, "mm"),
        axis.text.y = element_text(size = 20),
        axis.text.x = element_text(size = 20, angle = 45, hjust = 1, vjust = 1),
        # axis.text.x = element_text(size = 20),
        strip.background = element_blank(),
        strip.text.x = element_blank()
        ) + xlab("") + ylab("")
    
    plot(p1)
    
    # ggsave(file.path(anaDir, "heat", paste0("clusters_heat_", i, ".pdf")),
    ggsave(file.path(anaDir, "heat", paste0("clusters_heat_", exp.info[uID == j]$Treatment, "-", j, "-", i, ".pdf")),
           width = plotSize[[1]], height = plotSize[[2]])
  }
}
```

```{r}
# can you cluster and check if there is a difference in R and N?
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  # pID = pID, uID = "diRenc", versionID = versionID, initReactivity = FALSE # Set
  pID = pID, uID = "pN5QJi", versionID = versionID, initReactivity = FALSE # Set
)

regionUIDs <- names(cciaObj$cciaObjects())
# regionUIDs <- regionUIDs[!regionUIDs %in% c(
#   # "84kEUg", "LQFVZD", "6zeSiS"
#   "84kEUg"
# )]

# get pops
pops <- cciaObj$popPaths(uIDs = regionUIDs[[1]], popType = "clust", includeFiltered = TRUE)
pops <- pops[is.na(stringr::str_match(pops, "/"))]
# pops <- pops[!pops %in% c("none", "other")]

# run task
funParams <- list(
  uIDs = names(cciaObj$cciaObjects()),
  popType = c("clust"),
  # pops = c("root"),
  # pops = c("CD4"),
  pops = pops,
  objectMeasures = c(),
  # objectMeasures = c("nuc_Ki.67"),
  includeCounts = FALSE,
  # includeCounts = TRUE,
  includeTotalCount = FALSE,
  # includeTotalCount = TRUE,
  colName = "default",
  valueName = "default",
  excludeSelfConnections = FALSE,
  # useClusters = TRUE,
  useClusters = FALSE,
  expandPops = FALSE,
  # numClusters = 10,
  numClusters = 5,
  uIDs = regionUIDs
)

task <- cciaObj$runTask(
  funName = "clustRegions.kmeansClust",
  funParams = funParams,
  env = "local",
  runInplace = TRUE,
  taskID = 1
)
```

```{r}
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  # pID = pID, uID = "diRenc", versionID = versionID, initReactivity = FALSE # Set
  pID = pID, uID = "akAcC1", versionID = versionID, initReactivity = FALSE # Set
  # pID = pID, uID = "pN5QJi", versionID = versionID, initReactivity = FALSE # Set
)

regionUIDs <- names(cciaObj$cciaObjects())
# regionUIDs <- regionUIDs[!regionUIDs %in% c(
#   # "84kEUg", "LQFVZD", "6zeSiS"
#   "84kEUg"
# )]

# get experimental info
exp.info <- as.data.table(cciaObj$summary(
  withSelf = FALSE, fields = c("Attr")
))

# regionUIDs <- names(cciaObj$cciaObjects())
# regionUIDs <- exp.info[Include == "Y"]$uID
# regionUIDs <- exp.info[Stain == "A"]$uID
# regionUIDs <- exp.info[Include == "Y" & Stain == "B"]$uID
regionUIDs <- exp.info[Include == "Y" & Stain == "A"]$uID

# get objects from selected set
regionsDT <- cciaObj$popDT(
  popType = "region",
  uIDs = regionUIDs,
  includeFiltered = TRUE,
  completeDT = TRUE, replaceNA = TRUE,
  filterMeasures = c("regions")
)
```

```{r fig_regions_im, fig.height=15, fig.width=40}
# plot out regions
p1 <- ggplot(regionsDT,
             aes(centroid_x, -centroid_y, color = as.factor(regions))) +
  theme_classic() +
  plotThemeDark(angle = 0) +
  geom_point(size = 0.2) +
  scale_color_brewer(palette = "Set1", name = "Region") +
  # scale_color_brewer(palette = "Accent", name = "Region") +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    # legend.position = "none"
    legend.position = "bottom"
  ) + coord_fixed() +
  guides(colour = guide_legend(override.aes = list(size=10))) +
  facet_wrap(.~uID, ncol = 4)

ggsave(file.path(anaDir, "maps", paste0("regions_SPLEEN.png")), width = 24, height = 12, plot = p1)
```

```{r fig_regions_heat, fig.height=5, fig.width=8}
colsToShow <- colnames(regionsDT)
colsToShow <- colsToShow[!colsToShow %in% c("uID", "label", "regions", "pop", "centroid_x", "centroid_y")]
colsToPlot <- regionsDT %>%
  pivot_longer(
    cols = all_of(colsToShow), names_to = "col", values_to = "value"
    )
colsToPlot$col <- as.factor(colsToPlot$col)
colsToPlot$value <- as.numeric(colsToPlot$value)

# show heatmap for clusters
colsSummary <- colsToPlot %>%
  drop_na() %>%
  group_by(col, regions) %>%
  # group_by(regions, col) %>%
  summarise(mean = mean(value, rm.na = TRUE)) %>%
  mutate(norm = (mean - min(mean))/(max(mean) - min(mean))) %>%
  replace_na(list(norm = 0))

# find order for clusters and markers
colsMatrix <- colsSummary %>%
  pivot_wider(names_from = col, values_from = norm, id_cols = regions)

hcClusters <- hclust(dist(as.matrix(colsMatrix %>% dplyr::select(-all_of("regions")))), "ave")
hcMarkers <- hclust(dist(t(as.matrix(colsMatrix %>% dplyr::select(-all_of("regions"))))), "ave")

clustersOrder <- colsMatrix[["regions"]][hcClusters$order]
colsOrder <- colnames(colsMatrix)[2:length(colnames(colsMatrix))][hcMarkers$order]

colsSummary[["regions"]] <- factor(colsSummary[["regions"]], levels = clustersOrder)

p1 <- ggplot(colsSummary %>%
         mutate(
           # clusters = factor(clusters, levels = clustersOrder),
           # clusters = factor(clusters, levels = sort(as.numeric(levels(clustersOrder)))),
           col = factor(col, levels = colsOrder)
         ), aes(col, regions)) +
  theme_classic() +
  # geom_tile(aes(fill = mean), colour = "white", size = 0.5) +
  geom_tile(aes(fill = norm), colour = "white", size = 0.5) +
  viridis::scale_fill_viridis(
    breaks = c(0, 1),
    labels = c(0, 1)
  ) +
  theme(
    legend.position = "none",
    legend.title = element_blank(),
    legend.text = element_text(size = 20),
    legend.key.size = unit(10, "mm"),
    axis.text.y = element_text(size = 20),
    # axis.text.x = element_text(size = 20),
    axis.text.x = element_text(size = 20, angle = 45, hjust = 1, vjust = 1),
    strip.background = element_blank(),
    strip.text.x = element_blank()
    ) +
  xlab("") + ylab("")
plot(p1)

ggsave(file.path(anaDir, "regions_heat.pdf"), width = 6, height = 3)
```

```{r}
# run clustering
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "diRenc", versionID = versionID, initReactivity = FALSE # Set
)

uIDs <- names(cciaObj$cciaObjects())
imChannels <- cciaObj$cciaObjects(uIDs = uIDs)[[1]]$imChannelNames()

# run task
funParams <- list(
  popType = "clust",
  popsToCluster = c(),
  # popsToCluster = c("CD4"),
  savePops = FALSE,
  # savePops = TRUE,
  keepPops = FALSE,
  # keepPops = TRUE,
  resolution = clusterRes$init,
  # resolution = clusterRes$CD4,
  mergeUmap = TRUE,
  normalisePercentile = 99.5,
  # normalisePercentile = 99.0,
  normalisePercentileBottom = 0,
  normaliseAxis = "channels",
  normaliseToMedian = FALSE,
  normaliseIndividually = FALSE,
  # transformation = "log",
  transformation = "NONE",
  logBase = 0,
  valueName = "default",
  uIDs = uIDs,
  clusterChannels = list(
    # base = list(channels = unname(sapply(baseChannels$init, function(x) which(imChannels == x))) - 1),
    base = list(channels = list()),
    nuc = list(channels = unname(sapply(nucChannels$init, function(x) which(imChannels == x))) - 1)
    # base = list(channels = unname(sapply(baseChannels$CD4, function(x) which(imChannels == x))) - 1),
    # nuc = list(channels = unname(sapply(nucChannels$CD4, function(x) which(imChannels == x))) - 1)
  ),
  objectMeasures = objectMeasures$init
  # objectMeasures = objectMeasures$CD4
)

task <- cciaObj$runTask(
  funName = "clustPopulations.leidenClustering",
  funParams = funParams,
  env = "local",
  runInplace = TRUE,
  taskID = 1
)
```


```{r}
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "diRenc", versionID = versionID, initReactivity = FALSE # Set
)

exp.info <- as.data.table(
  cciaObj$summary(withSelf = FALSE, fields = c("Attr")) %>%
    mutate(GroupSeq = as.numeric(GroupSeq)) %>%
    left_join(
      read.csv(file.path(anaDir, "sample.info.csv")),
      by = c("GroupSeq" = "TMA.spot")) %>%
    left_join(read.csv(file.path(anaDir, "patient.info.csv")))
  )
exp.info$Biopsy <- "Post"
exp.info[
  !is.na(stringr::str_match(exp.info$Biopsy.timepoint, "Pre_.*"))[,1], ]$Biopsy <- "Pre"

uIDs <- names(cciaObj$cciaObjects())

pops <- cciaObj$popPaths(uIDs = uIDs[[1]], popType = "clust", includeFiltered = TRUE)
pops <- pops[is.na(stringr::str_match(pops, "/"))]

# get cluster information
popDT <- cciaObj$popDT(
  popType = "clust", includeFiltered = TRUE,
  # pops = c("root", "non.debris", "Th.ae", "Th.malignant", "Treg"), uIDs = uIDs,
  pops = c("root"), completePops = FALSE,
  # pops = pops,
  # uIDs = uIDs,
  completeDT = FALSE, replaceX = FALSE)
  # completeDT = TRUE, replaceX = TRUE)
  # pops = c("Th.ae", "Th.malignant", "Treg"), uIDs = uIDs)

popDT[, clusters := as.numeric(as.character(clusters))]
unique(popDT$clusters)
unique(popDT$pop)
```

```{r}
# define populations
clust <- list(
  init = c("root")
  # CD4 = c("CD4.Tumor", "CD4.ae", "CD4.Treg", "CD4.naive", "CD4.GrB+", "CD4.ICOS+")
  # CD4 = c("CD4")
)

markers <- list(
  init = c(
    # baseChannels$init, paste0("nuc_", nucChannels$init), objectMeasures$init
    paste0("nuc_", nucChannels$init), objectMeasures$init
  )
  # CD4 = c(
  #   baseChannels$CD4, paste0("nuc_", nucChannels$CD4), objectMeasures$CD4
  # )
)

groupVars <- list(
  init = "clusters"
  # CD4 = "clusters"
)

combinePops <- list(
  init = list(
    # CD4 = c("CD4.Tumor", "CD4.ae", "CD4.Treg", "CD4.naive", "CD4.GrB+", "CD4.ICOS+")
  )
  # CD4 = list()
)

plotSizes <- list(
  init = c(12, 20)
  # CD4 = c(4.5, 12)
)
```

```{r fig_clusters_heat, fig.height=4, fig.width=4}
for (i in names(clust)) {
  popsToShow <- clust[[i]]
  markersToShow <- .flowCorrectChannelNames(markers[[i]])
  groupVar <- groupVars[[i]]
  popsToCombine <- combinePops[[i]]
  plotSize <- plotSizes[[i]]
  
  markersToShowNoDye <- stringr::str_extract(markersToShow, ".*(?=_)")
  
  colsToShow <- c("label", "clusters", "pop", markersToShow)
  colsToShow <- colsToShow[colsToShow %in% colnames(popDT)]
  
  # show heatmap
  popLymphs <- popDT[pop %in% popsToShow, ..colsToShow]
  # popLymphs <- popDT[, ..colsToShow]
  
  # combine pops
  if (length(popsToCombine) > 0) {
    for (j in names(popsToCombine))
      popLymphs[pop %in% popsToCombine[[j]], pop := j]
  }
  
  markersToPlot <- popLymphs %>%
    pivot_longer(
      cols = all_of(markersToShow), names_to = "marker", values_to = "value"
      )
  markersToPlot$marker <- as.factor(markersToPlot$marker)
  
  # show heatmap for clusters
  markersSummary <- markersToPlot %>%
    drop_na() %>%
    group_by(marker, get(groupVar)) %>%
    summarise(mean = mean(value, rm.na = TRUE)) %>%
    mutate(norm = (mean - min(mean))/(max(mean) - min(mean))) %>%
    replace_na(list(norm = 0)) %>%
    rename_with(~ c(groupVar), all_of(c("get(groupVar)")))
  
  # find order for clusters and markers
  markersMatrix <- markersSummary %>%
    pivot_wider(names_from = marker, values_from = norm, id_cols = groupVar)
  
  hcClusters <- hclust(dist(as.matrix(markersMatrix %>% dplyr::select(-all_of(groupVar)))), "ave")
  hcMarkers <- hclust(dist(t(as.matrix(markersMatrix %>% dplyr::select(-all_of(groupVar))))), "ave")
  
  clustersOrder <- markersMatrix[[groupVar]][hcClusters$order]
  markerOrder <- colnames(markersMatrix)[2:length(colnames(markersMatrix))][hcMarkers$order]
  
  markersSummary[[groupVar]] <- factor(markersSummary[[groupVar]], levels = clustersOrder)
  
  p1 <- ggplot(markersSummary %>%
           mutate(
             # clusters = factor(clusters, levels = clustersOrder),
             # clusters = factor(clusters, levels = sort(as.numeric(levels(clustersOrder)))),
             marker = factor(marker, levels = markerOrder)
           ), aes(get(groupVar), marker)) +
    theme_classic() +
    # geom_tile(aes(fill = mean), colour = "white", size = 0.5) +
    geom_tile(aes(fill = norm), colour = "white", size = 0.5) +
    viridis::scale_fill_viridis(
      breaks = c(0, 1),
      labels = c(0, 1)
    ) +
    theme(
      legend.position = "none",
      legend.title = element_blank(),
      legend.text = element_text(size = 20),
      legend.key.size = unit(10, "mm"),
      axis.text.y = element_text(size = 20),
      # axis.text.x = element_text(size = 20),
      axis.text.x = element_text(size = 20, angle = 45, hjust = 1, vjust = 1),
      strip.background = element_blank(),
      strip.text.x = element_blank()
      ) +
    xlab("") + ylab("")
  plot(p1)
  
  ggsave(file.path(anaDir, paste0("clusters_heat_", i, ".pdf")),
         width = plotSize[[1]], height = plotSize[[2]])
}
```

```{r fig_pops_umap, fig.height=5, fig.width=5}
library(tidyverse)

dfToPlot <- popDT %>%
  mutate(clusters = as.character(clusters)) %>%
  mutate(clusters = as.numeric(clusters)) %>%
  # dplyr::filter(clusters >= 14) %>%
  left_join(exp.info)

# get mean positions of clusters
meanClusterPos <- dfToPlot %>%
  group_by(clusters) %>%
  # group_by(pop) %>%
  summarise(
    UMAP_1 = mean(UMAP_1),
    UMAP_2 = mean(UMAP_2)
  )

clustersPalette <- randomcoloR::distinctColorPalette(
  length(unique(dfToPlot$clusters)))

# plot UMAP with paga connections
ggplot(
  dfToPlot,
    # mutate(clusters = factor(clusters, levels = clustersOrder)),
  aes(UMAP_1, UMAP_2)
  ) +
  theme_classic() +
  geom_point(aes(color = as.factor(clusters)), size = 0.5) +
  # geom_point(aes(color = as.factor(pop)), size = 0.5) +
  # geom_point(aes(color = interaction(Biopsy, Response)), size = 0.5) +
  # scale_color_brewer(name = NULL, palette = "Set2") +
  scale_color_manual(values = clustersPalette) +
  theme(
    axis.text = element_text(size = 15),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    legend.justification = "center",
    legend.direction = "horizontal",
    # legend.position = "bottom",
    legend.position = "none",
    legend.title = element_blank(),
    legend.key.size = unit(5, "mm"),
    legend.text = element_text(size = 18),
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
    ) +
  guides(colour = guide_legend(override.aes = list(size = 4))) +
  geom_label(
    data = meanClusterPos,
    aes(label = clusters),
    # aes(label = pop),
    label.size = 0.25,
    color = "black"
    )

ggsave(file.path(anaDir, "T_umap.tiff"), width = 5, height = 5)
# ggsave(file.path(anaDir, "T_umap_wo_labels.tiff"), width = 5, height = 5)
```

```{r fig_umap_props, fig.height=4, fig.width=4}
datToPlot <- popDT %>%
  # dplyr::filter(pop == "CD4") %>%
  # dplyr::filter(pop %in% clust$CD4) %>%
  left_join(exp.info)

channels <- c("nuc_FOXP3", "nuc_Ki.67", "ICOS", "CD45RA",
              "CD45RO", "PD.1", "GranzymeB", "nuc_area")

for (i in channels) {
  # plot UMAP
  p1 <- ggplot(datToPlot, aes(UMAP_1, UMAP_2)) +
    theme_classic() +
    geom_point(aes(color = get(i)), size = 0.5, alpha = 1) +
    viridis::scale_color_viridis() +
    theme(
      axis.text = element_text(size = 15),
      axis.text.x = element_blank(),
      axis.ticks.x = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks.y = element_blank(),
      legend.justification = "right",
      # legend.position = "bottom"
      legend.position = "none",
      axis.line = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank()
      )
  
  print(p1)
  ggsave(file.path(anaDir, paste0("umap_", i ,".tiff")), width = 2, height = 2)
}
```

```{r fig_umap, fig.height=4, fig.width=4}
plot.p1 <- function(df.all, df.subset) {
  plot.colours <- .flowColours(df.subset$UMAP_1, df.subset$UMAP_2, colramp = viridisLite::inferno)
  
  ggplot(df.all, aes(UMAP_1, UMAP_2)) +
    theme_classic() +
    geom_point(color = "black", size = 0.5) +
    # geom_point(color = "#B3BCC2", size = 0.5) +
    geom_point(data = df.subset, color = plot.colours, size = 0.5, alpha = 1) +
    # geom_point(data = df.subset, color = "#AA1F5E", size = 0.5, alpha = 1) +
    theme(
      axis.text = element_text(size = 15),
      axis.text.x = element_blank(),
      axis.ticks.x = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks.y = element_blank(),
      legend.justification = "right",
      # legend.position = "bottom"
      legend.position = "none",
      axis.line = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank()
      )
}

datToPlot <- popDT %>%
  # dplyr::filter(pop %in% clust$CD4) %>%
  left_join(exp.info)

# plot UMAP
for (i in unique(unique(datToPlot$Biopsy))) {
  for (j in unique(unique(datToPlot$Response))) {
    plot.p1(datToPlot, datToPlot %>% dplyr::filter(Biopsy == i, Response == j))
    ggsave(file.path(anaDir, paste0("umap_", i, "_", j, ".tiff")), width = 2, height = 2)
  }
}
```

```{r fig_pops_freq, fig.height=5, fig.width=15}
# show frequency of populations for treatments
summaryDF <- popDT %>%
  # dplyr::filter(pop == "CD4") %>%
  # dplyr::filter(!clusters %in% c(0,6)) %>%
  group_by(uID, clusters) %>%
  # group_by(uID, pop) %>%
  summarise(n = n()) %>%
  mutate(freq = n/sum(n) * 100) %>%
  ungroup() %>%
  complete(uID, clusters, fill = list(freq = 0)) %>%
  # complete(uID, pop, fill = list(freq = 0)) %>%
  left_join(exp.info)

# plot
ggplot(summaryDF,
       aes(interaction(Biopsy, Response), freq, fill = interaction(Biopsy, Response))) +
  theme_classic() +
  geom_boxplot(outlier.alpha = 0) +
  geom_jitter(
    position = position_jitterdodge(jitter.width = 0.10), alpha = 1.0) +
    # width = 0.3, alpha = 0.6) +
  # stat_summary(fun=mean, geom="point", size=3, shape=18, color="Magenta") +
  theme(
    axis.text.x = element_text(size = 12, angle = 45, hjust = 1, vjust = 1),
  ) +
  ylab("Population (%)") + xlab("") +
  facet_wrap(.~clusters, ncol = 10, scales = "free_y") +
  # facet_wrap(.~clusters, ncol = 10) +
  # facet_wrap(.~pop, ncol = 10, scales = "free_y") +
  scale_fill_brewer(name = NULL, palette = "Set2")

summaryDF$Response.ID <- 0
summaryDF[summaryDF$Biopsy == "Post" & summaryDF$Response == "n",]$Response.ID <- 1
summaryDF[summaryDF$Biopsy == "Pre" & summaryDF$Response == "y",]$Response.ID <- 2
summaryDF[summaryDF$Biopsy == "Post" & summaryDF$Response == "y",]$Response.ID <- 3

# pops <- unique(summaryDF$pop)
# summaryDF$pop.ID <- unname(sapply(summaryDF$pop, function(x) which(pops == x)))

# save as csv for plotting in Prism
write.csv(summaryDF, file.path(anaDir, "clusters_freq.csv"))
# write.csv(summaryDF, file.path(anaDir, "pop_freq.csv"))
```


```{r fig_pops_freq, fig.height=4, fig.width=12}
# show frequency of regions for treatments
summaryDF <- regionsDT %>%
  # dplyr::filter(regions != 3) %>%
  group_by(uID, regions) %>%
  # group_by(Patient.ID, Biopsy, Response, regions) %>%
  summarise(n = n()) %>%
  mutate(freq = n/sum(n) * 100) %>%
  ungroup() %>%
  complete(uID, regions, fill = list(freq = 0)) %>%
  # complete(Patient.ID, Biopsy, Response, regions, fill = list(freq = 0))
  left_join(exp.info)

# plot
ggplot(summaryDF,
       aes(interaction(Biopsy, Response), freq, fill = interaction(Biopsy, Response))) +
  theme_classic() +
  geom_boxplot(outlier.alpha = 0) +
  geom_jitter(
    position = position_jitterdodge(jitter.width = 0.10), alpha = 1.0) +
    # width = 0.3, alpha = 0.6) +
  # stat_summary(fun=mean, geom="point", size=3, shape=18, color="Magenta") +
  theme(
    axis.text.x = element_text(size = 12, angle = 45, hjust = 1, vjust = 1),
  ) +
  ylab("Region (%)") + xlab("") +
  facet_wrap(.~regions, ncol = 10, scales = "free_y") + ylim(0,100) +
  scale_fill_brewer(name = NULL, palette = "Set3")

summaryDF$Response.ID <- 0
summaryDF[summaryDF$Biopsy == "Post" & summaryDF$Response == "n",]$Response.ID <- 1
summaryDF[summaryDF$Biopsy == "Pre" & summaryDF$Response == "y",]$Response.ID <- 2
summaryDF[summaryDF$Biopsy == "Post" & summaryDF$Response == "y",]$Response.ID <- 3

# save as csv for plotting in Prism
write.csv(summaryDF, file.path(anaDir, "regions_freq.csv"))
```

```{r fig_regions_heat, fig.height=5, fig.width=8}
colsToShow <- colnames(regionsDT)
colsToShow <- colsToShow[!colsToShow %in% c("uID", "label", "regions", "pop")]
colsToPlot <- regionsDT %>%
  pivot_longer(
    cols = all_of(colsToShow), names_to = "col", values_to = "value"
    )
colsToPlot$col <- as.factor(colsToPlot$col)
colsToPlot$value <- as.numeric(colsToPlot$value)

# show heatmap for clusters
colsSummary <- colsToPlot %>%
  drop_na() %>%
  group_by(col, regions) %>%
  # group_by(regions, col) %>%
  summarise(mean = mean(value, rm.na = TRUE)) %>%
  mutate(norm = (mean - min(mean))/(max(mean) - min(mean))) %>%
  replace_na(list(norm = 0))

# find order for clusters and markers
colsMatrix <- colsSummary %>%
  pivot_wider(names_from = col, values_from = norm, id_cols = regions)

hcClusters <- hclust(dist(as.matrix(colsMatrix %>% dplyr::select(-all_of("regions")))), "ave")
hcMarkers <- hclust(dist(t(as.matrix(colsMatrix %>% dplyr::select(-all_of("regions"))))), "ave")

clustersOrder <- colsMatrix[["regions"]][hcClusters$order]
colsOrder <- colnames(colsMatrix)[2:length(colnames(colsMatrix))][hcMarkers$order]

colsSummary[["regions"]] <- factor(colsSummary[["regions"]], levels = clustersOrder)

p1 <- ggplot(colsSummary %>%
         mutate(
           # clusters = factor(clusters, levels = clustersOrder),
           # clusters = factor(clusters, levels = sort(as.numeric(levels(clustersOrder)))),
           col = factor(col, levels = colsOrder)
         ), aes(col, regions)) +
  theme_classic() +
  # geom_tile(aes(fill = mean), colour = "white", size = 0.5) +
  geom_tile(aes(fill = norm), colour = "white", size = 0.5) +
  viridis::scale_fill_viridis(
    breaks = c(0, 1),
    labels = c(0, 1)
  ) +
  theme(
    legend.position = "none",
    legend.title = element_blank(),
    legend.text = element_text(size = 20),
    legend.key.size = unit(10, "mm"),
    axis.text.y = element_text(size = 20),
    # axis.text.x = element_text(size = 20),
    axis.text.x = element_text(size = 20, angle = 45, hjust = 1, vjust = 1),
    strip.background = element_blank(),
    strip.text.x = element_blank()
    ) +
  xlab("") + ylab("")
plot(p1)

ggsave(file.path(anaDir, "regions_heat.pdf"), width = 8, height = 6)
```

```{r}
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "diRenc", versionID = versionID, initReactivity = FALSE # Set
)

uIDs <- names(cciaObj$cciaObjects())

# get cluster information
popDT <- cciaObj$popDT(
  popType = "clust", includeFiltered = TRUE,
  pops = c("Treg", "Th", "malignant", "Treg"), uIDs = uIDs)

# rename inverted
# popDT[is.na(`clust.cell.min_distance#NOT_clust.Th.malignant`),
#       `clust.cell.min_distance#NOT_clust.Th.malignant` := `NOT_clust.cell.min_distance#clust.Th.malignant`]
# popDT[is.na(`clust.cell.min_distance#NOT_clust.Treg`),
#       `clust.cell.min_distance#NOT_clust.Treg` := `NOT_clust.cell.min_distance#clust.Treg`]
# popDT[is.na(`clust.cell.contact#NOT_clust.Th.malignant`),
#       `clust.cell.contact#NOT_clust.Th.malignant` := `NOT_clust.cell.contact#clust.Th.malignant`]
# popDT[is.na(`clust.cell.contact#NOT_clust.Treg`),
#       `clust.cell.contact#NOT_clust.Treg` := `NOT_clust.cell.contact#clust.Treg`]
```

```{r fig_pop_contact, fig.height=4, fig.width=8}
# show frequency of regions for treatments
summaryDF <- popDT %>%
  dplyr::filter(pop == "Th") %>%
  # group_by(uID, `clust.cell.contact#clust.malignant`) %>%
  group_by(uID, `clust.cell.contact#clust.Treg`) %>%
  summarise(n = n()) %>%
  mutate(freq = n/sum(n) * 100) %>%
  # dplyr::filter(`clust.cell.contact#clust.malignant` == FALSE) %>%
  dplyr::filter(`clust.cell.contact#clust.Treg` == FALSE) %>%
  mutate(freq = 100 - freq) %>%
  left_join(exp.info)

# plot
ggplot(summaryDF,
       aes(interaction(Biopsy, Response), freq, fill = interaction(Biopsy, Response))) +
  theme_classic() +
  geom_boxplot(outlier.alpha = 0) +
  geom_jitter(
    position = position_jitterdodge(jitter.width = 0.10), alpha = 1.0) +
    # width = 0.3, alpha = 0.6) +
  # stat_summary(fun=mean, geom="point", size=3, shape=18, color="Magenta") +
  theme(
    axis.text.x = element_text(size = 12, angle = 45, hjust = 1, vjust = 1),
  ) +
  ylab("Contact (%)") + xlab("") +
  scale_fill_brewer(name = NULL, palette = "Set3")

# summaryDF$Response.ID <- 0
# summaryDF[summaryDF$Biopsy == "Post" & summaryDF$Response == "n",]$Response.ID <- 1
# summaryDF[summaryDF$Biopsy == "Pre" & summaryDF$Response == "y",]$Response.ID <- 2
# summaryDF[summaryDF$Biopsy == "Post" & summaryDF$Response == "y",]$Response.ID <- 3

# save as csv for plotting in Prism
# write.csv(summaryDF, file.path(anaDir, "clusters_freq.csv"))
```

```{r fig_pop_contact, fig.height=1.5, fig.width=4}
dtToPlot <- copy(popDT)

# sample root cells
popList <- list()
for (i in uIDs) {
  n.pop <- nrow(dtToPlot[uID == i & pop == "Th.ae"])
  n.root <- nrow(dtToPlot[uID == i & pop == "root"])
  
  n.pop <- if (n.pop < n.root) n.pop else n.root
  
  for (j in seq(10)) {
    n.labels <- sample(dtToPlot[uID == i & pop == "root"]$label, n.pop)
    popList[[paste0(i, ".", j)]] <- dtToPlot[uID == i & label %in% n.labels]
    popList[[paste0(i, ".", j)]][, pop := "inv"]
  }
}

popList[["Th.ae"]] <- dtToPlot[pop == "Th.ae"]

datToPlot <- rbindlist(popList) %>%
  mutate(
    Treg.contact = `clust.cell.contact#clust.Treg`,
    Th.malignant.contact = `clust.cell.contact#clust.Th.malignant`,
    Treg.contact.inv = `clust.cell.contact#NOT_clust.Treg`,
    Th.malignant.contact.inv = `clust.cell.contact#NOT_clust.Th.malignant`
  ) %>%
  pivot_longer(
    cols = c(Treg.contact, Th.malignant.contact, Treg.contact.inv, Th.malignant.contact.inv),
    names_to = "pop.contact",
    values_to = "contact"
  ) %>%
  group_by(uID, pop, pop.contact, contact) %>%
  summarise(n = n()) %>%
  mutate(freq = n/sum(n) * 100) %>%
  dplyr::filter(contact == TRUE) %>%
  left_join(exp.info)

# plot out Th.ae v Treg / Th.malignant
ggplot(datToPlot, aes(
  pop.contact, freq, fill = interaction(Biopsy, Response))) +
  theme_classic() +
  geom_boxplot(outlier.alpha = 0) +
  geom_jitter(
    position = position_jitterdodge(jitter.width = 0.10), alpha = 1.0) +
    # width = 0.3, alpha = 0.6) +
  ylab("Contact (%)") + xlab("") +
  scale_fill_brewer(name = NULL, palette = "Set3")

datToPlot$Response.ID <- 0
datToPlot[datToPlot$Biopsy == "Post" & datToPlot$Response == "n",]$Response.ID <- 1
datToPlot[datToPlot$Biopsy == "Pre" & datToPlot$Response == "y",]$Response.ID <- 2
datToPlot[datToPlot$Biopsy == "Post" & datToPlot$Response == "y",]$Response.ID <- 3

datToPlot$pop.ID <- 0
datToPlot[datToPlot$pop.contact == "Th.malignant.contact", ]$pop.ID <- 1

write.csv(datToPlot, file.path(anaDir, "Th_ae_contact.csv"))
```

```{r fig_pop_dist, fig.height=1.5, fig.width=4}
# sample root cells
popList <- list()
for (i in uIDs) {
  n.pop <- nrow(dtToPlot[uID == i & pop == "Th.ae"])
  n.root <- nrow(dtToPlot[uID == i & pop == "root"])
  
  n.pop <- if (n.pop < n.root) n.pop else n.root
  
  for (j in seq(10)) {
    n.labels <- sample(dtToPlot[uID == i & pop == "root"]$label, n.pop)
    popList[[paste0(i, ".", j)]] <- dtToPlot[uID == i & label %in% n.labels]
    popList[[paste0(i, ".", j)]][, pop := "inv"]
  }
}

popList[["Th.ae"]] <- dtToPlot[pop == "Th.ae"]

datToPlot <- rbindlist(popList) %>%
  mutate(
    Treg.dist = `clust.cell.min_distance#clust.Treg`,
    Th.malignant.dist = `clust.cell.min_distance#clust.Th.malignant`,
    Treg.dist.inv = `clust.cell.min_distance#NOT_clust.Treg`,
    Th.malignant.dist.inv = `clust.cell.min_distance#NOT_clust.Th.malignant`
  ) %>%
  pivot_longer(
    cols = c(Treg.dist, Th.malignant.dist, Treg.dist.inv, Th.malignant.dist.inv),
    names_to = "pop.dist",
    values_to = "dist"
  ) %>%
  group_by(uID, pop.dist) %>%
  summarise(median.dist = median(dist)) %>%
  left_join(exp.info)

# plot out Th.ae v Treg / Th.malignant
ggplot(datToPlot, aes(
  pop.dist, median.dist, fill = interaction(Biopsy, Response))) +
  theme_classic() +
  geom_boxplot(outlier.alpha = 0) +
  geom_jitter(
    position = position_jitterdodge(jitter.width = 0.10), alpha = 1.0) +
    # width = 0.3, alpha = 0.6) +
  ylab("Distance (um)") + xlab("") +
  scale_fill_brewer(name = NULL, palette = "Set3")

datToPlot$Response.ID <- 0
datToPlot[datToPlot$Biopsy == "Post" & datToPlot$Response == "n",]$Response.ID <- 1
datToPlot[datToPlot$Biopsy == "Pre" & datToPlot$Response == "y",]$Response.ID <- 2
datToPlot[datToPlot$Biopsy == "Post" & datToPlot$Response == "y",]$Response.ID <- 3

datToPlot$pop.ID <- 0
datToPlot[datToPlot$pop.dist == "Th.malignant.dist", ]$pop.ID <- 1

write.csv(datToPlot, file.path(anaDir, "Th_ae_dist.csv"))
```

```{r fig_pop_dist_ratio, fig.height=1.5, fig.width=4}
datToPlot <- popDT %>%
  # normalise distances
  mutate(
    Treg.dist = `clust.cell.min_distance#clust.Treg`,
    Th.malignant.dist = `clust.cell.min_distance#clust.Th.malignant`
  ) %>%
  mutate(
    Treg.v.Th.malignant = Treg.dist/Th.malignant.dist,
    Th.malignant.v.Treg = Th.malignant.dist/Treg.dist
  ) %>%
  pivot_longer(
    cols = c(Th.malignant.v.Treg, Treg.v.Th.malignant),
    names_to = "pop.ratio",
    values_to = "ratio"
  ) %>%
  group_by(uID, pop.ratio) %>%
  summarise(median.ratio = median(ratio)) %>%
  left_join(exp.info)

# plot out Th.ae v Treg / Th.malignant
ggplot(datToPlot, aes(
  pop.ratio, median.ratio, fill = interaction(Biopsy, Response))) +
  theme_classic() +
  geom_boxplot(outlier.alpha = 0) +
  geom_jitter(
    position = position_jitterdodge(jitter.width = 0.10), alpha = 1.0) +
    # width = 0.3, alpha = 0.6) +
  ylab("Distance Ratio") + xlab("") +
  scale_fill_brewer(name = NULL, palette = "Set3") +
  facet_wrap(.~pop.ratio, scales = "free")

datToPlot$Response.ID <- 0
datToPlot[datToPlot$Biopsy == "Post" & datToPlot$Response == "n",]$Response.ID <- 1
datToPlot[datToPlot$Biopsy == "Pre" & datToPlot$Response == "y",]$Response.ID <- 2
datToPlot[datToPlot$Biopsy == "Post" & datToPlot$Response == "y",]$Response.ID <- 3

datToPlot$pop.ID <- 0
datToPlot[datToPlot$pop.ratio == "Th.malignant.v.Treg", ]$pop.ID <- 1

write.csv(datToPlot, file.path(anaDir, "Th_ae_dist_ratio.csv"))
```

```{r}
devtools::load_all("../")
cciaUse("~/cecelia/dev", initConda = FALSE)

# init ccia object
cciaObj <- initCciaObject(
  pID = pID, uID = "diRenc", versionID = versionID, initReactivity = FALSE # Set
)

uIDs <- names(cciaObj$cciaObjects())

# get populations
popDT <- cciaObj$popDT(
  popType = "clust", includeFiltered = TRUE,
  pops = cciaObj$popPaths(popType = "clust", includeFiltered = TRUE), uIDs = uIDs)
  # pops = c("Th.ae", "Th.malignant", "Treg"), uIDs = uIDs)

# get cluster information
spatialDT <- cciaObj$spatialDT(valueName = "default", uIDs = uIDs)

# join pops
spatialDT[popDT[, c("uID", "label", "pop")],
          on = c("uID", "to" = "label"),
          pop.to := pop]
spatialDT[popDT[, c("uID", "label", "pop")],
          on = c("uID", "from" = "label"),
          pop.from := pop]

# filter same type associations
spatialDT <- spatialDT[pop.to != pop.from]
```

```{r fig_kmeans_clusters, fig.height=1.5, fig.width=4}
# get "i-niches" and determine a given number of defined clusters
freqRegions <- spatialDT %>%
  # group_by(uID, from, pop.from, pop.to) %>%
  group_by(uID, pop.from, pop.to) %>%
  # group_by(uID, to, pop.to, pop.from) %>%
  summarise(n = n()) %>%
  mutate(freq = n/sum(n) * 100) %>%
  drop_na() %>%
  ungroup() %>%
  complete(uID, pop.from, pop.to, fill = list(freq = 0)) %>%
  left_join(exp.info)
  # pivot_wider(
  #   id_cols = c(uID, pop.from),
  #   names_from = pop.to,
  #   values_from = freq,
  #   values_fill = 0
  # )

# show interactions of individuals
for (i in unique(freqRegions$pop.from)) {
  x <- freqRegions %>% dplyr::filter(pop.from == i) 
  
  p1 <- ggplot(x,
               aes(pop.to, freq, fill = interaction(Biopsy, Response))) +
    theme_classic() +
    geom_boxplot(outlier.alpha = 0) +
    geom_jitter(
      position = position_jitterdodge(jitter.width = 0.10), alpha = 1.0) +
      # width = 0.3, alpha = 0.6) +
    # stat_summary(fun=mean, geom="point", size=3, shape=18, color="Magenta") +
    ylab("Population (%)") + xlab("") +
    # ylim(0, 1) +
    scale_fill_brewer(name = NULL, palette = "Set3") +
    ggtitle(i)

  print(p1)
  
  x$Response.ID <- 0
  x[x$Biopsy == "Post" & x$Response == "n",]$Response.ID <- 1
  x[x$Biopsy == "Pre" & x$Response == "y",]$Response.ID <- 2
  x[x$Biopsy == "Post" & x$Response == "y",]$Response.ID <- 3
  
  x$pop.to.ID <- 0
  x[x$pop.to == "Th.malignant", ]$pop.to.ID <- 1
  x[x$pop.to == "Treg", ]$pop.to.ID <- 2
  
  # save as csv for plotting in Prism
  write.csv(x, file.path(anaDir, "interactions", paste0("interactions_freq_", i, ".csv")))
}
```