in_silico_benchmark.Rmd

---
title: "Application and benchmark of pathway and TF activity inference tools on simulated scRNA-seq pathway and TF perturbation data"
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(scran)
library(GSEABase)
library(dplyr)
library(tidyr)
library(purrr)
library(readr)
library(tibble)
library(ggplot2)
library(stringr)
library(forcats)
library(yardstick)
library(tidylog)
library(furrr)
library(edgeR)
library(viper)
library(Matrix)
library(ggrepel)
library(Rtsne)
library(cluster)
library(biobroom)
library(broom)
library(pheatmap)
library(cluster)
library(umap)
library(AUCell)


source("src/simulation.R")
source("src/my_ggplot_themes.R")

theme_set(theme_cowplot())
options("tidylog.display" = list(print))
options(future.globals.maxSize= 2097152000) # corresponds to 2 GB

plan(multiprocess)
```

### Constructing required meta files
#### Parse benchmark meta data from ChEA3 and manual collection of pathway perturbation
```{r "parse benchmark meta data from ChEA3 and manual collection of pw pert"}
# This file is downloaded from https://amp.pharm.mssm.edu/ChEA3
df = read_delim("data/in_silico_benchmark/meta_data/Single-TF_Perturbations.gmt", 
                delim="\t", col_names = F) %>%
  select(id = X1)

# load DoRothEA TFs
dorothea_tfs = read_csv(
  "data/regulons/dorothea/dorothea_regulon_human_v1.csv"
  ) %>%
  distinct(tf, confidence) %>%
  rename(target = tf)

# define aliases for knockdown and overexpression
knockdown_alias = c("SIRNA", "KO", "SHRNA", "KD")
overexpression_alias = c("ACTIVATION", "OE", "HYPERACTIVE")

# parse the metadata of TF perturbation experiments from ChEA3
dorothea_meta_df = df %>%
  filter(str_detect(id, "RNASEQ")) %>%
  pull(id) %>%
  str_match(
    "^([0-9A-Z]*)_([A-Z]*)_([0-9A-Z]*)_HUMAN_(GSE\\d*)_([0-9A-Z]*).*RNASEQ"
    ) %>%
  as.data.frame() %>%
  setNames(
    ., c("description", "target", "perturbation", "molecule", "accession", "other")
    ) %>%
  na_if("") %>%
  as_tibble() %>%
  mutate(perturbation = case_when(perturbation %in% knockdown_alias ~ "inhibition",
                   perturbation %in% overexpression_alias ~ "activation")) %>%
  mutate(id = paste0("archs:", str_pad(row_number(),4, pad=0))) %>%
  mutate(class = "single_tf_perturbation") %>%
  inner_join(dorothea_tfs)

# load metadata from pathway perturbation experiments
progeny_meta_df = read_csv("data/in_silico_benchmark/meta_data/single_pw_perturbations.csv") %>%
  mutate_if(is.logical, as.character)

# save metadata of pathway and TF perturbation experiments
meta_df = bind_rows(dorothea_meta_df, progeny_meta_df)

saveRDS(meta_df, "output/in_silico_benchmark/meta_data/meta_df.rds")
```

#### Unzip and rename achives manually downloaded from ARCHS4
```{r "unzip and rename achives"}
# load count matrices of all perturbation experiments. Those experiments were 
# manually downloaded from ARCHS4 (https://amp.pharm.mssm.edu/archs4/)
zips = list.files("data/in_silico_benchmark/raw_count_data", full.names = T, 
                  pattern = "*zip") %>%
  map(function(archive) {
    accession = str_match(archive, "GSE\\d*")
    path = file.path("output", "in_silico_benchmark", "raw_count_data")
    unzip(archive, exdir = path)
    file.rename(paste0(path, "/gene_count.tsv"),
                paste0(path, "/", accession, ".tsv"))
  })
```

#### Setup simulation parameters
```{r "setup simulation parameters"}
# setup parameter combinations
library_sizes = c(1000, 2000, 5000, 10000, 20000)
num_cells = c(1, 10, 20, 30, 50, 100)
# libsize : mean of library size distribution for single cells
# num_cells: number of cells simulated from single bulk sample
# reps: how often the simulation should be repeated
param_combi = expand.grid(library_sizes, num_cells) %>%
  as_tibble() %>%
  rename(libsize = Var1, num_cells = Var2) %>%
  mutate(reps = 30)

saveRDS(param_combi, "output/in_silico_benchmark/simulation/paramter_combinations.rds")
```

#### Tidy manually curated sample annotation file
```{r "tidy manually curated sample annotation file"}
param_combi = readRDS("output/in_silico_benchmark/simulation/paramter_combinations.rds")

# This file contains manually curated annotation of samples (control or perturbation sample)
raw_anno = read_delim("data/in_silico_benchmark/meta_data/sample_annotations.csv", delim = ";") %>%
  gather(group, sample, -id) %>%
  separate_rows(sample, sep="[|]") %>%
  arrange(id, sample) %>%
  drop_na(sample)

anno_bulk = raw_anno %>%
  mutate(type = "bulk_sample")

# make sample annotation file for the simulated single cells genereated with 
# different combination of simulation parameters
anno_sc = map_df(unique(param_combi$num_cells), .f = function(i) {
  raw_anno %>%
    mutate(n = i) %>% 
    uncount(n) %>%
    group_by(id, sample) %>%
    mutate(cell_num = 1:i) %>%
    ungroup() %>%
    transmute(id, group, sample, 
              cell = str_c(sample, "_cell_", str_pad(cell_num, nchar(i)+1, pad=0))) %>%
    mutate(type = "single_cells") %>%
    mutate(num_cells = i)
})

bind_rows(anno_bulk, anno_sc) %>%
  saveRDS("output/in_silico_benchmark/meta_data/tidy_sample_annotations.rds")
```

#### Load count data, map samples to perturbation experiments and select only protein coding genes
```{r "load count data and sample samples"}
# load sample annotations
anno = readRDS("output/in_silico_benchmark/meta_data/tidy_sample_annotations.rds") %>%
  filter(type == "bulk_sample")
# load metadata of perturbation experiments and combine it with sample annotations
meta = readRDS("output/in_silico_benchmark/meta_data/meta_df.rds")
meta_df = inner_join(meta, anno, by="id") %>%
  select(id, accession, group, sample)
# list of protein coding genes
ensembl = biomaRt::useMart("ensembl",dataset="hsapiens_gene_ensembl")
pcg = biomaRt::getBM(attributes = c("hgnc_symbol","gene_biotype"), mart = ensembl) %>%
  drop_na() %>%
  filter(gene_biotype == "protein_coding") %>%
  distinct(gene = hgnc_symbol)

# read all count matrices and convert them to tidy format
tidy_counts = list.files("output/in_silico_benchmark/raw_count_data", pattern="*tsv", full.names = T) %>%
  map_df(function(path) {
    accession = path %>%
      basename() %>%
      tools::file_path_sans_ext()
    
    read_delim(path, delim="\t") %>%
      rename(gene = X1) %>%
      gather(sample, count, -gene) %>%
      mutate(accession = accession)
  }) %>%
  inner_join(meta_df, by = c("sample", "accession")) %>%
  select(id, sample, group, gene, count) %>%
  semi_join(pcg, by="gene") %>%
  mutate(gene = as_factor(gene))

# double check whether the mapping of samples is correct
a = meta_df %>% 
  count(id, group) %>%
  mutate_if(is.factor, as.character)
b = tidy_counts %>% 
  distinct(id,  group, sample) %>% 
  count(id, group) %>%
  mutate_if(is.factor, as.character)

stopifnot(identical(a,b))

tidy_counts %>%
  select(-group) %>%
  saveRDS("output/in_silico_benchmark/expression_data/mapped_count_data_tidy.rds")
```

#### Conversion from counts to TPM
```{r "Conversion from counts to TPM}
# extract gene length
ensembl = biomaRt::useMart("ensembl",dataset="hsapiens_gene_ensembl")
gene_length = biomaRt::getBM(attributes = c("hgnc_symbol","transcript_length"), mart = ensembl) %>%
  na_if("") %>%
  drop_na() %>%
  rename(gene = hgnc_symbol, length = transcript_length) %>%
  group_by(gene) %>%
  summarise(length = median(length))

# load count data
counts = readRDS("output/in_silico_benchmark/expression_data/mapped_count_data_tidy.rds")

# convert counts to TPM
tpm_df = counts %>%
  #semi_join(exps, by="id") %>%
  inner_join(gene_length, by="gene") %>%
  # divide by gene lengths in kb
  mutate(tmp = count / (length/1000)) %>%
  group_by(sample) %>%
  mutate(lib_size = sum(tmp)) %>%
  ungroup() %>%
  # divide by library size and one million
  mutate(tpm = tmp / (lib_size/1000000)) %>%
  select(id, sample, gene, tpm) %>%
  mutate(gene = as_factor(gene))

saveRDS(tpm_df, "output/in_silico_benchmark/expression_data/mapped_tpm_data_tidy.rds")
```

#### Create design matrices required for bulk RNAseq normalization and limma
```{r "create design matrices required for bulk rnaseq normalization"}
anno = readRDS("output/in_silico_benchmark/meta_data/tidy_sample_annotations.rds") %>%
  mutate(group = as_factor(group))

design_df = anno %>% 
  nest(-c(id, type, num_cells)) %>%
  transmute(id, type, num_cells, design = pmap(., .f=function(data, type, ...) {
    if (type == "bulk_sample") {
      design_matrix = data %>%
        arrange(sample) %>%
        model.matrix(~0+group, data=.)
    
      rownames(design_matrix) = data$sample
      colnames(design_matrix) = levels(data$group)
    } else if (type == "single_cells") {
      design_matrix = data %>%
        arrange(cell) %>%
        model.matrix(~0+group, data=.)
    
      rownames(design_matrix) = data$cell
      colnames(design_matrix) = levels(data$group)
    }
    return(design_matrix)
  }))

saveRDS(design_df, "output/in_silico_benchmark/meta_data/design_df.rds")
```

### Normalization
#### TF perturbation experiments
```{r "tf perturbation experiments - dorothea"}
# extract only tf-perturbation experiments
meta = readRDS("output/in_silico_benchmark/meta_data/meta_df.rds") %>%
  filter(class == "single_tf_perturbation")

# load bulk expression data
tpm_df = readRDS("output/in_silico_benchmark/expression_data/mapped_tpm_data_tidy.rds") %>%
  nest(-c(id, sample), .key="bulk_sample") %>%
  semi_join(meta, by="id")

count_df = readRDS("output/in_silico_benchmark/expression_data/mapped_count_data_tidy.rds") %>%
  nest(-c(id, sample), .key="raw_counts") %>%
  semi_join(meta, by="id")

# load simulation parameters
param_combi = readRDS("output/in_silico_benchmark/simulation/paramter_combinations.rds") %>%
  mutate(reps = 10) %>%
  arrange(-libsize, -num_cells)

# files required for the pipeline
design_df = readRDS("output/in_silico_benchmark/meta_data/design_df.rds")

# pipeline for simulated single cells
# set.seed 
set.seed(123)
param_combi %>%
  # loop over param combinations
  mutate(col = pmap(., .f = function(libsize, num_cells, reps, ...) {
    bulk_df = tpm_df %>%
        mutate(libsize = libsize, num_cells = num_cells, reps = reps)
    # loop over number of simulation replicates
    sim_res = map_dfr(seq_len(reps), function(run) {
      # single cell simulation
      single_sim_run = bulk_df %>% 
        # loop over single bulk samples with the given simulation setup
        mutate(run = run,
               type = "single_cells",
               raw_counts = future_pmap(., .f = sc_simulation, .progress = T)) %>%
        select(-bulk_sample) %>%
        merge_single_cells()
      
      # save raw counts
      save_path = construct_save_path("dorothea", "raw_counts", "sc",num_cells, libsize, run)
      saveRDS(single_sim_run, save_path)
      
      # normalization 
      if (!exists("design_df")) {
        design_df = readRDS("output/in_silico_benchmark/meta_data/design_df.rds")
      }
      norm_single_sim_run = normalize_sc_or_bulk(single_sim_run, design_df)
      save_path = construct_save_path("dorothea", "norm_expression", "sc",num_cells, libsize, run)
      saveRDS(norm_single_sim_run, save_path)
      rm(single_sim_run)
      
      return(as_tibble(NULL))
    })
  }))

# pipeline for bulk
raw_count_bulk = count_df %>%
  mutate(libsize = NA_integer_, num_cells = NA_integer_, reps = NA_integer_, run=NA_integer_,
         type = "bulk_sample") %>%
  unnest(raw_counts) %>%
  nest(sample, gene, count) %>%
  mutate(raw_counts = data %>% map(function(i) {
    i %>%
      spread(sample, count) %>%
      data.frame(row.names=1, check.names = F, stringsAsFactors = F) %>%
      as.matrix()
  })) %>%
  select(-data)

save_path = construct_save_path("dorothea", "raw_counts", "bulk", NULL, NULL, NULL)
saveRDS(raw_count_bulk, save_path)

# normliatzion
norm_expression_bulk = normalize_sc_or_bulk(raw_count_bulk, design_df)
save_path = construct_save_path("dorothea", "norm_expression", "bulk", NULL, NULL, NULL)
saveRDS(norm_expression_bulk, save_path)
```
#### Pathway perturbation experiments
```{r "normalization - pathway perturbation experiments}
### DO NOT RUN LOCALLY
# extract only pathway-perturbation experiments
meta = readRDS("output/in_silico_benchmark/meta_data/meta_df.rds") %>%
  filter(class == "single_pw_perturbation")

# load bulk expression data
tpm_df = readRDS("output/in_silico_benchmark/expression_data/mapped_tpm_data_tidy.rds") %>%
  nest(-c(id, sample), .key="bulk_sample") %>%
  semi_join(meta, by="id")

count_df = readRDS("output/in_silico_benchmark/expression_data/mapped_count_data_tidy.rds") %>%
  nest(-c(id, sample), .key="raw_counts") %>%
  semi_join(meta, by="id")

# load simulation parameters
param_combi = readRDS("output/in_silico_benchmark/simulation/paramter_combinations.rds") %>%
  mutate(reps = 10) %>%
  arrange(-libsize, -num_cells) %>%
  filter(libsize == 1000 & num_cells == 10) %>%
  mutate(reps = 3)

# files required for the pipeline
design_df = readRDS("output/in_silico_benchmark/meta_data/design_df.rds")


# pipeline for simulated single cells
# set.seed 
set.seed(123)
param_combi %>%
  # loop over param combinations
  mutate(col = pmap(., .f = function(libsize, num_cells, reps, ...) {
    bulk_df = tpm_df %>%
        mutate(libsize = libsize, num_cells = num_cells, reps = reps)
    # loop over number of simulation replicates
    sim_res = map_dfr(seq_len(reps), function(run) {
      # single cell simulation
      single_sim_run = bulk_df %>% 
        # loop over single bulk samples with the given simulation setup
        mutate(run = run,
               type = "single_cells",
               raw_counts = future_pmap(., .f = sc_simulation, .progress = T)) %>%
        select(-bulk_sample) %>%
        merge_single_cells()
      
      # save raw counts
      save_path = construct_save_path("progeny", "raw_counts", "sc",num_cells, libsize, run)
      saveRDS(single_sim_run, save_path)
      
      # normalization 
      if (!exists("design_df")) {
        design_df = readRDS("output/in_silico_benchmark/meta_data/design_df.rds")
      }
      norm_single_sim_run = normalize_sc_or_bulk(single_sim_run, design_df)
      save_path = construct_save_path("progeny", "norm_expression", "sc",num_cells, libsize, run)
      saveRDS(norm_single_sim_run, save_path)
      rm(single_sim_run)
      
      return(as_tibble(NULL))
    })
  }))

# pipeline for bulk
raw_count_bulk = count_df %>%
  mutate(libsize = NA_integer_, num_cells = NA_integer_, reps = NA_integer_, run=NA_integer_,
         type = "bulk_sample") %>%
  unnest(raw_counts) %>%
  nest(sample, gene, count) %>%
  mutate(raw_counts = data %>% map(function(i) {
    i %>%
      spread(sample, count) %>%
      data.frame(row.names=1, check.names = F, stringsAsFactors = F) %>%
      as.matrix()
  })) %>%
  select(-data)

save_path = construct_save_path("progeny", "raw_counts", "bulk", NULL, NULL, NULL)
saveRDS(raw_count_bulk, save_path)

# normliatzion
norm_expression_bulk = normalize_sc_or_bulk(raw_count_bulk, design_df)
save_path = construct_save_path("progeny", "norm_expression", "bulk", NULL, NULL, NULL)
saveRDS(norm_expression_bulk, save_path)
```

### 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)

meta = readRDS("output/in_silico_benchmark/meta_data/meta_df.rds") %>%
  filter(class == "single_tf_perturbation") %>%
  rename(tf = target) %>%
  distinct(tf)

sub_m_dorothea = dorothea %>%
  semi_join(meta)

sub_m_gtex = gtex_regulons %>%
  semi_join(meta)

# those tfs are in dorothea but not in gtex -> need to be removed for in silico benchmark
anti_join(sub_m_dorothea, sub_m_gtex) %>%
  saveRDS("output/in_silico_benchmark/meta_data/missing_tfs_dorothea_vs_gtex.rds")
```


### Single sample analysis
#### VIPER
```{r "single sample analysis - viper}
# DO NOT RUN LOCALLY
# only for those tfs there are perturbation experiments
target_tfs = readRDS("output/in_silico_benchmark/meta_data/meta_df.rds") %>%
  filter(class == "single_tf_perturbation") %>%
  rename(tf = target) %>%
  distinct(tf)

# we need to remove these tfs as they are not available in gtex regulons
missing_tfs = readRDS("output/in_silico_benchmark/meta_data/missing_tfs_dorothea_vs_gtex.rds")

# load dorothea regulons and subset them to targeted and common tfs.
human_dorothea_viper_format = read_csv("data/regulons/dorothea/dorothea_regulon_human_v1.csv") %>%
  semi_join(target_tfs, by="tf") %>%
  anti_join(missing_tfs, by="tf") %>%
  df2regulon()

# load design matrix
design_df = readRDS("output/in_silico_benchmark/meta_data/design_df.rds")

# load normalized expression data, apply viper on each individual sample/cell 
# and build contrast on TF activity level for each perturbatione experiment
list.files("output/in_silico_benchmark/simulation/dorothea/norm_expression", full.names = T) %>%
  future_map(possibly(function(path) {
    target_path = str_replace_all(path, "norm_expression", "ss_viper_contrast")
    if (file.exists(target_path)) return(as_tibble(NULL))
    df = readRDS(path)
    
    
    ss_viper_scores = df %>%
      mutate(
        ss_viper_result = emat %>% map(
          ~viper(
            eset = as.matrix(.x),
            regulon = human_dorothea_viper_format,
            nes = T, 
            method = "scale",
            minsize = 4,
            eset.filter = F, 
            adaptive.size = F))) %>%
      select(-emat)
    
    save_path = str_replace_all(path, "norm_expression", "ss_viper_scores")
    if (!dir.exists(dirname(save_path))) dir.create(dirname(save_path))
    saveRDS(ss_viper_scores, save_path)
    
    ss_viper_contrast = ss_calc_contrast(ss_viper_scores, design_df, "ss_viper_result")
    save_path = str_replace_all(path, "norm_expression", "ss_viper_contrast")
    if (!dir.exists(dirname(save_path))) dir.create(dirname(save_path))
    saveRDS(ss_viper_contrast, save_path)
    
    rm(ss_viper_scores, ss_viper_contrast, df)
    
    return(as_tibble(NULL))
  }, otherwise = as_tibble(NULL), quiet=F), .progress = T)
```
#### metaVIPER
```{r "single sample analysis - metaviper}
# DO NOT RUN LOCALLY
job_id = as.numeric(commandArgs(trailingOnly = T)[1])
# only for those tfs there are perturbation experiments
target_tfs = readRDS("output/in_silico_benchmark/meta_data/meta_df.rds") %>%
  filter(class == "single_tf_perturbation") %>%
  rename(tf = target) %>%
  distinct(tf) %>%
  pull(tf)

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% target_tfs)
    sub_tissue_regulon = tissue_regulon[keep_ix]
    return(sub_tissue_regulon)
  })
design_df = readRDS("output/in_silico_benchmark/meta_data/design_df.rds")

# extract system path pointing to normalized expression data
all_files = list.files("output/in_silico_benchmark/simulation/dorothea/norm_expression", full.names = T)

all_files[((job_id - 1) * 16 + 1) : (job_id * 16)] %>%
  future_map(possibly(function(path) {
    if (file.exists(target_path)) return(as_tibble(NULL))
    df = readRDS(path)
    
    ss_metaviper_scores = df %>%
      mutate(
        ss_metaviper_result = emat %>% map(
          ~viper(
            eset = as.matrix(.x),
            regulon = gtex_regulons,
            nes = T, 
            method = "scale",
            minsize = 4,
            eset.filter = F, 
            adaptive.size = F))) %>%
      select(-emat)
    
    save_path = str_replace_all(path, "norm_expression", "ss_metaviper_scores")
    if (!dir.exists(dirname(save_path))) dir.create(dirname(save_path))
    saveRDS(ss_metaviper_scores, save_path)
    
    ss_metaviper_contrast = ss_calc_contrast(ss_metaviper_scores, design_df, "ss_metaviper_result")
    save_path = str_replace_all(path, "norm_expression", "ss_metaviper_contrast")
    if (!dir.exists(dirname(save_path))) dir.create(dirname(save_path))
    saveRDS(ss_metaviper_contrast, save_path)
    
    rm(ss_metaviper_scores, ss_metaviper_contrast, df)
    
    return(as_tibble(NULL))
  }, otherwise = as_tibble(NULL), quiet = F), .progress = T)
```

#### PROGENy
```{r "single sample analysis - progeny"}
# DO NOT RUN LOCALLY
# files required for the pipeline
design_df = readRDS("output/in_silico_benchmark/meta_data/design_df.rds")
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)

# all 
m_all = model %>%
  select(gene, pathway, weight=zscore)

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


list.files("output/in_silico_benchmark/simulation/progeny/norm_expression", full.names = T, pattern="bulk") %>%
  future_map_dfr(function(path) {
    print(path)
    df = readRDS(path)

    # calcuation of single cell (sample) progeny scores
    ss_progeny_scores = df %>%
      mutate(ss_progeny_scores = pmap(., ss_calc_progeny_scores)) %>%
      unnest(ss_progeny_scores)
    
    save_path = str_replace_all(path, "norm_expression", "ss_progeny_scores")
    if (!dir.exists(dirname(save_path))) dir.create(dirname(save_path))
    saveRDS(ss_progeny_scores, save_path)
    
    ss_progeny_contrast = ss_calc_contrast(ss_progeny_scores, design_df, "ss_progeny_result")
    save_path = str_replace_all(path, "norm_expression", "ss_progeny_contrast")
    if (!dir.exists(dirname(save_path))) dir.create(dirname(save_path))
    saveRDS(ss_progeny_contrast, save_path)
    
    rm(ss_progeny_scores, ss_progeny_contrast)
    
    return(as_tibble(NULL))
  }, .progress = T)
```
#### AUCell - DoRothEA
```{r "single sample analysis d-aucell"}
# only for those tfs there are perturbation experiments
target_tfs = readRDS("output/in_silico_benchmark/meta_data/meta_df.rds") %>%
  filter(class == "single_tf_perturbation") %>%
  rename(tf = target) %>%
  distinct(tf)

# we need to remove this tf as it is not available in gtex regulons
missing_tfs = readRDS("output/in_silico_benchmark/meta_data/missing_tfs_dorothea_vs_gtex.rds")

# subset dorothea
human_dorothea = read_csv("data/regulons/dorothea/dorothea_regulon_human_v1.csv") %>%
  semi_join(target_tfs, by="tf") %>%
  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()


design_df = readRDS("output/in_silico_benchmark/meta_data/design_df.rds")

list.files("output/in_silico_benchmark/simulation/dorothea/norm_expression", full.names = T)[1] %>%
  future_map(possibly(function(path) {
    df = readRDS(path)
    
    ss_daucell_scores = df %>%
      mutate(ss_daucell_result = emat %>% map(function(emat) {
        # g = subsetGeneSets(genesets, rownames(emat)) 
        obj = AUCell_buildRankings(emat, nCores=1, plotStats = F, verbose = F) %>%
          AUCell_calcAUC(genesets, ., verbose=F)
        
        res = obj@assays$data@listData$AUC
      })) %>%
      select(-emat)
      
    save_path = str_replace_all(path, "norm_expression", "ss_daucell_scores")
    if (!dir.exists(dirname(save_path))) dir.create(dirname(save_path))
    saveRDS(ss_daucell_scores, save_path)
    
    ss_daucell_contrast = ss_calc_contrast(ss_daucell_scores, design_df, "ss_daucell_result")
    save_path = str_replace_all(path, "norm_expression", "ss_daucell_contrast")
    if (!dir.exists(dirname(save_path))) dir.create(dirname(save_path))
    saveRDS(ss_daucell_contrast, save_path)
    
    rm(ss_daucell_scores, ss_daucell_contrast, df)
    
    return(as_tibble(NULL))
  },otherwise = as_tibble(NULL), quiet = F), .progress = T)
```

#### AUCell - PROGENy
```{r "single sample analysis p-aucell}
# files required for the pipeline
design_df = readRDS("output/in_silico_benchmark/meta_data/design_df.rds")
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))

design_df = readRDS("output/in_silico_benchmark/meta_data/design_df.rds")

list.files("output/in_silico_benchmark/simulation/progeny/norm_expression", full.names = T) %>%
  future_map(possibly(function(path) {
    
    df = readRDS(path) 
    
    input = models %>%
      mutate(emat = list(df)) %>%
      unnest(emat, .preserve = genesets)
    
    
    ss_paucell_scores = input %>%
      mutate(ss_paucell_result = pmap(., .f=function(emat, genesets, ...) {
        
        g = subsetGeneSets(genesets, rownames(emat))
        obj = AUCell_buildRankings(emat, nCores=1, plotStats = F, verbose = F) %>%
          AUCell_calcAUC(genesets, ., verbose=F)
        
        res = obj@assays$data@listData$AUC
      })) %>%
      select(-emat, -genesets)
      
    save_path = str_replace_all(path, "norm_expression", "ss_paucell_scores")
    
    if (!dir.exists(dirname(save_path))) dir.create(dirname(save_path))
    saveRDS(ss_paucell_scores, save_path)
    
    ss_paucell_contrast = ss_calc_contrast(ss_paucell_scores, design_df, "ss_paucell_result")
    save_path = str_replace_all(path, "norm_expression", "ss_paucell_contrast")
    if (!dir.exists(dirname(save_path))) dir.create(dirname(save_path))
    saveRDS(ss_paucell_contrast, save_path)
    
    rm(ss_paucell_scores, ss_paucell_contrast)
    
    return(as_tibble(NULL))
  }, otherwise = as_tibble(NULL), quiet = F), .progress = T)
```


### Performance evaluation
#### VIPER
```{r "performance evaluation - viper"}
meta = readRDS("output/in_silico_benchmark/meta_data/meta_df.rds") %>%
  select(target, perturbation, id)

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

basic_setup = tribble(
  ~confidence, ~group,
  # "A", "single_conf",
  # "B", "single_conf",
  # "C", "single_conf",
  # "D", "single_conf",
  # "E", "single_conf",
  "A", "multi_conf",
  "AB", "multi_conf",
  "ABC", "multi_conf",
  "ABCD", "multi_conf",
  "ABCDE", "multi_conf"
) %>%
    mutate(confidence = factor(confidence, levels=c("A", "B", "C", "D", "E",
                                                    "AB", "ABC", "ABCD", 
                                                    "ABCDE")))

# path = list.files("output/in_silico_benchmark/simulation/dorothea/ss_viper_contrast", full.names = T)[4]
performance = list.files("output/in_silico_benchmark/simulation/dorothea/ss_viper_contrast", full.names = T) %>%
  future_map_dfr(function(path) {
    df = readRDS(path) %>%
      unnest() %>%
      rename(tf = gene, activity = logfc) %>%
      left_join(human_dorothea, by="tf") %>%
      left_join(meta)
    
    setup = basic_setup %>%
      mutate(viper_scores = list(df))
    
    input_performance = setup %>%
      mutate(input = pmap(., .f = function(confidence, viper_scores, ...) {
          conf_levels = str_split(confidence, "") %>% 
            pluck(1)
          
          # extract number of uniquely covered tfs
          coverage = viper_scores %>% 
            filter(confidence %in% conf_levels) %>%
            filter(tf == target) %>%
            summarise(coverage = n_distinct(tf)) %>%
            pull()
          
          viper_scores %>%
            filter(confidence %in% conf_levels) %>%
            mutate(response = case_when(tf == target ~ 1,
                                        TRUE ~ 0),
                   response = factor(response, levels = c("1", "0")),
                   predictor = case_when(perturbation == "inhibition" ~ -1 * activity,
                                         perturbation == "activation" ~ activity),
                   coverage = coverage) %>%
            select(libsize, num_cells, reps, run, type, response, predictor, coverage)
      })) %>%
      select(-viper_scores) %>%
      unnest(input)
    
    meta_performance = input_performance %>%
      dplyr::count(group, confidence, libsize, num_cells, reps, run, type, response, coverage) %>%
      spread(response, n) %>%
      dplyr::rename(positive = `1`, negative = `0`) %>%
      mutate(random = positive/(positive + negative))
    
    # acess metric coordinates
    roc_coords = input_performance %>%
      group_by(confidence, group, libsize, num_cells, reps, run, type) %>%
      roc_curve(response, predictor, options=list(direction = "<", transpose = F)) %>%
      ungroup() %>%
      mutate(metric = "roc_curve") %>%
      left_join(meta_performance, by=c("group", "confidence", "libsize", "num_cells", "reps", "run", "type")) %>%
      mutate(random = 0.5)
    
    roc_path = path %>% 
      str_replace("ss_viper_contrast", "ss_viper_performance") %>%
      str_replace("ss_viper_contrast", "roc_curve")
    saveRDS(roc_coords, roc_path)
    
    pr_coords = input_performance %>%
      group_by(confidence, group, libsize, num_cells, reps, run, type) %>%
      pr_curve(response, predictor) %>%
      ungroup() %>%
      mutate(metric = "pr_curve") %>%
      left_join(meta_performance, by=c("group", "confidence", "libsize", "num_cells", "reps", "run", "type"))
    
    pr_path = path %>% 
      str_replace("ss_viper_contrast", "ss_viper_performance") %>%
      str_replace("ss_viper_contrast", "pr_curve")
    saveRDS(pr_coords, pr_path)
    
    
    
    # acess summary metrics
    performance_result_roc = input_performance %>% 
      group_by(confidence, group, libsize, num_cells, reps, run, type) %>%
      roc_auc(response, predictor, options=list(direction = "<", transpose=T))
    
    performance_result_pr = input_performance %>%
      group_by(confidence, group, libsize, num_cells, reps, run, type) %>%
      pr_auc(response, predictor)
      
    
    perfomance_result = bind_rows(performance_result_roc, performance_result_pr) %>%
      left_join(meta_performance, by = c("confidence", "group", "libsize", "num_cells", "reps", "run", "type")) %>%
      rename(metric = .metric, auc = .estimate) %>%
      select(-.estimator) %>%
      mutate(random = case_when(metric == "roc_auc" ~ 0.5,
                                TRUE ~  random),
             label = case_when(type == "bulk_sample" ~ str_c("bulk", confidence, sep = "_"),
                               type == "single_cells" ~ str_c("sc", confidence, sep="_"))) %>%
      nest(-c(libsize, num_cells, reps, run, type))
      }, .progress = T)

saveRDS(performance, "output/in_silico_benchmark/simulation/dorothea/ss_viper_performance/performance_result.rds")
```

#### PROGENy
```{r "performance evaluation - progeny"}
meta = readRDS("output/in_silico_benchmark/meta_data/meta_df.rds") %>%
  select(target, perturbation, id)

# path = list.files("output/in_silico_benchmark/simulation/progeny/ss_progeny_contrast", full.names = T)[1]
performance = list.files("output/in_silico_benchmark/simulation/progeny/ss_progeny_contrast", full.names = T) %>%
  future_map_dfr(function(path) {
    print(path)
    df = readRDS(path) %>%
      unnest() %>%
      rename(pathway = gene, activity = logfc) %>%
      left_join(meta)
    
    setup = df %>%
      nest(-footprints, .key="progeny_scores")
    
    input_performance = setup %>%
      mutate(input = pmap(., .f = function(progeny_scores, footprints, ...) {
          
          # extract number of uniquely covered pathway
          coverage = progeny_scores %>% 
            filter(pathway == target) %>%
            summarise(coverage = n_distinct(pathway)) %>%
            pull()
          
          progeny_scores %>%
            mutate(response = case_when(pathway == target ~ 1,
                                        TRUE ~ 0),
                   response = factor(response, levels = c("1", "0")),
                   predictor = case_when(perturbation == "inhibition" ~ -1 * activity,
                                         perturbation == "activation" ~ activity),
                   coverage = coverage) %>%
            select(libsize, num_cells, reps, run, type, response, predictor, coverage)
      })) %>%
      select(-progeny_scores) %>%
      unnest(input)
    
    meta_performance = input_performance %>%
      count(footprints, libsize, num_cells, reps, run, type, response, coverage) %>%
      spread(response, n) %>%
      rename(positive = `1`, negative = `0`) %>%
      mutate(random = positive/(positive + negative))
    
    # acess metric coordinates
    roc_coords = input_performance %>%
      group_by(footprints, libsize, num_cells, reps, run, type) %>%
      roc_curve(response, predictor, options=list(direction = "<", transpose = F)) %>%
      ungroup() %>%
      mutate(metric = "roc_curve") %>%
      left_join(meta_performance, by=c("footprints", "libsize", "num_cells", "reps", "run", "type")) %>%
      mutate(random = 0.5)
    
    roc_path = path %>% 
      str_replace("ss_progeny_contrast", "ss_progeny_performance") %>%
      str_replace("ss_progeny_contrast", "roc_curve")
    saveRDS(roc_coords, roc_path)
    
    pr_coords = input_performance %>%
      group_by(footprints, libsize, num_cells, reps, run, type) %>%
      pr_curve(response, predictor) %>%
      ungroup() %>%
      mutate(metric = "pr_curve") %>%
      left_join(meta_performance, by=c("footprints", "libsize", "num_cells", "reps", "run", "type"))
    
    pr_path = path %>% 
      str_replace("ss_progeny_contrast", "ss_progeny_performance") %>%
      str_replace("ss_progeny_contrast", "pr_curve")
    saveRDS(pr_coords, pr_path)
    
    performance_result_roc = input_performance %>% 
      group_by(footprints, libsize, num_cells, reps, run, type) %>%
      roc_auc(response, predictor, options=list(direction = "<"))
    
    
    performance_result_pr = input_performance %>%
      group_by(footprints, libsize, num_cells, reps, run, type) %>%
      pr_auc(response, predictor)
      
    
    performance_result = bind_rows(performance_result_roc, performance_result_pr) %>%
      left_join(meta_performance, by = c("footprints", "libsize", "num_cells", "reps", "run", "type")) %>%
      rename(metric = .metric, auc = .estimate) %>%
      select(-.estimator) %>%
      mutate(random = case_when(metric == "roc_auc" ~ 0.5,
                                TRUE ~  random),
             label = case_when(type == "bulk_sample" ~ str_c("bulk", footprints, sep = "_"),
                               type == "single_cells" ~ str_c("sc", footprints, sep="_"))) %>%
      nest(-c(libsize, num_cells, reps, run, type))
      }, .progress = T)

saveRDS(performance, "output/in_silico_benchmark/simulation/progeny/ss_progeny_performance/performance_result.rds")
```

#### AUCell-PROGENy
```{r performance evaluation - p-aucell"}
meta = readRDS("output/in_silico_benchmark/meta_data/meta_df.rds") %>%
  select(target, perturbation, id)

# path = list.files("output/in_silico_benchmark/simulation/progeny/ss_paucell_contrast", full.names = T)[2]
performance = list.files("output/in_silico_benchmark/simulation/progeny/ss_paucell_contrast", full.names = T) %>%
  future_map_dfr(function(path) {
    print(path)
    df = readRDS(path) %>%
      unnest() %>%
      rename(pathway = gene, activity = logfc) %>%
      left_join(meta)
    
    setup = df %>%
      nest(-footprints, .key="progeny_scores")
    
    input_performance = setup %>%
      mutate(input = pmap(., .f = function(progeny_scores, footprints, ...) {
          
          # extract number of uniquely covered pathway
          coverage = progeny_scores %>% 
            filter(pathway == target) %>%
            summarise(coverage = n_distinct(pathway)) %>%
            pull()
          
          progeny_scores %>%
            mutate(response = case_when(pathway == target ~ 1,
                                        TRUE ~ 0),
                   response = factor(response, levels = c("1", "0")),
                   predictor = case_when(perturbation == "inhibition" ~ -1 * activity,
                                         perturbation == "activation" ~ activity),
                   coverage = coverage) %>%
            select(libsize, num_cells, reps, run, type, response, predictor, coverage)
      })) %>%
      select(-progeny_scores) %>%
      unnest(input)
    
    meta_performance = input_performance %>%
      count(footprints, libsize, num_cells, reps, run, type, response, coverage) %>%
      spread(response, n) %>%
      rename(positive = `1`, negative = `0`) %>%
      mutate(random = positive/(positive + negative))
    
    # acess metric coordinates
    roc_coords = input_performance %>%
      group_by(footprints, libsize, num_cells, reps, run, type) %>%
      roc_curve(response, predictor, options=list(direction = "<", transpose = T)) %>%
      ungroup() %>%
      mutate(metric = "roc_curve") %>%
      left_join(meta_performance, by=c("footprints", "libsize", "num_cells", "reps", "run", "type")) %>%
      mutate(random = 0.5)
    
    roc_path = path %>% 
      str_replace("ss_paucell_contrast", "ss_paucell_performance") %>%
      str_replace("ss_paucell_contrast", "roc_curve")
    saveRDS(roc_coords, roc_path)
    
    pr_coords = input_performance %>%
      group_by(footprints, libsize, num_cells, reps, run, type) %>%
      pr_curve(response, predictor) %>%
      ungroup() %>%
      mutate(metric = "pr_curve") %>%
      left_join(meta_performance, by=c("footprints", "libsize", "num_cells", "reps", "run", "type"))
    
    pr_path = path %>% 
      str_replace("ss_paucell_contrast", "ss_paucell_performance") %>%
      str_replace("ss_paucell_contrast", "pr_curve")
    saveRDS(pr_coords, pr_path)
    
    performance_result_roc = input_performance %>% 
      group_by(footprints, libsize, num_cells, reps, run, type) %>%
      roc_auc(response, predictor, options=list(direction = "<", transpose = T))
    
    
    performance_result_pr = input_performance %>%
      group_by(footprints, libsize, num_cells, reps, run, type) %>%
      pr_auc(response, predictor)
      
    
    performance_result = bind_rows(performance_result_roc, performance_result_pr) %>%
      left_join(meta_performance, by = c("footprints", "libsize", "num_cells", "reps", "run", "type")) %>%
      rename(metric = .metric, auc = .estimate) %>%
      select(-.estimator) %>%
      mutate(random = case_when(metric == "roc_auc" ~ 0.5,
                                TRUE ~  random),
             label = case_when(type == "bulk_sample" ~ str_c("bulk", footprints, sep = "_"),
                               type == "single_cells" ~ str_c("sc", footprints, sep="_"))) %>%
      nest(-c(libsize, num_cells, reps, run, type))
      }, .progress = T)

saveRDS(performance, "output/in_silico_benchmark/simulation/progeny/ss_paucell_performance/performance_result.rds")
```

#### AUCell-DoRothEA
```{r "performance evaluation - d-aucell"}
meta = readRDS("output/in_silico_benchmark/meta_data/meta_df.rds") %>%
  select(target, perturbation, id)

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

basic_setup = tribble(
  ~confidence, ~group,
  "A", "multi_conf",
  "AB", "multi_conf",
  "ABC", "multi_conf",
  "ABCD", "multi_conf",
  "ABCDE", "multi_conf"
) %>%
    mutate(confidence = factor(confidence, levels=c("A", "B", "C", "D", "E",
                                                    "AB", "ABC", "ABCD", 
                                                    "ABCDE")))

# path = list.files("output/in_silico_benchmark/simulation/dorothea/ss_daucell_contrast", full.names = T)[4]
performance = list.files("output/in_silico_benchmark/simulation/dorothea/ss_daucell_contrast", full.names = T) %>%
  future_map_dfr(function(path) {
    df = readRDS(path) %>%
      unnest() %>%
      rename(tf = gene, activity = logfc) %>%
      left_join(human_dorothea, by="tf") %>%
      left_join(meta)
    
    setup = basic_setup %>%
      mutate(viper_scores = list(df))
    
    input_performance = setup %>%
      mutate(input = pmap(., .f = function(confidence, viper_scores, ...) {
          conf_levels = str_split(confidence, "") %>% 
            pluck(1)
          
          # extract number of uniquely covered tfs
          coverage = viper_scores %>% 
            filter(confidence %in% conf_levels) %>%
            filter(tf == target) %>%
            summarise(coverage = n_distinct(tf)) %>%
            pull()
          
          viper_scores %>%
            filter(confidence %in% conf_levels) %>%
            mutate(response = case_when(tf == target ~ 1,
                                        TRUE ~ 0),
                   response = factor(response, levels = c("1", "0")),
                   predictor = case_when(perturbation == "inhibition" ~ -1 * activity,
                                         perturbation == "activation" ~ activity),
                   coverage = coverage) %>%
            select(libsize, num_cells, reps, run, type, response, predictor, coverage)
      })) %>%
      select(-viper_scores) %>%
      unnest(input)
    
    meta_performance = input_performance %>%
      dplyr::count(group, confidence, libsize, num_cells, reps, run, type, response, coverage) %>%
      spread(response, n) %>%
      dplyr::rename(positive = `1`, negative = `0`) %>%
      mutate(random = positive/(positive + negative))
    
    # acess metric coordinates
    roc_coords = input_performance %>%
      group_by(confidence, group, libsize, num_cells, reps, run, type) %>%
      roc_curve(response, predictor, options=list(direction = "<", transpose = F)) %>%
      ungroup() %>%
      mutate(metric = "roc_curve") %>%
      left_join(meta_performance, by=c("group", "confidence", "libsize", "num_cells", "reps", "run", "type")) %>%
      mutate(random = 0.5)
    
    roc_path = path %>% 
      str_replace("ss_daucell_contrast", "ss_daucell_performance") %>%
      str_replace("ss_daucell_contrast", "roc_curve")
    saveRDS(roc_coords, roc_path)
    
    pr_coords = input_performance %>%
      group_by(confidence, group, libsize, num_cells, reps, run, type) %>%
      pr_curve(response, predictor) %>%
      ungroup() %>%
      mutate(metric = "pr_curve") %>%
      left_join(meta_performance, by=c("group", "confidence", "libsize", "num_cells", "reps", "run", "type"))
    
    pr_path = path %>% 
      str_replace("ss_daucell_contrast", "ss_daucell_performance") %>%
      str_replace("ss_daucell_contrast", "pr_curve")
    saveRDS(pr_coords, pr_path)
    
    
    
    # acess summary metrics
    performance_result_roc = input_performance %>% 
      group_by(confidence, group, libsize, num_cells, reps, run, type) %>%
      roc_auc(response, predictor, options=list(direction = "<", transpose=T))
    
    performance_result_pr = input_performance %>%
      group_by(confidence, group, libsize, num_cells, reps, run, type) %>%
      pr_auc(response, predictor)
      
    
    perfomance_result = bind_rows(performance_result_roc, performance_result_pr) %>%
      left_join(meta_performance, by = c("confidence", "group", "libsize", "num_cells", "reps", "run", "type")) %>%
      rename(metric = .metric, auc = .estimate) %>%
      select(-.estimator) %>%
      mutate(random = case_when(metric == "roc_auc" ~ 0.5,
                                TRUE ~  random),
             label = case_when(type == "bulk_sample" ~ str_c("bulk", confidence, sep = "_"),
                               type == "single_cells" ~ str_c("sc", confidence, sep="_"))) %>%
      nest(-c(libsize, num_cells, reps, run, type))
      }, .progress = T)

saveRDS(performance, "output/in_silico_benchmark/simulation/dorothea/ss_daucell_performance/performance_result.rds")
```

#### metaVIPER
```{r "performance evaluation - metaviper"}
meta = readRDS("output/in_silico_benchmark/meta_data/meta_df.rds") %>%
  select(target, perturbation, id)

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

basic_setup = tribble(
  ~confidence, ~group,
  "A", "multi_conf",
  "AB", "multi_conf",
  "ABC", "multi_conf",
  "ABCD", "multi_conf",
  "ABCDE", "multi_conf"
) %>%
    mutate(confidence = factor(confidence, levels=c("A", "B", "C", "D", "E",
                                                    "AB", "ABC", "ABCD", 
                                                    "ABCDE")))

# path = list.files("output/in_silico_benchmark/simulation/dorothea/ss_metaviper_contrast", full.names = T)[4]
performance = list.files("output/in_silico_benchmark/simulation/dorothea/ss_metaviper_contrast", full.names = T) %>%
  future_map_dfr(function(path) {
    df = readRDS(path) %>%
      unnest() %>%
      rename(tf = gene, activity = logfc) %>%
      left_join(human_dorothea, by="tf") %>%
      left_join(meta)
    
    setup = basic_setup %>%
      mutate(viper_scores = list(df))

    input_performance = setup %>%
      mutate(input = pmap(., .f = function(confidence, viper_scores, ...) {
          conf_levels = str_split(confidence, "") %>% 
            pluck(1)
          
          # extract number of uniquely covered tfs
          coverage = viper_scores %>% 
            filter(confidence %in% conf_levels) %>%
            filter(tf == target) %>%
            summarise(coverage = n_distinct(tf)) %>%
            pull()
          
          viper_scores %>%
            filter(confidence %in% conf_levels) %>%
            mutate(response = case_when(tf == target ~ 1,
                                        TRUE ~ 0),
                   response = factor(response, levels = c("1", "0")),
                   predictor = case_when(perturbation == "inhibition" ~ -1 * activity,
                                         perturbation == "activation" ~ activity),
                   coverage = coverage) %>%
            select(libsize, num_cells, reps, run, type, response, predictor, coverage)
      })) %>%
      select(-viper_scores) %>%
      unnest(input)
    
    meta_performance = input_performance %>%
      dplyr::count(group, confidence, libsize, num_cells, reps, run, type, response, coverage) %>%
      spread(response, n) %>%
      dplyr::rename(positive = `1`, negative = `0`) %>%
      mutate(random = positive/(positive + negative))
    
    # acess metric coordinates
    roc_coords = input_performance %>%
      group_by(confidence, group, libsize, num_cells, reps, run, type) %>%
      roc_curve(response, predictor, options=list(direction = "<", transpose = F)) %>%
      ungroup() %>%
      mutate(metric = "roc_curve") %>%
      left_join(meta_performance, by=c("group", "confidence", "libsize", "num_cells", "reps", "run", "type")) %>%
      mutate(random = 0.5)
    
    roc_path = path %>% 
      str_replace("ss_metaviper_contrast", "ss_metaviper_performance") %>%
      str_replace("ss_metaviper_contrast", "roc_curve")
    saveRDS(roc_coords, roc_path)
    
    pr_coords = input_performance %>%
      group_by(confidence, group, libsize, num_cells, reps, run, type) %>%
      pr_curve(response, predictor) %>%
      ungroup() %>%
      mutate(metric = "pr_curve") %>%
      left_join(meta_performance, by=c("group", "confidence", "libsize", "num_cells", "reps", "run", "type"))
    
    pr_path = path %>% 
      str_replace("ss_metaviper_contrast", "ss_metaviper_performance") %>%
      str_replace("ss_metaviper_contrast", "pr_curve")
    saveRDS(pr_coords, pr_path)
    
    
    
    # acess summary metrics
    performance_result_roc = input_performance %>% 
      group_by(confidence, group, libsize, num_cells, reps, run, type) %>%
      roc_auc(response, predictor, options=list(direction = "<", transpose=T))
    
    performance_result_pr = input_performance %>%
      group_by(confidence, group, libsize, num_cells, reps, run, type) %>%
      pr_auc(response, predictor)
      
    
    perfomance_result = bind_rows(performance_result_roc, performance_result_pr) %>%
      left_join(meta_performance, by = c("confidence", "group", "libsize", "num_cells", "reps", "run", "type")) %>%
      rename(metric = .metric, auc = .estimate) %>%
      select(-.estimator) %>%
      mutate(random = case_when(metric == "roc_auc" ~ 0.5,
                                TRUE ~  random),
             label = case_when(type == "bulk_sample" ~ str_c("bulk", confidence, sep = "_"),
                               type == "single_cells" ~ str_c("sc", confidence, sep="_"))) %>%
      nest(-c(libsize, num_cells, reps, run, type))
      }, .progress = T)

saveRDS(performance, "output/in_silico_benchmark/simulation/dorothea/ss_metaviper_performance/performance_result.rds")
```