hca_data_analysis.Rmd

---
title: "Analysis of PBMC + HEK data profiled with various scRNA-seq protocols"
author: "Christian Holland"
date: "24/09/2018"
output: html_document
---
```{r "knitr config", cache=FALSE, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
knitr::knit(..., quiet = TRUE)
```

### Libraries and sources
These libraries and sources are used in this analysis 
```{r "setup", message=F}
library(GSEABase)
library(scran)
library(tidyverse)
library(tidylog)
library(ggrepel)
library(cowplot)
library(furrr)
library(pheatmap)
library(AUCell)
library(Rtsne)
library(cluster)
library(Matrix)
library(umap)
library(reticulate)
library(furrr)
library(viper)
library(scMerge)
library(ggsignif)
library(ggpubr)
library(rstatix)
library(ggdendro)
library(AUCell)
library(msigdf)
library(fgsea)
library(Seurat)
library(corrr)

source("src/my_ggplot_themes.R")
source("src/dorothea_analysis.R")
theme_set(theme_cowplot())
plan(multiprocess, workers = 4)
options("tidylog.display" = list(print))
options(future.globals.maxSize= 2097152000) # corresponds to 2 GB

```

### Data pre-processing
```{r "pre-process data"}
sce = get(load("data/hca_data/expression_data/sce.all.technologies.RData"))

raw = assay(sce) %>%
  data.frame(check.names = F, stringsAsFactors = F)

meta = colData(sce) %>%
  data.frame(check.names = F, stringsAsFactors = F) %>%
  rownames_to_column("cell") %>%
  as_tibble() %>%
  dplyr::rename(technology = batch) %>%
  mutate(technology = as_factor(technology))

raw_expr = meta %>% 
  separate(cell, into = c("id"), remove = F) %>%
  distinct(id, technology) %>%
  select(technology, id) %>%
  nest(-technology, .key="id") %>%
  transmute(technology, raw = id %>% map(function(i) {
    # subset raw matrix to cells belonging to given technology
    sub_raw = raw %>% 
      select(contains(pull(i, id)))
    
    # keep only genes which are expressed in more than 5 cells
    keep_genes = apply(sub_raw,1, function(row) {
      sum(row != 0) >= 5
    })
    
    processed_raw = sub_raw[keep_genes,] %>%
      as.matrix() %>%
      Matrix(sparse=T)
  }))


saveRDS(raw_expr, "output/hca_data/expression/raw_preprocessed.rds")
saveRDS(meta, "output/hca_data/meta/meta_df.rds")
```

### Exploraroty analysis
#### Library sizes
```{r "library sizes"}
count_mat = readRDS("output/hca_data/expression/raw_preprocessed.rds")
df = count_mat %>%
  transmute(technology, x = raw %>% map(function(raw) {
    raw %>% 
      as.matrix() %>% 
      colSums() %>% 
      enframe("cell", "libsize")
  }))

num_cells = df %>% 
  unnest() %>% 
  dplyr::count(technology)

df %>% 
  unnest() %>%
  ggplot(aes(x=fct_reorder(technology, libsize, .fun = median), y=libsize)) +
  geom_jitter(aes(color=technology), alpha=0.5) +
  geom_violin(fill="white", aes(color=technology)) +
  geom_text(data = num_cells, aes(x=technology, y = 20000, label = n)) +
  labs(x="Technology", y= "Library size") +
  theme(legend.position = "none") +
  my_theme() +
  coord_flip()
```

#### Numbers of different celltype per technology
```{r "numbers of different celltype per technology"}
meta = readRDS("output/hca_data/meta/meta_df.rds")

meta %>%
  #filter(technology == "Quartz-Seq2") %>%
  dplyr::count(technology, nnet2) %>%
  ggplot(aes(x=nnet2, y=log2(n), fill=nnet2)) +
  geom_col() +
  coord_flip() +
  facet_wrap(~technology, scales="free_x") +
  theme(legend.position = "none") +
  my_theme() +
  geom_hline(yintercept = log2(50)) +
  labs(y = "log2(count)", x="Celltype")
```

#### Gene coverage per technology
```{r "gene coverage per technology"}
count_mat = readRDS("output/hca_data/expression/raw_preprocessed.rds")
df = count_mat %>%
  transmute(technology, x = raw %>% map(function(raw) {
    raw %>% 
      summary() %>%
      as_tibble() %>%
      dplyr::rename(cell=j) %>%
      dplyr::count(cell)
  }))

num_cells = df %>% 
  unnest() %>%
  distinct(technology, cell) %>%
  dplyr::count(technology)

df %>%
  unnest(x) %>%
  ggplot(aes(x=fct_reorder(technology, n, .fun=median), y=n, color=technology)) +
  geom_boxplot() +
  geom_text(data = num_cells, aes(x=technology, y = 7000, label = n)) +
  labs(x = "Technology", y="Number of non-zero genes") +
  my_theme() +
  theme(legend.position = "none") +
  coord_flip()
```

#### Relationship between number of non-zero genes and number of cells
```{r "relationship between number of non-zero genes and number of cells}
count_mat = readRDS("output/hca_data/expression/raw_preprocessed.rds")
df = count_mat %>%
  transmute(technology, x = raw %>% map(function(raw) {
    raw %>% 
      summary() %>%
      as_tibble() %>%
      dplyr::rename(cell=j) %>%
      dplyr::count(cell, name = "num_non_zero_genes") %>%
      add_tally(name = "num_cells") %>%
      select(-cell)
  }))

df %>%
  unnest() %>%
  group_by(technology) %>%
  summarise(num_non_zero_genes = mean(num_non_zero_genes),
            num_cells = mean(num_cells)) %>%
  ggplot(aes(x=num_cells, y=num_non_zero_genes, label=technology)) +
  geom_point() +
  geom_smooth(method="lm") +
  my_theme() +
  geom_text_repel()
```

### Normalization
#### Individual normalization with srcan
```{r "Individual normalization with scran}
count_mat = readRDS("output/hca_data/expression/raw_preprocessed.rds")

expr_norm = count_mat %>%
  transmute(technology, norm = raw %>% map(function(raw) {
    SingleCellExperiment(list(counts=raw)) %>%
      computeSumFactors() %>%
      normalize() %>%
      exprs()
  }))

saveRDS(expr_norm, "output/hca_data/expression/norm.rds")
```

### Compare DoRothEA with GTEx regulons to find common TFs
```{r}
gtex_regulons = list.files("data/regulons/gtex", full.names = T) %>%
  map(function(regulon_path) {
    tissue_regulon = get(load(regulon_path))
    tf_names = names(tissue_regulon) %>% str_split(pattern = " ", simplify = T)
    names(tissue_regulon) = tf_names[,1]
    return(tissue_regulon)
  }) %>%
  map(names) %>% 
  unlist() %>% 
  unique() %>% 
  enframe(name = NULL, value = "tf")


dorothea = read_csv("data/regulons/dorothea/dorothea_regulon_human_v1.csv") %>%
  distinct(tf)

anti_join(dorothea, gtex_regulons, by="tf") %>%
  saveRDS("output/hca_data/meta/tfs_available_in_dorothea_but_not_gtex.rds")
  
anti_join(gtex_regulons, dorothea, by="tf") %>%
  saveRDS("output/hca_data/meta/tfs_available_in_gtex_but_not_dorothea.rds")
```

### Run Footprint methods
#### Run PROGENy with various footprint number
```{r "run progeny with various footprint number"}
expr_norm = readRDS("output/hca_data/expression/norm.rds") %>%
  mutate(count = 6) %>% 
  uncount(count) %>%
  group_by(technology) %>%
  mutate(footprints = c(100,200,300,500,1000, "all")) %>%
  ungroup() %>%
  dplyr::rename(expr = norm)

model = readRDS("data/progeny_models/full_matrix/full_progeny_matrix.rds")
m_100 = model %>% group_by(pathway) %>% top_n(100, -adj.p) %>% ungroup() %>%
  select(gene, pathway, weight=zscore) %>%
  spread(pathway, weight, fill=0) %>%
  data.frame(row.names = 1, check.names = F, stringsAsFactors = F)

# top 200 
m_200 = model %>% group_by(pathway) %>% top_n(200, -adj.p) %>% ungroup() %>%
  select(gene, pathway, weight=zscore) %>%
  spread(pathway, weight, fill=0) %>%
  data.frame(row.names = 1, check.names = F, stringsAsFactors = F)

# top 300 
m_300 = model %>% group_by(pathway) %>% top_n(300, -adj.p) %>% ungroup() %>%
  select(gene, pathway, weight=zscore) %>%
  spread(pathway, weight, fill=0) %>%
  data.frame(row.names = 1, check.names = F, stringsAsFactors = F)

# top 500 
m_500 = model %>% group_by(pathway) %>% top_n(500, -adj.p) %>% ungroup() %>%
  select(gene, pathway, weight=zscore) %>%
  spread(pathway, weight, fill=0) %>%
  data.frame(row.names = 1, check.names = F, stringsAsFactors = F)

# top 1000 
m_1000 = model %>% group_by(pathway) %>% top_n(1000, -adj.p) %>% ungroup() %>%
  select(gene, pathway, weight=zscore) %>%
  spread(pathway, weight, fill=0) %>%
  data.frame(row.names = 1, check.names = F, stringsAsFactors = F)

# all 
m_all = model %>%
  select(gene, pathway, weight=zscore) %>%
  spread(pathway, weight, fill=0) %>%
  data.frame(row.names = 1, check.names = F, stringsAsFactors = F)

models = tribble(
  ~M, ~footprints,
  m_100, 100,
  m_200, 200,
  m_300, 300,
  m_500, 500,
  m_1000, 1000,
  m_all, "all"
) 

setup = expr_norm %>%
  left_join(models, by="footprints") %>%
  mutate(footprints = as_factor(footprints))

sc_progeny = function(expr, M, ...) {
  # expr = t(scale(t(expr)))
  common_genes = intersect(rownames(expr), rownames(M))
  expr_matched = expr[common_genes,,drop=FALSE] %>%
    t()
  M_matched = M[common_genes,, drop=FALSE] %>%
    data.matrix()
  
  stopifnot(names(expr_matched) == rownames(M_matched))
  
  scores = expr_matched %*% M_matched
  res = scale(scores)
  return(res)
}

progeny_scores = setup %>%
  transmute(technology, footprints, progeny = future_pmap(., .f=sc_progeny, .progress = T))

saveRDS(progeny_scores, "output/hca_data/progeny/progeny_scores.rds")
```

#### Run VIPER
```{r "run viper"}
expr_norm = readRDS("output/hca_data/expression/norm.rds")
missing_tfs = readRDS("output/hca_data/meta/tfs_available_in_dorothea_but_not_gtex.rds")

regulon = read_csv("data/regulons/dorothea/dorothea_regulon_human_v1.csv") %>%
  anti_join(missing_tfs, by="tf") %>%
  df2regulon()

sc_viper = function(expr, regulon, ...) {
  res = viper(eset = as.matrix(expr), regulon = regulon, nes = T, method = "scale",
              minsize = 4, eset.filter = F)
  return(res)
}

viper_scores = expr_norm %>%
  transmute(technology, viper = norm %>% future_map(sc_viper, regulon=regulon, .progress = T))

saveRDS(viper_scores, "output/hca_data/viper/viper_scores.rds")
```

#### Run MetaVIPER
```{r}
# we need to remove this tf as it is not available in dorothea regulons
missing_tfs = readRDS("output/hca_data/meta/tfs_available_in_gtex_but_not_dorothea.rds") %>%
  pull(tf)


expr_norm = readRDS("output/hca_data/expression/norm.rds")
gtex_regulons = list.files("data/regulons/gtex", full.names = T) %>%
  map(function(regulon_path) {
    tissue_regulon = get(load(regulon_path))
    tf_names = names(tissue_regulon) %>% str_split(pattern = " ", simplify = T)
    names(tissue_regulon) = tf_names[,1]
    
    # index of tfs we want to keep
    keep_ix = which(!names(tissue_regulon) %in% missing_tfs)
    sub_tissue_regulon = tissue_regulon[keep_ix]
    return(sub_tissue_regulon)
  })

ss_metaviper_result = expr_norm %>%
  mutate(ss_metaviper_result = norm %>% map(
    ~viper(eset = as.matrix(.x),
           regulon = gtex_regulons,
           nes = T, method = "scale",minsize = 4,eset.filter = F, 
           adaptive.size = F, cores = 4))) %>%
  select(-norm)

saveRDS(ss_metaviper_result, "output/hca_data/metaviper/metaviper_scores.rds")
```

#### Run AUCell - DoRothEA
```{r}
# we need to remove this tf as it is not available in gtex regulons
missing_tfs = readRDS("output/hca_data/meta/tfs_available_in_dorothea_but_not_gtex.rds")

# subset dorothea
human_dorothea = read_csv("data/regulons/dorothea/dorothea_regulon_human_v1.csv") %>%
  # filter(confidence %in% c("A", "B")) %>%
  anti_join(missing_tfs, by="tf") %>%
  distinct(tf, target, confidence) %>%
  add_count(tf) %>%
  filter(n>=4) %>%
  select(-n)

genesets = human_dorothea %>%
  group_by(tf) %>%
  summarise(geneset = list(GeneSet(target))) %>%
  transmute(tf, geneset2 = pmap(., .f=function(tf, geneset, ...) {
    setName(geneset) = tf
    return(geneset)
  })) %>%
  deframe() %>%
  GeneSetCollection()

df = readRDS("output/hca_data/expression/norm.rds")

ss_daucell_scores = df %>%
  mutate(ss_daucell_result = norm %>% map(function(norm) {
    print("hi")
    # g = subsetGeneSets(genesets, rownames(emat)) 
    obj = AUCell_buildRankings(norm, nCores=4, plotStats = F, verbose = F) %>%
      AUCell_calcAUC(genesets, ., verbose=F)
    
    res = obj@assays$data@listData$AUC
  })) %>%
  select(-norm)

saveRDS(ss_daucell_scores, "output/hca_data/daucell/daucell_result.rds")
```

#### Run AUCell - PROGENy
```{r}
# files required for the pipeline
model = readRDS("data/progeny_models/full_matrix/full_progeny_matrix.rds")

# filter based on top x genes
# top 100 
m_100 = model %>% group_by(pathway) %>% top_n(100, -adj.p) %>% ungroup() %>%
  select(gene, pathway, weight=zscore)

# top 200 
m_200 = model %>% group_by(pathway) %>% top_n(200, -adj.p) %>% ungroup() %>%
  select(gene, pathway, weight=zscore)

# top 300 
m_300 = model %>% group_by(pathway) %>% top_n(300, -adj.p) %>% ungroup() %>%
  select(gene, pathway, weight=zscore)

# top 500 
m_500 = model %>% group_by(pathway) %>% top_n(500, -adj.p) %>% ungroup() %>%
  select(gene, pathway, weight=zscore)

# top 1000 
m_1000 = model %>% group_by(pathway) %>% top_n(1000, -adj.p) %>% ungroup() %>%
  select(gene, pathway, weight=zscore)

models = tribble(
  ~M, ~footprints,
  m_100, 100,
  m_200, 200,
  m_300, 300,
  m_500, 500,
  m_1000, 1000,
) %>%
  transmute(footprints, genesets = M %>% map(function(m) {
    m %>%
      group_by(pathway) %>%
      summarise(geneset = list(GeneSet(gene))) %>%
      transmute(pathway, geneset2 = pmap(., .f=function(pathway, geneset, ...) {
        setName(geneset) = pathway
        return(geneset)
        })) %>%
      deframe() %>%
      GeneSetCollection()
  })) %>%
  mutate(footprints = as_factor(footprints))


df = readRDS("output/hca_data/expression/norm.rds")
    
input = models %>%
  mutate(emat = list(df)) %>%
  unnest(emat, .preserve = genesets)

ss_paucell_scores = input %>%
  mutate(ss_paucell_result = pmap(., .f=function(norm, genesets, footprints, technology, ...) {
    message(footprints,"-",technology)
    
    g = subsetGeneSets(genesets, rownames(norm))
    obj = AUCell_buildRankings(norm, nCores=4, plotStats = F, verbose = F) %>%
      AUCell_calcAUC(genesets, ., verbose=F)
    
    res = obj@assays$data@listData$AUC
  })) %>%
  select(-norm, -genesets)

saveRDS(ss_paucell_scores, "output/hca_data/paucell/paucell_result.rds")
```

#### Run SCENIC
##### Build protocol dependent regulatory networks
```{r}
# path = "output/hca_data/scenic/complete/CEL-Seq2/int/2.6_regulons_asGeneSet.Rds"
scenic_regulons = list.files("output/hca_data/scenic/complete", recursive = T, 
                             full.names = T, pattern = "2.6_") %>%
  map_dfr(function(path) {
    technology = path %>% str_split(pattern = "/") %>% pluck(1,5)
    readRDS(path) %>%
      enframe("tf", "target") %>%
      mutate(technology = technology) %>%
      filter(!str_detect(tf, "extended")) %>%
      unnest(target) %>%
      select(technology, tf, target)
  }) %>%
  nest(-technology, .key="regulons") %>%
  mutate(technology = factor(
    technology, levels = c("Chromium","inDrop","C1HT-medium","C1HT-small",
                           "CEL-Seq2","ddSEQ","Drop-Seq","ICELL8","MARS-Seq",
                           "Chromium (sn)","Quartz-Seq2","SCRB-Seq",
                           "Smart-Seq2")))

saveRDS(scenic_regulons, "output/hca_data/scenic/complete/scenic_regulons.rds")
```

##### Run SCENIC
```{r}
scenic_regulons = readRDS("output/hca_data/scenic/complete/scenic_regulons.rds")

genesets = scenic_regulons %>%
  transmute(technology, geneset = regulons %>% map(function(regulons) {
    regulons %>%
      group_by(tf) %>%
      summarise(geneset = list(GeneSet(target, setName = unique(tf)))) %>%
      ungroup() %>%
      deframe() %>%
      GeneSetCollection()
  }))

df = readRDS("output/hca_data/expression/norm.rds")

input = df %>%
  inner_join(genesets, by="technology")


ss_scenic_scores = input %>%
  transmute(technology, ss_scenic_result = pmap(., .f= function(norm, geneset, ...) {
    obj = AUCell_buildRankings(norm, nCores=4, plotStats = F, verbose = F) %>%
      AUCell_calcAUC(geneset, ., verbose=F)
    
    res = obj@assays$data@listData$AUC
  }))


saveRDS(ss_scenic_scores, "output/hca_data/scenic/complete/scenic_result.rds")
```

#### Summarize footprint scores
```{r "summarize footprint scores}
meta = readRDS("output/hca_data/meta/meta_df.rds")
tf_confidence_mapping = read_csv("data/regulons/dorothea/dorothea_regulon_human_v1.csv") %>%
  distinct(feature=tf, confidence)
progeny = readRDS("output/hca_data/progeny/progeny_scores.rds") %>%
  mutate(tool = "progeny") %>%
  transmute(technology, footprints, tool, activity = progeny %>% map(function(progeny) {
    progeny %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("cell") %>%
      as_tibble() %>%
      gather(feature, activity, -cell)
  })) 

viper = readRDS("output/hca_data/viper/viper_scores.rds") %>%
  mutate(tool = "viper") %>%
  transmute(technology, tool, activity = viper %>% map(function(dorothea) {
    dorothea %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("feature") %>%
      as_tibble() %>%
      gather(cell, activity, -feature) %>%
      inner_join(tf_confidence_mapping, by="feature")
  })) 

metaviper = readRDS("output/hca_data/metaviper/metaviper_scores.rds") %>%
  mutate(tool = "metaviper") %>%
  transmute(technology, tool, activity = ss_metaviper_result %>% map(function(metaviper) {
    metaviper %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("feature") %>%
      as_tibble() %>%
      gather(cell, activity, -feature) %>%
      inner_join(tf_confidence_mapping, by="feature")
  })) 

daucell = readRDS("output/hca_data/daucell/daucell_result.rds") %>%
  mutate(tool = "daucell") %>%
  transmute(technology, tool, activity = ss_daucell_result %>% map(function(daucell) {
    daucell %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("feature") %>%
      as_tibble() %>%
      gather(cell, activity, -feature) %>%
      inner_join(tf_confidence_mapping, by="feature")
  })) 

paucell = readRDS("output/hca_data/paucell/paucell_result.rds") %>%
  mutate(tool = "paucell") %>%
  transmute(technology, footprints, tool, activity = ss_paucell_result %>% map(function(paucell) {
    paucell %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("feature") %>%
      as_tibble() %>%
      gather(cell, activity, -feature)
  })) 

scenic = readRDS("output/hca_data/scenic/complete/scenic_result.rds") %>%
  mutate(tool = "scenic") %>%
  transmute(technology,tool, activity = ss_scenic_result %>% map(function(scenic) {
    scenic %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("feature") %>%
      as_tibble() %>%
      gather(cell, activity, -feature)
  })) 

summarized_footprint_scores = bind_rows(progeny,viper, paucell, daucell, metaviper, scenic) %>%
  unnest() %>%
  inner_join(meta, by=c("cell", "technology")) %>%
  rename(celltype = nnet2) %>%
  group_by(technology, footprints, tool, feature, confidence, celltype) %>%
  summarise(m = mean(activity),
            s = sd(activity)) %>%
  ungroup()

saveRDS(summarized_footprint_scores, 
        "output/hca_data/summarized_footprint_scores/summarized_footprints_scores.rds")
write_csv(summarized_footprint_scores, 
          "output/hca_data/summarized_footprint_scores/summarized_footprints_scores.csv")
```

### Clustering
#### Setup hierarchy structure wo Megakaryocytes
```{r, "setup and visualize hierarchy structure}
# setup hierarchy structure
hierarchy = readRDS("output/hca_data/meta/meta_df.rds") %>%
  filter(!nnet2 %in% c("unclassified", "Megakaryocytes")) %>%
  distinct(cell, celltype = nnet2) %>%
  arrange(cell) %>%
  mutate(
    # hierarchy level 0
    hrchy0_id = as.integer(as_factor(celltype)),
    hrchy0_label = celltype,
    # hierarchy level 1
    hrchy1_id = case_when(
      celltype == "B cells" ~ 1,
      celltype == "NK cells" ~ 2,
      celltype %in% c("CD4 T cells","CD8 T cells") ~ 3,
      celltype == "Dendritic cells" ~ 4,
      celltype %in% c("CD14+ Monocytes","FCGR3A+ Monocytes") ~ 5,
      celltype == "HEK cells" ~ 6),
    hrchy1_label = case_when(
      celltype == "B cells" ~ "B cells",
      celltype == "NK cells" ~ "NK cells",
      celltype %in% c("CD4 T cells","CD8 T cells") ~ "T cells",
      celltype == "Dendritic cells" ~ "Dendritic cells",
      celltype %in% c("CD14+ Monocytes","FCGR3A+ Monocytes") ~ "Monocytes",
      celltype == "HEK cells" ~ "HEK cells"),
    # hierarchy level 2
    hrchy2_id = case_when(
      celltype == "B cells" ~ 1,
      celltype %in% c("CD4 T cells","CD8 T cells","NK cells") ~ 2,
      celltype %in% c("CD14+ Monocytes","FCGR3A+ Monocytes","Dendritic cells") ~ 
        3,
      celltype == "HEK cells" ~ 4),
    hrchy2_label = case_when(
      celltype == "B cells" ~ "B cells",
      celltype %in% c("CD4 T cells","CD8 T cells","NK cells") ~ "T + NK cells",
      celltype %in% c("CD14+ Monocytes","FCGR3A+ Monocytes","Dendritic cells") ~
        "Monocytes + Dendritic",
      celltype == "HEK cells" ~ "HEK cells"),
    # hierarchy level 3
    hrchy3_id = case_when(
      celltype %in% c("CD4 T cells","CD8 T cells","NK cells", "B cells") ~ 1,
      celltype %in% c("CD14+ Monocytes","FCGR3A+ Monocytes","Dendritic cells") ~ 2,
      celltype == "HEK cells" ~ 3),
    hrchy3_label = case_when(
      celltype %in% c("CD4 T cells","CD8 T cells","NK cells", "B cells") ~ 
        "unknown",
      celltype %in% c("CD14+ Monocytes","FCGR3A+ Monocytes","Dendritic cells") ~ "Myeloid",
      celltype == "HEK cells" ~ "HEK cells"),
    # hierarchy level 4
    hrchy4_id = case_when(
      celltype == "HEK cells" ~ 1,
      TRUE ~ 2),
    hrchy4_label = case_when(
      celltype == "HEK cells" ~ "cell line",
      TRUE ~ "cell type")) %>%
  gather(key, value, -cell, -celltype) %>%
  separate(key, into=c("hrchy_level","key"), sep = "_") %>% 
  spread(key, value) %>%
  mutate(id = as.integer(id)) %>%
  nest(-hrchy_level, .key="cluster_vec")

saveRDS(hierarchy, "output/hca_data/umap/cluster_hierarchy.rds")
```

#### Prepare different input matrices for clustering with UMAP
```{r "prepare different input matrices for clustering with UMAP}
# set seed
set.seed(123)

missing_tfs = readRDS("output/hca_data/meta/tfs_available_in_dorothea_but_not_gtex.rds")
d = read_csv("data/regulons/dorothea/dorothea_regulon_human_v1.csv")

tfs_A = d %>%
  filter(confidence %in% c("A")) %>%
  distinct(tf) %>%
  anti_join(missing_tfs)

tfs_AB = d %>%
  filter(confidence %in% c("A", "B")) %>%
  distinct(tf) %>%
  anti_join(missing_tfs)

tfs_ABC = d %>%
  filter(confidence %in% c("A", "B", "C")) %>%
  distinct(tf) %>%
  anti_join(missing_tfs)

tfs_ABCD = d %>%
  filter(confidence %in% c("A", "B", "C", "D")) %>%
  distinct(tf) %>%
  anti_join(missing_tfs)

tfs_ABCDE = d %>%
  distinct(tf) %>%
  anti_join(missing_tfs)

df = readRDS("output/hca_data/expression/norm.rds")


norm_expr = df %>%
  mutate(class = "normalized_expression") %>%
  rename(mat = norm)

saveRDS(norm_expr,"output/hca_data/umap/umap_input/norm_expr.rds")

most_var_norm = df %>%
  transmute(technology, mat = norm %>% future_map(function(norm) {
    most_var_genes = norm %>% 
      apply(.,1,var) %>%
      enframe("gene", "variance") %>%
      top_n(2000, variance) %>%
      pull(gene)
    
    norm[most_var_genes, ]
  }, .progress = T)) %>%
  mutate(class = "most_var_normalized_expr")

saveRDS(most_var_norm,"output/hca_data/umap/umap_input/most_var_norm_expr.rds")

cv2 = df %>%
  transmute(technology, mat = norm %>% future_map(function(norm) {
    sd_df = norm %>% 
          apply(.,1,sd) %>%
          enframe("gene", "sd")
    
    mean_df = norm %>%
      apply(., 1, mean) %>%
      enframe("gene", "mean")
    
    genes = left_join(sd_df, mean_df, by="gene") %>%
      mutate(cv = (sd/mean)**2) %>%
      top_n(2000, -cv) %>%
      pull(gene)
    
    norm[genes, ]
  }, .progress = T)) %>%
  mutate(class = "most_cv_normalized_expr")

saveRDS(cv2,"output/hca_data/umap/umap_input/most_cv_norm_expr.rds")

seurat_2000_vst = df %>%
  transmute(technology, mat = future_pmap(., .f= function(norm, technology, ...) {
    print(technology)
    rownames(norm) = str_replace_all(rownames(norm), pattern = "_", "-")
    
    # NOTE feature names are changed: "_" to "-"
    genes = CreateSeuratObject(norm) %>%
      FindVariableFeatures(., selection.method = "vst", nfeatures = 2000) %>%
      VariableFeatures()
  
    norm[genes, ]
  }, .progress = T)) %>%
  mutate(class = "top_2000_seurat_vst")

saveRDS(seurat_2000_vst,"output/hca_data/umap/umap_input/top_2000_seurat_vst.rds")

seurat_2000_mvp = df %>%
  transmute(technology, mat = future_pmap(., .f= function(norm, technology, ...) {
    print(technology)
    rownames(norm) = str_replace_all(rownames(norm), pattern = "_", "-")
    
    # NOTE feature names are changed: "_" to "-"
    genes = CreateSeuratObject(norm) %>%
      FindVariableFeatures(., selection.method = "mvp", nfeatures = 2000) %>%
      VariableFeatures()
  
    norm[genes, ]
  }, .progress = T)) %>%
  mutate(class = "top_2000_seurat_mvp")

saveRDS(seurat_2000_mvp,"output/hca_data/umap/umap_input/top_2000_seurat_mvp.rds")

seurat_2000_disp = df %>%
  transmute(technology, mat = future_pmap(., .f= function(norm, technology, ...) {
    print(technology)
    rownames(norm) = str_replace_all(rownames(norm), pattern = "_", "-")
    
    # NOTE feature names are changed: "_" to "-"
    genes = CreateSeuratObject(norm) %>%
      FindVariableFeatures(., selection.method = "disp", nfeatures = 2000) %>%
      VariableFeatures()
  
    norm[genes, ]
  }, .progress = T)) %>%
  mutate(class = "top_2000_seurat_disp")

saveRDS(seurat_2000_disp,"output/hca_data/umap/umap_input/top_2000_seurat_disp.rds")

set.seed(12345)
random_seurat_113 = df %>%
  transmute(technology, mat = future_pmap(., .f= function(norm, technology, ...) {
    print(technology)
    rownames(norm) = str_replace_all(rownames(norm), pattern = "_", "-")
    
    # NOTE feature names are changed: "_" to "-"
    genes = CreateSeuratObject(norm) %>%
      FindVariableFeatures(., selection.method = "vst", nfeatures = 2000) %>%
      VariableFeatures()
    
    map(seq_len(10), function(rep) {
      rand_genes = sample(genes, 113)
      norm[rand_genes,]
    }) %>%
      enframe("rep", "mat")
  }, .progress = T)) %>%
  mutate(class = "random_113_seurat_mvg") %>%
  unnest(mat)

saveRDS(random_seurat_113,"output/hca_data/umap/umap_input/random_113_seurat_mvg.rds")

set.seed(12345)
random_seurat_14 = df %>%
  transmute(technology, mat = future_pmap(., .f= function(norm, technology, ...) {
    print(technology)
    rownames(norm) = str_replace_all(rownames(norm), pattern = "_", "-")
    
    # NOTE feature names are changed: "_" to "-"
    genes = CreateSeuratObject(norm) %>%
      FindVariableFeatures(., selection.method = "vst", nfeatures = 2000) %>%
      VariableFeatures()
    
    map(seq_len(10), function(rep) {
      rand_genes = sample(genes, 14)
      norm[rand_genes,]
    }) %>%
      enframe("rep", "mat")
  }, .progress = T)) %>%
  mutate(class = "random_14_seurat_mvg") %>%
  unnest(mat)

saveRDS(random_seurat_14,"output/hca_data/umap/umap_input/random_14_seurat_mvg.rds")


seurat_113 = df %>%
  transmute(technology, mat = future_pmap(., .f= function(norm, technology, ...) {
    print(technology)
    rownames(norm) = str_replace_all(rownames(norm), pattern = "_", "-")
    
    # NOTE feature names are changed: "_" to "-"
    genes = CreateSeuratObject(norm) %>%
      FindVariableFeatures(., selection.method = "vst", nfeatures = 113) %>%
      VariableFeatures()
  
    norm[genes, ]
  }, .progress = T)) %>%
  mutate(class = "top_113_seurat_mvg")

saveRDS(seurat_113,"output/hca_data/umap/umap_input/top_113_seurat_mvg.rds")

seurat_14 = df %>%
  transmute(technology, mat = future_pmap(., .f= function(norm, technology, ...) {
    print(technology)
    rownames(norm) = str_replace_all(rownames(norm), pattern = "_", "-")
    
    # NOTE feature names are changed: "_" to "-"
    genes = CreateSeuratObject(norm) %>%
      FindVariableFeatures(., selection.method = "vst", nfeatures = 14) %>%
      VariableFeatures()
  
    norm[genes, ]
  }, .progress = T)) %>%
  mutate(class = "top_14_seurat_mvg")

saveRDS(seurat_14,"output/hca_data/umap/umap_input/top_14_seurat_mvg.rds")

progeny = readRDS("output/hca_data/progeny/progeny_scores.rds") %>%
  mutate(class = str_c("progeny",footprints)) %>%
  rename(mat = progeny) %>%
  select(-footprints)

saveRDS(progeny,"output/hca_data/umap/umap_input/progeny.rds")

# VIPER ------------------------------------------------------------------------
# A - contains activities for 93 TFs
viperA = readRDS("output/hca_data/viper/viper_scores.rds") %>%
  mutate(viper = viper %>% map(function(viper) {
    viper %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("tf") %>%
      inner_join(tfs_A, by="tf") %>%
      data.frame(row.names = 1, stringsAsFactors = F, check.names = F) %>%
      as.matrix()
  })) %>%
  mutate(class = "viperA") %>%
  rename(mat=viper)

saveRDS(viperA,"output/hca_data/umap/umap_input/viperA.rds")

# AB - contains activities for 113 TFs
viperAB = readRDS("output/hca_data/viper/viper_scores.rds") %>%
  mutate(viper = viper %>% map(function(viper) {
    viper %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("tf") %>%
      inner_join(tfs_AB, by="tf") %>%
      data.frame(row.names = 1, stringsAsFactors = F, check.names = F) %>%
      as.matrix()
  })) %>%
  mutate(class = "viperAB") %>%
  rename(mat=viper)

saveRDS(viperAB,"output/hca_data/umap/umap_input/viperAB.rds")

# ABC - contains activities for 278 TFs
viperABC = readRDS("output/hca_data/viper/viper_scores.rds") %>%
  mutate(viper = viper %>% map(function(viper) {
    viper %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("tf") %>%
      inner_join(tfs_ABC, by="tf") %>%
      data.frame(row.names = 1, stringsAsFactors = F, check.names = F) %>%
      as.matrix()
  })) %>%
  mutate(class = "viperABC") %>%
  rename(mat=viper)

saveRDS(viperABC,"output/hca_data/umap/umap_input/viperABC.rds")

# ABCD - contains activities for 369 TFs
viperABCD = readRDS("output/hca_data/viper/viper_scores.rds") %>%
  mutate(viper = viper %>% map(function(viper) {
    viper %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("tf") %>%
      inner_join(tfs_ABCD, by="tf") %>%
      data.frame(row.names = 1, stringsAsFactors = F, check.names = F) %>%
      as.matrix()
  })) %>%
  mutate(class = "viperABCD") %>%
  rename(mat=viper)

saveRDS(viperABCD,"output/hca_data/umap/umap_input/viperABCD.rds")

# ABCDE - contains activities for 1299 TFs
viperABCDE = readRDS("output/hca_data/viper/viper_scores.rds") %>%
  mutate(viper = viper %>% map(function(viper) {
    viper %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("tf") %>%
      inner_join(tfs_ABCDE, by="tf") %>%
      data.frame(row.names = 1, stringsAsFactors = F, check.names = F) %>%
      as.matrix()
  })) %>%
  mutate(class = "viperABCDE") %>%
  rename(mat=viper)

saveRDS(viperABCDE,"output/hca_data/umap/umap_input/viperABCDE.rds")

# metaVIPER --------------------------------------------------------------------
# A - contains activities for 93 TFs
metaviperA = readRDS("output/hca_data/metaviper/metaviper_scores.rds") %>%
  mutate(metaviper = ss_metaviper_result %>% map(function(metaviper) metaviper[tfs_A$tf,])) %>%
  select(-ss_metaviper_result) %>%
  mutate(class = "metaviperA") %>%
  rename(mat=metaviper)

saveRDS(metaviperA,"output/hca_data/umap/umap_input/metaviperA.rds")

# AB - contains activities for 113 TFs
metaviperAB = readRDS("output/hca_data/metaviper/metaviper_scores.rds") %>%
  mutate(metaviper = ss_metaviper_result %>% map(function(metaviper) metaviper[tfs_AB$tf,])) %>%
  select(-ss_metaviper_result) %>%
  mutate(class = "metaviperAB") %>%
  rename(mat=metaviper)

saveRDS(metaviperAB,"output/hca_data/umap/umap_input/metaviperAB.rds")

# ABC - contains activities for 278 TFs
metaviperABC = readRDS("output/hca_data/metaviper/metaviper_scores.rds") %>%
  mutate(metaviper = ss_metaviper_result %>% map(function(metaviper) metaviper[tfs_ABC$tf,])) %>%
  select(-ss_metaviper_result) %>%
  mutate(class = "metaviperABC") %>%
  rename(mat=metaviper)

saveRDS(metaviperABC,"output/hca_data/umap/umap_input/metaviperABC.rds")

# ABCD - contains activities for 369 TFs
metaviperABCD = readRDS("output/hca_data/metaviper/metaviper_scores.rds") %>%
  mutate(metaviper = ss_metaviper_result %>% map(function(metaviper) metaviper[tfs_ABCD$tf,])) %>%
  select(-ss_metaviper_result) %>%
  mutate(class = "metaviperABCD") %>%
  rename(mat=metaviper)

saveRDS(metaviperABCD,"output/hca_data/umap/umap_input/metaviperABCD.rds")

# ABCDE - contains activities for 1299 TFs
metaviperABCDE = readRDS("output/hca_data/metaviper/metaviper_scores.rds") %>%
  mutate(metaviper = ss_metaviper_result %>% map(function(metaviper) metaviper[tfs_ABCDE$tf,])) %>%
  select(-ss_metaviper_result) %>%
  mutate(class = "metaviperABCDE") %>%
  rename(mat=metaviper)

saveRDS(metaviperABCDE,"output/hca_data/umap/umap_input/metaviperABCDE.rds")


# AUCell with DoRothEA ---------------------------------------------------------
# A
daucellA = readRDS("output/hca_data/daucell/daucell_result.rds") %>%
  mutate(ss_daucell_result = ss_daucell_result %>% map(function(daucell) {
    daucell %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("tf") %>%
      inner_join(tfs_A, by="tf") %>%
      data.frame(row.names = 1, stringsAsFactors = F, check.names = F) %>%
      as.matrix()
  })) %>%
  mutate(class = "daucellA") %>%
  rename(mat=ss_daucell_result)

saveRDS(daucellA,"output/hca_data/umap/umap_input/daucellA.rds")

# AB
daucellAB = readRDS("output/hca_data/daucell/daucell_result.rds") %>%
  mutate(ss_daucell_result = ss_daucell_result %>% map(function(daucell) {
    daucell %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("tf") %>%
      inner_join(tfs_AB, by="tf") %>%
      data.frame(row.names = 1, stringsAsFactors = F, check.names = F) %>%
      as.matrix()
  })) %>%
  mutate(class = "daucellAB") %>%
  rename(mat=ss_daucell_result)

saveRDS(daucellAB,"output/hca_data/umap/umap_input/daucellAB.rds")

# ABC
daucellABC = readRDS("output/hca_data/daucell/daucell_result.rds") %>%
  mutate(ss_daucell_result = ss_daucell_result %>% map(function(daucell) {
    daucell %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("tf") %>%
      inner_join(tfs_ABC, by="tf") %>%
      data.frame(row.names = 1, stringsAsFactors = F, check.names = F) %>%
      as.matrix()
  })) %>%
  mutate(class = "daucellABC") %>%
  rename(mat=ss_daucell_result)

saveRDS(daucellABC,"output/hca_data/umap/umap_input/daucellABC.rds")

# ABCD
daucellABCD = readRDS("output/hca_data/daucell/daucell_result.rds") %>%
  mutate(ss_daucell_result = ss_daucell_result %>% map(function(daucell) {
    daucell %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("tf") %>%
      inner_join(tfs_ABCD, by="tf") %>%
      data.frame(row.names = 1, stringsAsFactors = F, check.names = F) %>%
      as.matrix()
  })) %>%
  mutate(class = "daucellABCD") %>%
  rename(mat=ss_daucell_result)

saveRDS(daucellABCD,"output/hca_data/umap/umap_input/daucellABCD.rds")

# ABCDE
daucellABCDE = readRDS("output/hca_data/daucell/daucell_result.rds") %>%
  mutate(ss_daucell_result = ss_daucell_result %>% map(function(daucell) {
    daucell %>%
      data.frame(check.names = F, stringsAsFactors = F) %>%
      rownames_to_column("tf") %>%
      inner_join(tfs_ABCDE, by="tf") %>%
      data.frame(row.names = 1, stringsAsFactors = F, check.names = F) %>%
      as.matrix()
  })) %>%
  mutate(class = "daucellABCDE") %>%
  rename(mat=ss_daucell_result)

saveRDS(daucellABCDE,"output/hca_data/umap/umap_input/daucellABCDE.rds")

# AUCell with PROGENy gene sets ------------------------------------------------
paucell = readRDS("output/hca_data/paucell/paucell_result.rds") %>%
  mutate(class = str_c("paucell",footprints)) %>%
  rename(mat = ss_paucell_result) %>%
  select(-footprints)

saveRDS(paucell,"output/hca_data/umap/umap_input/paucell.rds")

# SCENIC -----------------------------------------------------------------------
scenic = readRDS("output/hca_data/scenic/scenic_result.rds") %>%
  mutate(class = "scenic") %>%
  dplyr::rename(mat = ss_scenic_result)

# GSEA with GO terms -----------------------------------------------------------
gogsea = readRDS("output/hca_data/gogsea/gogsea_scores.rds") %>%
  mutate(class = "gogsea") %>%
  rename(mat = ss_gogsea_scores)

saveRDS(gogsea,"output/hca_data/umap/umap_input/gogsea.rds")

raw_data = readRDS("output/hca_data/expression/raw_preprocessed.rds") %>%
  mutate(class = "raw_expression") %>%
  rename(mat = raw)

saveRDS(raw_data,"output/hca_data/umap/umap_input/raw_data.rds")

footprint_genes = progeny_matrix_human_v1 %>% distinct(gene) %>% pull()
footprints = df %>%
  transmute(technology, mat = norm %>% map(function(norm) {
    z = intersect(rownames(norm), footprint_genes)
    norm[z,]
  })) %>%
  mutate(class = "footprint_genes")

saveRDS(footprints,"output/hca_data/umap/umap_input/footprints.rds")

tf_genes = dorothea_regulon_human_v1 %>% 
  filter(confidence %in% c("A", "B")) %>%
  anti_join(missing_tfs) %>%
  distinct(tf) %>% 
  pull()

tf_expression = df %>%
  transmute(technology, mat = norm %>% map(function(norm) {
    z = intersect(rownames(norm), tf_genes)
    norm[z,]
  })) %>%
  mutate(class = "tf_genes")

saveRDS(tf_expression,"output/hca_data/umap/umap_input/tf_expression.rds")

set.seed(1234)
random_1301_genes = df
random_1301_genes_res = map_df(seq_len(10), ~random_1301_genes) %>%
  mutate(rep = rep(seq_len(10), each=13)) %>%
  group_by(rep) %>%
  mutate(seed = sample(1:1e6,1)) %>% 
  ungroup() %>%
  transmute(technology, rep, mat = map2(norm, seed, function(norm, seed) {
    set.seed(seed)
    genes = rownames(norm) %>%
      sample(1301)
    norm[genes,]
  })) %>%
  mutate(class = "random_footprints")

saveRDS(random_1301_genes_res,"output/hca_data/umap/umap_input/random_1301_genes.rds")

set.seed(12345)
random_113_genes = df
random_113_genes_res = map_df(seq_len(10), ~random_113_genes) %>%
  mutate(rep = rep(seq_len(10), each=13)) %>%
  group_by(rep) %>%
  mutate(seed = sample(1:1e6,1)) %>% 
  ungroup() %>%
  transmute(technology, rep, mat = map2(norm, seed, function(norm, seed) {
    set.seed(seed)
    genes = rownames(norm) %>%
      sample(113)
    norm[genes,]
  })) %>%
  mutate(class = "random_113_features")

saveRDS(random_113_genes_res,"output/hca_data/umap/umap_input/random_113_genes.rds")

set.seed(123456)
random_14_genes = df
random_14_genes_res = map_df(seq_len(10), ~random_14_genes) %>%
  mutate(rep = rep(seq_len(10), each=13)) %>%
  group_by(rep) %>%
  mutate(seed = sample(1:1e6,1)) %>% 
  ungroup() %>%
  transmute(technology, rep, mat = map2(norm, seed, function(norm, seed) {
    set.seed(seed)
    genes = rownames(norm) %>%
      sample(14)
    norm[genes,]
  })) %>%
  mutate(class = "random_14_features")
saveRDS(random_14_genes_res,"output/hca_data/umap/umap_input/random_14_genes.rds")
```

#### Application of UMAP
```{r "clustering"}
do_umap = function(mat, meta, class, technology,...) {
  # browser()
  message(technology, " - ", class)
  if (str_detect(class, "progeny")) {
    mat = t(mat)
  }
  cells = meta %>% pull(cell)
  stopifnot(colnames(mat[,cells]) == meta$cell)
  
  mat[,cells] %>%
    as.matrix() %>%
    t() %>%
    umap(method="umap-learn") %>%
    pluck("layout") %>%
    data.frame(check.names = F, stringsAsFactors = F) %>%
    rownames_to_column("cell") %>%
    rename(x = `1`, y=`2`) %>%
    as_tibble() %>%
    left_join(hrchy, by="cell") %>%
    arrange(cell)
}
meta = readRDS("output/hca_data/meta/meta_df.rds") %>%
  filter(!nnet2 %in% c("unclassified", "Megakaryocytes")) %>%
  nest(-technology, .key="meta")

hrchy = readRDS("output/hca_data/umap/cluster_hierarchy.rds") %>%
  unnest(cluster_vec) %>%
  select(-id) %>%
  spread(hrchy_level, label)

df = list.files("output/hca_data/umap/umap_input", full.names = T) %>%
  future_map_dfr(function(path) {
    df = readRDS(path) %>%
      left_join(meta)
    coords = df %>%
      mutate(coords = pmap(., .f=do_umap)) %>%
      select(-mat, -meta)
    save_path = str_replace(path, "umap_input", "umap_coords")
    saveRDS(coords, save_path)
  }, .progress = T)

res = list.files("output/hca_data/umap/umap_coords", full.names = T) %>%
  map_df(readRDS)

saveRDS(res, "output/hca_data/umap/umap_coords.rds")
```

#### Silhouettes
```{r "silhouettes"}
do_sil = function(coords, technology, class,...) {
  # browser()
  message(technology, "-", class)
  
  res = h %>%
    transmute(hrchy_level, sil = pmap(., .f=function(cluster_vec, hrchy_level,...) {
      message(hrchy_level)
      cv = cluster_vec %>%
        filter(cell %in% coords$cell)
      meta_cluster = cv %>%
        distinct(label, id) %>%
        mutate_if(is.integer, as.character)
        
      stopifnot(coords$cell == cv$cell)
      dist_mat = coords %>%
        select_if(is.numeric) %>%
        dist()
  
      s = silhouette(cv$id, dist_mat)
      summary_s = s %>% summary()
    
    
      summary_s %>%
          pluck("clus.avg.widths") %>%
          enframe("id", "avg_width") %>%
          left_join(meta_cluster) %>%
          mutate(mean_avg_width = pluck(summary_s, "avg.width"),
                 size = as.vector(pluck(summary_s, "clus.sizes")),
                 obj = list(s))
    })) %>%
    unnest(sil)
}

umap_coords = readRDS("output/hca_data/umap/umap_coords.rds")
h = readRDS("output/hca_data/umap/cluster_hierarchy.rds")

sil_result = umap_coords %>%
  # filter(class == "scenic") %>%
  # mutate(sil = pmap(., .f=do_sil))
  mutate(sil = future_pmap(., .f=do_sil, .progress = T))

saveRDS(sil_result, "output/hca_data/umap/silhouette_results.rds")
```



#### Prepare silhouette correlation
```{r "prepare silhouette correlation"}
sil_res = readRDS("output/hca_data/umap/silhouette_results.rds")

cor = sil_res %>%
  filter(!class %in% c("raw_expression", "most_var_raw_expr")) %>%
  mutate(gg = case_when(class == "top_2000_seurat_vst" ~ "reference",
                        class != "top_2000_seurat_vst" ~ "case")) %>%
    # mutate(gg = case_when(class == "normalized_expression" ~ "reference",
    #                     class != "normalized_expression" ~ "case")) %>%
  mutate(label = case_when(
    class == "random_14_features" ~ "Old Neg. control",
    class == "footprint_genes" ~ "Footprint expression",
    class == "random_footprints" ~ "Old Neg. control",
    class == "viperA" ~ "DoRothEA (A) ",
    class == "viperAB" ~ "DoRothEA (AB)",
    class == "viperABC" ~ "DoRothEA (ABC)",
    class == "viperABCD" ~ "DoRothEA (ABCD)",
    class == "viperABCDE" ~ "DoRothEA (ABCDE)",
    class == "metaviperA" ~ "metaVIPER (A)",
    class == "metaviperAB" ~ "metaVIPER (AB)",
    class == "metaviperABC" ~ "metaVIPER (ABC)",
    class == "metaviperABCD" ~ "metaVIPER (ABCD)",
    class == "metaviperABCDE" ~ "metaVIPER (ABCDE)",
    class == "daucellA" ~ "D-AUCell (A)",
    class == "daucellAB" ~ "D-AUCell (AB)",
    class == "daucellABC" ~ "D-AUCell (ABC)",
    class == "daucellABCD" ~ "D-AUCell (ABCD)",
    class == "daucellABCDE" ~ "D-AUCell (ABCDE)",
    class == "random_113_features" ~ "Old Neg. control",
    class == "tf_genes" ~ "TF expression (AB)",
    class == "progeny100" ~ "PROGENy (100)",
    class == "progeny200" ~ "PROGENy (200)",
    class == "progeny300" ~ "PROGENy (300)",
    class == "progeny500" ~ "PROGENy (500)",
    class == "progeny1000" ~ "PROGENy (1000)",
    class == "progenyall" ~ "PROGENy (All)",
    # class == "gogsea" ~ "GSEA (GO-terms)",
    class == "paucell100" ~ "P-AUCell (100)",
    class == "paucell200" ~ "P-AUCell (200)",
    class == "paucell300" ~ "P-AUCell (300)",
    class == "paucell500" ~ "P-AUCell (500)",
    class == "paucell1000" ~ "P-AUCell (1000)",
    class == "most_var_normalized_expr" ~ "MVG",
    class == "most_cv_normalized_expr" ~ "CV",
    class == "top_113_seurat_mvg" ~ "Pos. control",
    class == "random_113_seurat_mvg" ~ "Neg. control",
    class == "top_14_seurat_mvg" ~ "Pos. control",
    class == "random_14_seurat_mvg" ~ "Neg. control",
    class == "top_2000_seurat_vst" ~ "Seurat (vst)",
    class == "top_2000_seurat_mvp" ~ "Seurat (mvp)",
    class == "top_2000_seurat_disp" ~ "Seurat (disp)",
    class == "normalized_expression" ~ "Normalized expression",
    class == "scenic" ~ "SCENIC")) %>%
  mutate(label = as_factor(label)) %>%
  mutate(label = fct_relevel(label, "Neg. control", after=Inf)) %>%
  unnest(sil) %>%
  group_by(technology, gg, class, label, hrchy_level) %>%
  summarise(m = mean(mean_avg_width)) %>%
  ungroup()


saveRDS(cor, "output/hca_data/umap/sil_correlation_input.rds")
```
#### Statistical evaluation of silhouette plot
```{r}
df = readRDS("output/hca_data/umap/sil_correlation_input.rds")
class_label_mapping = df %>% distinct(feature = class, label)
technologies = df %>% pull(technology) %>% unique()

tf_stats = df %>%
  filter(class %in% c("top_2000_seurat_vst","daucellAB","metaviperAB",
                      "random_113_seurat_mvg","tf_genes",
                      "top_113_seurat_mvg","viperAB", "scenic")) %>%
  semi_join(class_label_mapping) %>%
  mutate(class = label) %>%
  select(-label) %>%
  mutate(class = fct_relevel(class, "Pos. control"),
         technology = fct_relevel(technology, "Quartz-Seq2")) %>%
  nest(-hrchy_level) %>%
  transmute(hrchy_level, model = data %>% map(~lm(m~technology + class, data = .))) %>%
  mutate(aov = model %>% map(~aov(.))) %>%
  mutate(tukey = aov %>% map(~TukeyHSD(.))) %>%
  mutate(stats = tukey %>% map(~tidy(.))) %>%
  mutate(tidy_lm = model %>% map(function(i) {
    i %>%
      tidy() %>%
      filter(term != "(Intercept)") %>%
      separate(term, into = c("tmp", "feature"), sep="technology") %>%
      separate(tmp, into = c("tmp1", "feature1"), sep="class") %>%
      mutate(feature = as_factor(coalesce(feature, feature1))) %>%
      select(-tmp1, -feature1) %>%
      mutate(group = case_when(feature %in% technologies ~ "technology",
                               TRUE ~ "class"))
  }))


saveRDS(tf_stats, "output/hca_data/umap/tf_stats.rds")


pathway_stats = df %>%
  filter(class %in% c("paucell200","progeny500","random_14_seurat_mvg",
                 "top_14_seurat_mvg","top_2000_seurat_vst")) %>%
  semi_join(class_label_mapping) %>%
  mutate(class = label) %>%
  select(-label) %>%
  mutate(class = fct_relevel(class, "Pos. control"),
         technology = fct_relevel(technology, "Quartz-Seq2")) %>%
  nest(-hrchy_level) %>%
  transmute(hrchy_level, model = data %>% map(~lm(m~technology + class, data = .))) %>%
  mutate(aov = model %>% map(~aov(.))) %>%
  mutate(tukey = aov %>% map(~TukeyHSD(.))) %>%
  mutate(stats = tukey %>% map(~tidy(.))) %>%
    mutate(tidy_lm = model %>% map(function(i) {
    i %>%
      tidy() %>%
      filter(term != "(Intercept)") %>%
      separate(term, into = c("tmp", "feature"), sep="technology") %>%
      separate(tmp, into = c("tmp1", "feature1"), sep="class") %>%
      mutate(feature = as_factor(coalesce(feature, feature1))) %>%
      select(-tmp1, -feature1) %>%
      mutate(group = case_when(feature %in% technologies ~ "technology",
                               TRUE ~ "class"))
  }))


saveRDS(pathway_stats, "output/hca_data/umap/pathway_stats.rds")

```

#### Correlation of TF activity scores across scRNA-seq protocols
```{r}
df = readRDS("output/hca_data/summarized_footprint_scores/summarized_footprints_scores.rds") %>%
  filter(!tool %in% c("progeny", "paucell")) %>%
  filter(!celltype %in% c("unclassified", "Megakaryocytes")) %>%
  select(-s, -footprints)

cor_input = df %>%
  unite(key, celltype, feature, remove = F) %>%
  group_by(technology, tool, feature) %>%
  mutate(z = scale(m)) %>%
  ungroup() 

cor_res_high_conf = cor_input %>%
  nest(-tool) %>%
  mutate(cor_df = data %>% map(function(data) {
    data %>%
      filter(confidence %in% c("A","B") | is.na(confidence)) %>%
      select(technology, key, z) %>%
      spread(technology, z, fill=NA)
  })) %>%
  mutate(r = cor_df %>% map(function(cor_df) {
    cor_df %>%
      select_if(is.numeric) %>%
      correlate() %>%
      shave() %>%
      stretch() %>%
      drop_na()
  })) 

cor_res_low_conf = cor_input %>%
  nest(-tool) %>%
  mutate(cor_df = data %>% map(function(data) {
    data %>%
      filter(confidence %in% c("C","D", "E") | is.na(confidence)) %>%
      select(technology, key, z) %>%
      spread(technology, z, fill=NA)
  })) %>%
  mutate(r = cor_df %>% map(function(cor_df) {
    cor_df %>%
      select_if(is.numeric) %>%
      correlate() %>%
      shave() %>%
      stretch() %>%
      drop_na()
  })) 

saveRDS(cor_res_high_conf, "output/hca_data/tf_correlations/tf_correlation_high_conf.rds")
saveRDS(cor_res_low_conf, "output/hca_data/tf_correlations/tf_correlation_low_conf.rds")
```
### Blood cancer analysis
#### TF activity vs gene essentiality
```{python}
import pandas as pd
import numpy as np
from scipy.stats import pearsonr as pcor

### read latest HCA data
hca_tf_activity = pd.read_csv('output/hca_data/summarized_footprint_scores/summarized_footprints_scores.csv', 
                              sep=',',header=0,index_col=None)

# remove HEK, Megakaryocytes and unlassified cells and remove the column footprints
hca_tf_activity = hca_tf_activity.drop(hca_tf_activity[(hca_tf_activity.celltype == "HEK cells") | (hca_tf_activity.celltype == "Megakaryocytes") | (hca_tf_activity.celltype == "unclassified") | (hca_tf_activity.tool == "progeny") | (hca_tf_activity.tool == "paucell")].index).drop(["footprints"], axis=1)


### get all tfs in our analysis
all_tfs = hca_tf_activity['feature'].unique()

### read achilles data
achilles = pd.read_csv('data/hca_data/blood_cancer/D2_combined_gene_dep_scores.csv',sep=',',header=0,index_col=0)

### set achilles index to gene names
achilles.index=pd.Series(achilles.index).apply(lambda x:x.split()[0]).values

### get ccle gene exspression values
gex = pd.read_csv('data/hca_data/blood_cancer/CCLE_expression.csv',sep=',',header=0,index_col=0)

### set ccle columns to gene names
gex.columns = pd.Series(gex.columns).apply(lambda x:x.split()[0]).values

### we concentrate on tfs shared in HCA analysis, Achilles and CCLE
shared_tfs = list(set(all_tfs) & set(gex.columns) & set(achilles.index))
gex = gex[shared_tfs]
achilles = achilles.loc[shared_tfs]

cell_anno = pd.read_csv('data/hca_data/blood_cancer/sample_info.csv',sep=',',header=0,index_col=0)

### set gex index to ccle name
gex.index = cell_anno['CCLE Name'][gex.index].values

### shared cell lines between achilles and gex
all_cells = list(set(achilles.columns) & set(gex.index))
gex = gex.loc[all_cells]
achilles = achilles[all_cells].T

gex_achilles_corr=[]
for tf in achilles.columns:
    cells=achilles.index[~achilles[tf].isnull()]
    gex_achilles_corr.append(pcor(achilles.loc[cells,tf],gex.loc[cells,tf])[0])
gex_achilles_corr=pd.DataFrame(pd.Series(gex_achilles_corr,index=achilles.columns))
gex_achilles_corr.columns=['correlation']

gex_achilles_corr.sort_values('correlation')

gex_achilles_corr.to_csv('output/hca_data/blood_cancer/gex_achilles_corr.csv',sep=',')

# multiply achilles with sign of correlation
achilles=achilles*np.sign(gex_achilles_corr['correlation'])
achilles.to_csv('output/hca_data/blood_cancer/tf_achilles_sign_corrected.csv',sep=',')

########
tf_ess = pd.read_csv('output/hca_data/blood_cancer/tf_achilles_sign_corrected.csv',sep=',',header=0,index_col=0)

# Load manually created annotation file for blood cancers. We will use the myeloid_lymphid columns to see the differences between myeloid and lymphid cells.
blood_anno = pd.read_csv('data/hca_data/blood_cancer/blood_cancer_anno.csv',
                         sep=',',header=0,index_col=None)
blood_anno.index=blood_anno['CCLE Name']
fil=blood_anno['myeloid_lymphoid']=='M'
blood_anno.loc[blood_anno.index[fil],'myeloid_lymphoid']='Myeloid'
blood_anno.loc[blood_anno.index[~fil],'myeloid_lymphoid']='Lymphoid'
blood_anno.head()

#make long version of achilles and add annotation
tf_ess = tf_ess.loc[blood_anno.index]
tf_ess['Cell'] = tf_ess.index
tf_ess = tf_ess.melt(id_vars='Cell',value_vars = tf_ess.columns[:-1])
tf_ess.columns=['cell','feature','essentiality']
tf_ess["MorL"] = blood_anno.loc[tf_ess['cell'].values,'myeloid_lymphoid'].values

### we create a similar annotation file for blood cells in HCA data
pbmc_anno = pd.Series({'B cells':'Lymphoid',
          'CD14+ Monocytes':'Myeloid',
          'CD4 T cells':'Lymphoid',
          'CD8 T cells':'Lymphoid',
          'Dendritic cells':'Myeloid',
          'FCGR3A+ Monocytes':'Myeloid',
          'NK cells':'Lymphoid'})
          
#add annotation to hca data
hca_tf_activity['MorL']=pbmc_anno[hca_tf_activity['celltype'].values].values

hca_tf_activity.to_csv('output/hca_data/blood_cancer/tf_activities.csv',sep=',')
tf_ess.to_csv('output/hca_data/blood_cancer/tf_essentialities.csv', sep=',')
```

#### Global analysis of difference in effect size of activity and essentiality
```{r}
tf_act = read_csv("output/hca_data/blood_cancer/tf_activities.csv") %>%
  select(-X1)
tf_ess = read_csv("output/hca_data/blood_cancer/tf_essentialities.csv") %>%
  select(-X1)

a = tf_act %>%
  group_by(technology, tool, feature, confidence, MorL) %>%
  summarize(mean_m = mean(m)) %>%
  ungroup() %>%
  spread(MorL, mean_m) %>%
  mutate(diff_activity = Lymphoid - Myeloid) %>%
  select(-Lymphoid, -Myeloid)

e = tf_ess %>%
  group_by(feature, MorL) %>%
  summarise(mean_ess = mean(essentiality, na.rm=T)) %>%
  spread(MorL, mean_ess) %>%
  mutate(diff_ess = Lymphoid - Myeloid) %>%
  select(-Lymphoid, -Myeloid)

res = a %>%
  inner_join(e, by="feature") %>%
  filter(confidence %in% c("A", "B") | is.na(confidence)) %>%
  nest(-c(technology, tool)) %>%
  mutate(r = data %>% map(function(data) {
    data %>%
      select_if(is.numeric) %>%
      cor.test(.$diff_activity, .$diff_ess, data=.) %>%
      tidy() %>%
      select(estimate, p.value, statistic)
  }))

saveRDS(res, "output/hca_data/blood_cancer/effect_size_diff.rds")
```