UCell_benchmarks.small.machines.Rmd

---
title: UCell demo
author: 
- Massimo Andreatta^[massimo.andreatta@unil.ch]
- Santiago Carmona^[santiago.carmona@unil.ch]
date: "15/02/2021"
knit: (function(input_file, encoding) {
  out_dir <- 'docs';
  rmarkdown::render(input_file,
 encoding=encoding,
 output_file=file.path(dirname(input_file), out_dir, 'index.html'))})
#output: html_notebook
---


```{r message=F, warning=F}
#renv::restore()
library(ggplot2)
library(plotly)
#remotes::install_github("mojaveazure/seurat-disk")
library(SeuratDisk)
library(Seurat)
#remotes::install_git("https://gitlab.unil.ch/carmona/UCell.git")
remotes::install_github("carmonalab/UCell")

#system("R CMD build ../UCell")
#system("R CMD INSTALL UCell_0.2.2.tar.gz")

library(UCell)
library(BiocManager)
library(AUCell)
```

# Limit amount of available memory (simulate smaller machine)

To limit memory, we can set a memory limit in .Renviron (local for this project) to e.g. R_MAX_VSIZE=8Gb
Then re-run code, to see the effect on performance
#NB! remember to comment out R_MAX_VSIZE after you're done!

Also consider variable R_GC_MEM_GROW:
"The strategy used for growth can be specified by setting the environment variable R_GC_MEM_GROW to an integer value between 0 and 3. This variable is read at start-up. Higher values grow the heap more aggressively, thus reducing garbage collection time but using more memory."

Run gc(): you will see that "gc trigger", the threshold of memory that prompts garbage collection, is equal to the limit in Mb for a small machine (e.g. 8GB). It's tricky to set a good value for R_GC_MEM_GROW because "gc trigger" depends on current load on the system and will change with time...


```{r}
max_ram <- Sys.getenv("R_MAX_VSIZE")
max_ram
Sys.getenv("R_GC_MEM_GROW")
```

# Prepare subsets for benchmark

Generate subsets of different sizes for the benchmark

First, load a large atlas for benchmarking (>100.000 cells)
```{r}
make_sets=F
sizes <- c(100, 200, 500, 1000, 2000, 5000, 1e4, 2e4, 5e4, 1e5)

if (make_sets) {
  pbmc.azimuth.full <- LoadH5Seurat("../Azimuth/data/pbmc_multimodal.h5seurat")
  exp.matrix <- pbmc.azimuth.full@assays$SCT@counts
  
  seurat <- CreateSeuratObject(counts = exp.matrix)
  saveRDS(seurat, file="data/pbmc.seurat.4bench.rds")
  saveRDS(exp.matrix, file="data/pbmc.expmat.rds")

  for (size in sizes) {
    matrix.sub <- sprintf("data/sub/pbmc.expmat.%i.rds",size)
    this <- exp.matrix[,1:size]
    saveRDS(this, matrix.sub)
  }
}
```

Define signatures for human T cell subtypes
```{r, eval=T}
HCA.markers.Hs.Tcell <- list()
HCA.markers.Hs.Tcell$Tcell_CD4 <- c("CD4","CD40LG")
HCA.markers.Hs.Tcell$Tcell_CD8 <- c("CD8A","CD8B")
HCA.markers.Hs.Tcell$Tcell_Treg <- c("FOXP3","IL2RA")
HCA.markers.Hs.Tcell$Tcell_MAIT <- c("KLRB1", "SLC4A10", "NCR3")
HCA.markers.Hs.Tcell$Tcell_gd <- c("TRDC", "TRGC1", "TRGC2", "TRDV1")
HCA.markers.Hs.Tcell$Tcell_NK <- c("FGFBP2", "SPON2", "KLRF1", "FCGR3A", "KLRD1", "TRDC")
```


```{r eval=T}
library(tidyr)

#testSamp <- c(100,1000,2000,3000)
#testSamp <- c(100,1000,2000,5000, 1e4)
testSamp <- c(100, 200, 500, 1000, 2000, 5000, 1e4, 2e4, 5e4, 1e5)
chunk.size <- 1000
force.gc <- FALSE
#gcinfo(TRUE)

features <- HCA.markers.Hs.Tcell
  
time_table <- matrix(NA,nrow = length(testSamp),ncol = 3)
colnames(time_table) <- c("size","AUCell","UCell")

memory_table <- matrix(NA,nrow = length(testSamp),ncol = 3)
colnames(memory_table) <- c("size","AUCell","UCell")
gc()  #check triggers

for (i in seq_along(testSamp)[]){
  size <- testSamp[i]
  time_table[i,1] <- size
  memory_table[i,1] <- size
  
  matrix.sub <- sprintf("data/sub/pbmc.expmat.%i.rds",size)
  this.data <- readRDS(matrix.sub)
  
  print(size)
  
  gc1 <- gc(reset = TRUE)
  t <- system.time({
    
    out <- tryCatch({
       scores_UCell <- ScoreSignatures_UCell(this.data, features = features, chunk.size = chunk.size, force.gc = force.gc)
       1
    },
    error=function(cond) {
      message(cond)
      return(NA)
    })
  })
  gc2 <- gc()
  if (is.na(out)) { #Out of memory
    memPeak <- NA
    time <- NA
  } else {
    memPeak <- sum(gc2[,7]) - sum(gc1[,7])
    time <- t[["elapsed"]]
  }
  
  time_table[i,3] <- time
  memory_table[i,3] <- memPeak
  
  
  gc1 <- gc(reset = TRUE)
  t <- system.time({
    
    out <- tryCatch({
      cells_rankings <- AUCell_buildRankings(this.data, nCores = 1, plotStats = F)
      cells_AUC <- AUCell_calcAUC(features, cells_rankings, aucMaxRank=1000)
      scores_AUC <- as.data.frame(t(getAUC(cells_AUC)))
      1
    },
    error=function(cond) {
      message(cond)
      return(NA)
    })
    
  })
  gc2 <- gc()
  if (is.na(out)) { #Out of memory
    memPeak <- NA
    time <- NA
  } else {
    memPeak <- sum(gc2[,7]) - sum(gc1[,7])
    time <- t[["elapsed"]]
  }
  
  time_table[i,2] <- time
  memory_table[i,2] <- memPeak
  
}
#gcinfo(FALSE)

#Print results
time_table.df <- time_table %>% as.data.frame() %>% pivot_longer(-size,values_to="time") 
memory_table.df <- memory_table %>% as.data.frame() %>% pivot_longer(-size,values_to="memory") 
memory_table.df$memory <- memory_table.df$memory/1000  #Convert to Gb

max_ram_int <- as.numeric(gsub("(\\d+)Gb", "\\1", max_ram))
min_ram_int <- min(memory_table.df$memory, na.rm = T)

colors <- c("#eba223","#23a5eb")
ggplot(time_table.df, aes(x=size, y = time, name)) +  geom_point(aes(color = name)) + 
  scale_x_log10() + scale_y_log10() + xlab("Size (# cells)") + ylab("Time (seconds)") + 
  scale_color_manual(values=colors) + ggtitle(sprintf("Machine with %s RAM", max_ram)) + theme_bw()

ggsave(sprintf("plots/benchmark_time_direct.%iGB.png", max_ram_int), height=4, width=5)


ggplot(memory_table.df, aes(x=size, y = memory, name)) +  geom_point(aes(color = name)) + 
  scale_x_log10() + scale_y_log10(limits = c(min_ram_int,max_ram_int)) + xlab("Size (# cells)") + ylab("Memory (GB)") + 
  scale_color_manual(values=colors) + ggtitle(sprintf("Machine with %s RAM", max_ram)) + theme_bw()

ggsave(sprintf("plots/benchmark_mem_direct.%iGB.png", max_ram_int), height=4, width=5)

```

#Even better, run benchmark externally using Rscript: wrapper_benchmark_UCell.R

Then we load the results here and prepare the plots.

#Machine with 8GB
```{r}
res.8GB.UCell <- read.table(file="bench_results/UCell.8Gb.bench.test.txt", header = T)
res.8GB.AUCell <- read.table(file="bench_results/AUCell.8Gb.bench.test.txt", header = T)
res<- rbind(res.8GB.UCell, res.8GB.AUCell)
res$PeakMemory <- res$PeakMemory/1000   #Convert to GB

max_ram_int <- 8
min_ram_int <- min(res$PeakMemory, na.rm = T)
colors <- c("#eba223","#23a5eb")

p1 <- ggplot(res, aes(x=Size, y = PeakMemory, Method)) +  geom_point(aes(color = Method)) + 
  scale_x_log10() + scale_y_continuous(trans='log2', limits = c(min_ram_int,max_ram_int), breaks=c(1/4, 1, 4, 8)) +
  xlab("Size (# cells)") + ylab("Memory (GB)") + scale_color_manual(values=colors) + ggtitle("Machine with 8GB RAM") + 
  theme_bw() + geom_hline(yintercept=max_ram_int, linetype="dashed", color = "black", size=1)

p2 <- ggplot(res, aes(x=Size, y = Time, Method)) +  geom_point(aes(color = Method)) + 
  scale_x_log10() + scale_y_log10(breaks=c(0.3, 1, 3, 10, 30, 100)) + xlab("Size (# cells)") + ylab("Time (seconds)") + 
  scale_color_manual(values=colors) + ggtitle("Machine with 8GB RAM") + theme_bw()


p1
p2

ggsave(sprintf("plots/benchmark_memory.%iGB.pdf", max_ram_int), plot=p1, height=3, width=5)
ggsave(sprintf("plots/benchmark_time.%iGB.pdf", max_ram_int), plot=p2, height=3, width=5)

```


#Machine with 128GB
```{r}
res.128GB.UCell <- read.table(file="bench_results/UCell.128Gb.bench.test.txt", header = T)
res.128GB.AUCell <- read.table(file="bench_results/AUCell.128Gb.bench.test.txt", header = T)
res <- rbind(res.128GB.UCell, res.128GB.AUCell)
res$PeakMemory <- res$PeakMemory/1000   #Convert to GB

max_ram_int <- 128
min_ram_int <- min(res$PeakMemory, na.rm = T)
colors <- c("#eba223","#23a5eb")

p1 <- ggplot(res, aes(x=Size, y = PeakMemory, Method)) +  geom_point(aes(color = Method)) + 
  scale_x_log10() + scale_y_continuous(trans='log2', limits = c(min_ram_int,max_ram_int), breaks=c(1/4, 1, 4, 16, 64)) +
  xlab("Size (# cells)") + ylab("Peak Memory (GB)") + scale_color_manual(values=colors) + ggtitle("Machine with 128GB RAM") + 
  theme_bw()

p2 <- ggplot(res, aes(x=Size, y = Time, Method)) +  geom_point(aes(color = Method)) + 
  scale_x_log10() + scale_y_log10(breaks=c(0.3, 1, 3, 10, 30, 100, 300, 1000)) + xlab("Size (# cells)") + ylab("Time (seconds)") + 
  scale_color_manual(values=colors) + ggtitle("Machine with 128GB RAM") + theme_bw()

p3 <- ggplot(res, aes(x=Size, y = Time, Method)) +  geom_point(aes(color = Method)) + 
  scale_x_log10() + xlab("Size (# cells)") + ylab("Time (seconds)") + 
  scale_color_manual(values=colors) + ggtitle("Machine with 128GB RAM") + theme_bw()



p1
p2
p3

ggsave(sprintf("plots/benchmark_memory.%iGB.pdf", max_ram_int), plot=p1, height=3, width=4)
ggsave(sprintf("plots/benchmark_time.%iGB.pdf", max_ram_int), plot=p2, height=3, width=4)

ggsave(sprintf("plots/benchmark_memory.%iGB.png", max_ram_int), plot=p1, height=3, width=4)
ggsave(sprintf("plots/benchmark_time.%iGB.png", max_ram_int), plot=p2, height=3, width=4)



```