add src/IDR.Rmd

src/IDR.Rmd

+---
+title: "Sister cell IDR"
+author: "Laurent Modolo"
+date: '`r Sys.Date()`'
+output: html_document
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+
+# helper funciton to compute an cell * sister, idr value matrix
+compute_matrix <- function(data){
+  data %>%
+    pivot_wider(id_cols = gene, names_from = sister, values_from = idr) %>%
+    (function(x){
+      res <- x %>% select(-gene) %>% as.matrix()
+      rownames(res) <-  x %>% pull(gene)
+      return(res)
+    })
+}
+
+# helper funciton to build an heatmap from the compute_matrix function output
+compute_heatmap <- function(data){
+  data %>%
+    t() %>%
+    log() %>%
+    ggheatmap(
+      legendName = "log(IDR)",
+      cluster_rows = T,
+      cluster_cols = T
+    ) %>% print()
+}
+```
+
+We load the necessary packages
+
+```{r library}
+library(idr) # devtools::install_git("https://gitbio.ens-lyon.fr/LBMC/sbdm/idr")
+library(tidyverse)
+library(parallel)
+library(ggheatmap)
+```
+
+# data loading
+
+```{r load_data}
+getwd()
+data <- read_tsv("../../data/CF47_CS.csv", show_col_types = F) %>%
+  pivot_longer(-c(gene), values_to = "count") %>%
+  mutate(
+    condition = str_replace(name, "^([LD]M\\d+)\\..*", "\\1") %>% as.factor(),
+    cell = str_replace(name, "^[LD]M\\d+\\.(.*)", "\\1") %>%  as.factor(),
+    sister = str_replace(name, "^[LD]M\\d+\\.(.*)[AB]", "\\1") %>%  as.factor()
+  ) %>%
+  select(-name)
+data %>% summary
+```
+
+# Archimedean IDR computation
+
+We compute the for each gene between each pair of cells
+
+```{r idr_arch, cache=T, dependson="load_data"}
+idr_data <- data %>%
+  group_by(sister) %>%
+  nest(data = c(count, gene, cell)) %>%
+  mutate(
+    data_pair = map(data, function(x){
+      x %>%
+        pivot_wider(id_cols = gene, names_from = cell, values_from = count) %>%
+        arrange(gene)
+    }),
+    idr = mclapply(data_pair, function(x){
+      x %>%
+        select(-c(gene)) %>%
+        compute_idr()
+    }, mc.cores = 4)
+  )
+```
+We plot the IDR values for each pairs of cells
+
+```{r idr_arch_plots}
+idr_data %>%
+  mutate(
+    plot = map2(data_pair, idr, function(x, y){
+      x %>% select(-c(gene)) %>%
+      pairs_plot(idr = y, alpha = 0) %>% print()
+    })
+  )
+```
+
+We may want to rewrite the IDR for gene, for which the expression is zero in the pair of cells, to set it to `1`.
+
+```{r idr_arch_zeros}
+idr_data <- idr_data %>%
+  mutate(
+    idr = map2(data_pair, idr, function(x, y){
+      all_zero <- x %>%
+        select(-gene) %>%
+        rowSums()
+      y[all_zero == 0] <- 1
+      return(y)
+    })
+  )
+```
+
+
+We display the IDR heatmap for genes x pairs of cells
+
+```{r idr_arch_heatmap}
+idr_matrix <-
+  idr_data %>%
+  mutate(
+    gene = map(data_pair, function(x){
+      x %>% pull(gene)
+    })
+  ) %>%
+  unnest(c(gene, idr)) %>%
+  group_by(condition) %>%
+  nest() %>%
+  mutate(
+    matrix = map(data, compute_matrix),
+    heatmap = map(matrix, compute_heatmap)
+  )
+```
+
+# Gaussian IDR computation
+
+We compute the for each gene between each pair of cells
+
+```{r idr_gauss, cache=T, dependson="load_data"}
+idr_g_data <- data %>%
+  group_by(sister) %>%
+  nest(data = c(count, gene, cell)) %>%
+  mutate(
+    data_pair = map(data, function(x){
+      x %>%
+        pivot_wider(id_cols = gene, names_from = cell, values_from = count) %>%
+        arrange(gene)
+    }),
+    idr = mclapply(data_pair, function(x){
+      x %>%
+        select(-c(gene)) %>%
+        compute_idr(samic = F)
+    }, mc.cores = 4)
+  )
+```
+
+We plot the IDR values for each pairs of cells
+
+```{r idr_gauss_plots}
+idr_g_data %>%
+  mutate(
+    plot = map2(data_pair, idr, function(x, y){
+      x %>% select(-c(gene)) %>%
+      pairs_plot(idr = y, alpha = 0) %>% print()
+    })
+  )
+```
+
+We may want to rewrite the IDR for gene, for which the expression is zero in the pair of cells, to set it to `1`.
+
+```{r idr_gauss_zeros}
+idr_g_data <- idr_g_data %>%
+  mutate(
+    idr = map2(data_pair, idr, function(x, y){
+      all_zero <- x %>%
+        select(-gene) %>%
+        rowSums()
+      y[all_zero == 0] <- 1
+      return(y)
+    })
+  )
+```
+
+
+We display the IDR heatmap for genes x pairs of cells
+
+```{r idr_gauss_heatmap}
+idr_g_matrix <-
+  idr_g_data %>%
+  mutate(
+    gene = map(data_pair, function(x){
+      x %>% pull(gene)
+    })
+  ) %>%
+  unnest(c(gene, idr)) %>%
+  group_by(condition) %>%
+  nest() %>%
+  mutate(
+    matrix = map(data, compute_matrix),
+    heatmap = map(matrix, compute_heatmap)
+  )
+```