---

title: "Neighborhood analysis"

author: Tobias Roider

date: "Last compiled on `r format(Sys.time(), '%d %B, %Y, %X')`"

output: 

  rmdformats::readthedown:

editor_options: 

  chunk_output_type: console

---

```{r options, include=FALSE, warning = FALSE}

library(knitr)

opts_chunk$set(echo=TRUE, tidy=FALSE, include=TRUE, message=FALSE,

               dpi = 100, cache = FALSE, warning = FALSE)

opts_knit$set(root.dir = "../")

options(bitmapType = "cairo")

```

# Load packages, functions and data

```{r Load packages and functions}

library(Seurat)

library(tidyverse)

library(readxl)

source("R/ReadPackages.R")

source("R/Functions.R")

# Read CODEX data (only meta data) 

# Available at BioStudies database (https://www.ebi.ac.uk/biostudies/) under accession number S-BIAD565

codex_annotation <- data.table::fread("data/cells_annotation.csv") %>% tibble() %>% 

  filter(Merged_final!="na")

```

# Function

```{r function}

# Run nearest neighbor analysis

run_NNanalysis <- function(data=NULL, plan_session="sequential", n_workers=1, nn=25, regions=NULL, 

                           add.prop=TRUE){

  require(future)

  require(future.apply)

  require(FNN)

  if(plan_session=="sequential"){

    plan(strategy = plan_session)

  }

  if(plan_session=="multisession"){

    plan(strategy = plan_session, workers=n_workers)

  }

  options(future.globals.maxSize= 8000*1024^2)

  output <- future_lapply(regions, function(r){

    df <- data %>% 

      filter(unique_region==r) 

    df_nn <- get.knn(select(df,x,y), nn)

    df_nn <- df_nn$nn.index %>% 

      `rownames<-`(df$unique_cell_id) %>% 

      data.frame() %>% 

      rownames_to_column("unique_cell_id") %>% 

      pivot_longer(cols = 2:ncol(.)) %>% 

      left_join(., df %>% select(unique_cell_id_nn=unique_cell_id, Merged_final) %>% mutate(value=1:nrow(.)), by="value") %>% 

      select(-value) 

    if(add.prop==TRUE){

      mat <- df_nn %>% 

        add_prop(vars = c("unique_cell_id", "Merged_final"), group.vars = 1) %>% 

        pivot_wider(names_from = "Merged_final", values_from = "Prop", values_fill = 0) %>% 

        column_to_rownames("unique_cell_id")

      mat <- mat[names(which(rowSums(mat)>0)), ]

    } else {

      mat <- df_nn 

      }

    return(mat)

  })

  plan(strategy = "sequential")

  mat <- do.call(rbind, output)

  return(mat)

}

```

# Run Analysis

```{r run analysis}

# Read and process codex data (uses future package)

nn <- run_NNanalysis(data = codex_annotation, regions = unique(codex_annotation$unique_region), 

                     plan_session = "multisession", 

                     add.prop=TRUE,

                     n_workers = 10, 

                     nn = 25)

# Use k-means clustering to determine neighborhood classes

set.seed(1)

hclust <- kmeans(nn, centers = 10, iter.max = 50, nstart = 1)

nn_classes <- data.frame(Region=hclust$cluster)  %>% 

  mutate(Region=as.character(Region)) %>% 

  rownames_to_column("unique_cell_id")

# Use PCA on nearst neighbour matrix

pca_codex <- prcomp(nn, scale. = T, center = T)

# Remove objects

rm(list = setdiff(ls(), c("nn", "nn_classes", "pca_codex", "run_NNanalysis")))

```

# Save

```{r save}

save.image("output/Neighborhood_results.RData")

```

# Session Info

```{r session info}

sessionInfo()

```