Figure3-rna-fm_seperate.Rmd

---
title: "Figure3 - RNAseq"
author: "E Onur Karakaslar"
date: "1/15/2020"
output: html_document
---

```{r setup, include=FALSE}
require(knitr)
knitr::opts_chunk$set(echo = TRUE)
opts_knit$set(root.dir = "/Users/karako/Dropbox (JAX)/MouseAging_clean/") #set root dir!
```

```{r library}
library(useR)  # for clear function, you can delete this one.
library(edgeR) # finding Differentially Expressed genes
library(limma) # for quantile normalization
library(dplyr) # using pipe (%>%) and select
library(ggpubr)
library(ggplot2) 
library(writexl) 
library(tidyverse)
library(preprocessCore)
```


```{r save_as_RData}
save_as_RData <- function(){
  
  count_matrix <- read.csv("data/RNAseq/F3_input/rna_count_matrix.csv")
  cols <- colnames(count_matrix[,-1]) %>%strsplit(".", fixed = T)
  
  tissues <- sapply(cols, function(x){
    x[4]
  })
  
  bm     <- count_matrix [ , c(1 , which(tissues == "BM") + 1)]
  pbl    <- count_matrix [ , c(1 , which(tissues == "PBL") + 1)]
  naive  <- count_matrix [ , c(1 , which(tissues == "naive") + 1)]
  spleen <- count_matrix [ , c(1 , which(tissues == "spleen") + 1)]
  memory <- count_matrix [ , c(1 , which(tissues == "memory") + 1)]
  
  save(bm    , file = "data/RNAseq/F3_input/rnaseq-expcount-bm.RData")
  save(pbl   , file = "data/RNAseq/F3_input/rnaseq-expcount-pbl.RData")
  save(naive , file = "data/RNAseq/F3_input/rnaseq-expcount-naive.RData")
  save(spleen, file = "data/RNAseq/F3_input/rnaseq-expcount-spleen.RData")
  save(memory, file = "data/RNAseq/F3_input/rnaseq-expcount-memory.RData")
}
```



```{r preprocess-data}
preprocess_data <- function(tissue_cell_type){
    
    cat  (paste0("Loading RNAseq data for: ", toupper(tissue_cell_type), "\n"))
    name <- load (paste0("./data/RNAseq/F3_input/rnaseq-expcount-",tissue_cell_type,'.RData'))
    count.matrix <- get(name)
    gene.names <- count.matrix[,1]
    count.matrix <- count.matrix[,-1]
    
    sample_rna <- colnames(count.matrix) %>% strsplit(".", fixed = T)
    STRAIN = TYPE <- sapply(sample_rna, function(x){
      x[1]
    })
    
    AGE    <- sapply(sample_rna, function(x){
      as.numeric(gsub("([0-9]+).*$", "\\1", x[2] %>% trimws)) 
    }) 
    
    GENDER <- sapply(sample_rna, function(x){
      x[3]
    }) 
    
    TISSUE <- sapply(sample_rna, function(x){
      x[4]
    }) 
    
    SAMPLEMOUSEID <- sapply(sample_rna, function(x){
      x[5]
    }) 
    
    specs = NULL
    specs$STRAIN <- STRAIN
    specs$TYPE   <- TYPE
    specs$AGE    <- AGE
    specs$GENDER <- GENDER 
    specs$TISSUE <- TISSUE
    specs$SAMPLEMOUSEID <- SAMPLEMOUSEID
    
    result <- list(count.matrix = count.matrix,gene.names=gene.names, specs=specs)
    return (result)
}
```

```{r load genesets}
load_genesets <- function(){
  
  load('data/genesets/scRNA_and_DICE/geneset.info.RData')
  
  assign("selected_genesets", value = list(
    scrnaseq_tcells_specific_10x   = geneset.genes.scrnaseq_tcells_specific_10x,
    vp2008                         = geneset.genes.vp2008,
    wp                             = geneset.genes.wp,
    scrnaseq_tcells_expressed_10x  = geneset.genes.scrnaseq_tcells_expressed_10x,
    scrnaseq_pbmc_top              = geneset.genes.scrnaseq_pbmc_top,
    scrnaseq_pbmc_simple_exclusive = geneset.genes.scrnaseq_pbmc_simple_exclusive,
    scrnaseq_pbmc_simple_specific  = geneset.genes.scrnaseq_pbmc_simple_specific,
    #gobp                           = geneset.genes.gobp,
    #gomf                           = geneset.genes.gomf,
    dice_major                     = geneset.genes.dice_major
  ), envir = .GlobalEnv)
  
  selected_genesets_mice <- lapply(selected_genesets, function(gs){
    gs %>% convertHumanGeneList
  })
  assign("selected_genesets_mice", selected_genesets_mice, .GlobalEnv)
  
  assign("selected_genesets_labels", list(
    scrnaseq_tcells_specific_10x   = geneset.names.scrnaseq_tcells_specific_10x,
    vp2008                         = geneset.names.vp2008,
    wp                             = geneset.names.wp,
    scrnaseq_tcells_expressed_10x  = geneset.names.scrnaseq_tcells_expressed_10x,
    scrnaseq_pbmc_top              = geneset.names.scrnaseq_pbmc_top,
    scrnaseq_pbmc_simple_exclusive = geneset.names.scrnaseq_pbmc_simple_exclusive,
    scrnaseq_pbmc_simple_specific  = geneset.names.scrnaseq_pbmc_simple_specific,
    #gobp                           = geneset.names.gobp,
    #gomf                           = geneset.names.gomf,
    dice_major                     = geneset.names.dice_major
  ), envir = .GlobalEnv) 
  
  assign("union_size", 
         lapply(selected_genesets, function(gs){
           gs[,"GeneName"]
         }) %>% unlist(recursive = F) %>% unique %>% length, envir = .GlobalEnv)
    
}
```


```{r convert_genesets}
# Basic function to convert human gene names to mouse ensembl gene ids
convertHumanGeneList <- function(x){

x_genename <- x [,"GeneName"] %>% unique

require("biomaRt")
if (!exists("human")){
  assign( x     = "human", 
          value = useMart("ensembl", dataset = "hsapiens_gene_ensembl"),
          envir = .GlobalEnv)
  cat("Human Genes are imported...\n")
}
if (!exists("mouse")){
  assign( x     = "mouse",
          value = useMart("ensembl", dataset = "mmusculus_gene_ensembl"),
          envir = .GlobalEnv)
  cat("Mouse Genes are imported...\n")
}

# map from human to mice
genesV2 = getLDS(attributes = c("hgnc_symbol"), 
                 filters = "hgnc_symbol", 
                 values = x_genename , 
                 mart = human, 
                 attributesL = c("ensembl_gene_id"), 
                 martL = mouse, uniqueRows=T)
# each human genome should be unique, so 1 to 1 map should be possible
# genesV2 <- genesV2[ !duplicated(genesV2$HGNC.symbol),]
genesV2 <- genesV2[ !duplicated(genesV2$Gene.stable.ID),]

humanx <- merge(x, genesV2, by.x = "GeneName", by.y = "HGNC.symbol")

return(humanx)
}

```


```{r}
#' Differential Expression Analysis
#'
#' This function takes a count matrix, normalizes it with TMM and returns a fit matrix
#' 
#' @param data raw count matrix, rownames should be Gene IDs
#' @param specs specialities of mice: age, gender, strain
#' @return fit matrix which can be used later for differential analysis

DE_fit <- function(count.matrix, specs){
  y       <- DGEList(counts = count.matrix)  
  keep    <- filterByExpr(y, min.count = 1)
  y       <- y[keep,,keep.lib.sizes = F]
  y       <- calcNormFactors(object = y, method = "TMM")
  group   <- factor(paste(specs$GENDER,specs$AGE,specs$TYPE,sep="."))
  design  <- model.matrix(~0+group)
  colnames(design) <- levels(group)
  y       <- estimateDisp(y, design) 
  fit     <- glmQLFit(y, design, robust = TRUE) # recommended in edgeR manual 4.4.7
  return (fit)
} 
```

```{r DE_fit_quantile}

#' Differential Expression Analysis
#'
#' This function takes a count matrix normalizes it with quantile normalization and
#' return fit matrix.
#' 
#' @param data raw count matrix, rownames should be Gene IDs
#' @param specs specialities of mice: age, gender, strain
#' @return fit matrix which can be used later for differential analysis

DE_fit_quantile <- function(count.matrix, specs){
  y       <- DGEList(counts = count.matrix)  
  keep    <- filterByExpr(y, min.count=1)
  y       <- y[keep,,keep.lib.sizes=F]
  y       <- calcNormFactors(object = y, method = "none")
  group   <- factor(paste(specs$GENDER,specs$AGE,specs$TYPE,sep="."))
  design  <- model.matrix(~0+group)
  colnames(design) <- levels(group)
  rownames(design) <- colnames(count.matrix)
  v       <- voom(y,design,plot = F, normalize.method ="quantile")
  fit     <- lmFit(v, design)
  fit$aveLogCPM    <- aveLogCPM(count.matrix)
  return (fit)
} 
```

```{r DE_test}
#' @param fit matrix containing model parameters and design matrix
#' @param contrast DE contrast
#' @return qlm f-test in edgeR (qlf)
DE_test <- function(fit, contrast ="Gender_3.B6"){
  contrasts.age_sex_str <- makeContrasts(

    Gender_3.B6    = (M.3.B6   - F.3.B6),
    Gender_3.NZO   = (M.3.NZO  - F.3.NZO),
    Gender_18.B6    = (M.18.B6  - F.18.B6),
    Gender_18.NZO   = (M.18.NZO - F.18.NZO),

    Age18vs3_M.B6  = (M.18.B6  - M.3.B6),
    Age18vs3_M.NZO = (M.18.NZO - M.3.NZO),
    Age18vs3_F.B6  = (F.18.B6  - F.3.B6),
    Age18vs3_F.NZO = (F.18.NZO - F.3.NZO),
  levels = fit$design)
  
  qlf <- glmQLFTest(fit, contrast = contrasts.age_sex_str[, contrast])
  return(qlf)
}
```

```{r DE_test_quantile}
DE_test_quantile <- function (fit, contrast = "Gender_3.B6"){
  contrasts.age_sex_str <- makeContrasts(
    Gender_3.B6    = (M.3.B6   - F.3.B6),
    Gender_3.NZO   = (M.3.NZO  - F.3.NZO),
    Gender_18.B6   = (M.18.B6  - F.18.B6),
    Gender_18.NZO  = (M.18.NZO - F.18.NZO),
    Age18vs3_M.B6  = (M.18.B6  - M.3.B6),
    Age18vs3_M.NZO = (M.18.NZO - M.3.NZO),
    Age18vs3_F.B6  = (F.18.B6  - F.3.B6),
    Age18vs3_F.NZO = (F.18.NZO - F.3.NZO),
  levels = fit$design)
  
  tmp <- contrasts.fit(fit, contrasts.age_sex_str[, contrast])
  tmp <- eBayes(tmp)
  return (tmp)
}

```

```{r DE_toptags}
#' @param adjust.method default is BH, check p.adjust doc for more
#' @param p.value determines FDR threshold
#' @param n number of returned genes
DE_toptags <- function(qlf, p.value = 1, adjust.method = "BH", n = Inf){
  top.tags <- topTags(qlf, n = n, adjust.method = adjust.method, p.value = p.value)
  return (top.tags)
}
```

```{r DE_toptags_quantile}
DE_toptags_quantile <- function (tmp, sort.by = "p", n = Inf, p.value = 1){
  top.table <- topTable(tmp, sort.by = sort.by, n = Inf, p.value = p.value)
  return(top.table)
}
```

```{r DE_annotate_genes}

DE_annotate_genes <- function (top.table, gene.names){
  
  opening_gene_locs <- top.table[top.table$logFC > 0,] %>% rownames %>% as.numeric
  closing_gene_locs <- top.table[top.table$logFC < 0,] %>% rownames %>% as.numeric
  
  opening_genes <- gene.names[opening_gene_locs] %>% as.character
  closing_genes <- gene.names[closing_gene_locs] %>% as.character
  
  opening_genes <- cbind(Gene.Name = opening_genes, top.table[top.table$logFC > 0,])
  closing_genes <- cbind(Gene.Name = closing_genes, top.table[top.table$logFC < 0,])
  
  return(list(opening_genes = opening_genes, closing_genes = closing_genes))
}
```

```{r check_genesets}
check_genesets <- function (genes, tissue_cell_type, contrast, union_size = 20e3){
  
  # Up regulated genes specs
  genes_up_tbl   <- genes[[1]] 
  # Down regulated genes specs
  genes_down_tbl <- genes[[2]]
  
  # Change between TMM and Quantile fitting
  loc_up <- match(c("P.Value","adj.P.Val"), colnames(genes_up_tbl))
  if(!is.na(loc_up[1])) {
   colnames(genes_up_tbl)[c(loc_up)] <- c("PValue","FDR")
   colnames(genes_down_tbl)[c(loc_up)] <- c("PValue","FDR")
  }
  
  # These are mice gene names which are differentially expressed.
  genes_up   <- genes_up_tbl  [, "Gene.Name"] %>% as.character
  genes_down <- genes_down_tbl[, "Gene.Name"] %>% as.character
  

  # create an empty dataframe for enriched modules,
  # so that later we can sort them and prepare excel tables, yey!
  enriched_modules_df <- data.frame()
  
  # for each geneset, iterate each module
  for (i in 1:length(selected_genesets)){
    
    geneset_name <- names(selected_genesets_labels)[[i]]
    module_names <- selected_genesets_labels[[i]]
    modules      <- selected_genesets_mice  [[i]]
    
    modules <- merge(modules, module_names, by = "Module.ID")
    
    # here we create empty p values, geneset name vectors so that we can adjust p values later
    module_count        <- nrow(module_names)
    vector_p_up         <- vector_p_down <- numeric(module_count) 
    vector_module_names <- vector_geneset_names <- character(module_count)
    # iterate the modules
    
    cat ("Geneset Name:", geneset_name, "\n")
    for (j in 1:nrow(module_names)){
      
      module_ID   <- module_names[j, "Module.ID"]
      module_name <- module_names[j, "Module.Name"]
      module      <- modules[modules$Module.Name %in% module_name, "Gene.stable.ID"]
      
      # module gene count (constant for up/down)
      gene_count_module <- unique(module) %>% length
      
      # upregulated gene count
      n_up <- length(genes_up) 
      
      # overlapped up-regulated genes with the module
      q_up <- genes_up %in% module %>% sum

      # here we calculate the probability of having a bigger intersection
      # than the count of overlapping genes given the module size and the total gene count.
      # we substract 1 for removing the equality when the lower.tail = F, which changes P(X<x) to 1-P(X>=x).
      p_up <- phyper(q_up-1, gene_count_module, union_size - gene_count_module, n_up,
                     lower.tail = F, log.p = F)
      
      # upregulated gene count
      n_down <- length(genes_down) 
      
      # overlapped down-regulated genes with the module
      q_down <- genes_down %in% module %>% sum
      
      # downregulated gene count
      p_down <- phyper(q_down-1, gene_count_module, union_size - gene_count_module, n_down,
                       lower.tail = F, log.p = F)
      
      vector_p_up[j]          <- p_up
      vector_p_down[j]        <- p_down
      vector_geneset_names[j] <- names(selected_genesets)[[i]]
      vector_module_names[j]  <- module_name
      
    }
    
    df_up   <- data.frame(geneset.name = vector_geneset_names, 
                module.name = vector_module_names,
                p = vector_p_up,
                stringsAsFactors = F)
    
    df_down <- data.frame(geneset.name = vector_geneset_names, 
                module.name = vector_module_names,
                p = vector_p_down,
                stringsAsFactors = F)
    
    # adjust the p-values for each module
    df_up$adj.p   <- p.adjust(p = df_up$p  , method = "fdr")
    df_down$adj.p <- p.adjust(p = df_down$p, method = "fdr")
    
    # sort according to adjusted p-values and then to p-values
    df_up   <- df_up  [order(df_up$adj.p  , df_up$p)  ,]
    df_down <- df_down[order(df_down$adj.p, df_down$p),]
    
    # this is an important parameter since it affects the number of modules that are chosen
    fdr.threshold <- 0.05
    
    # check if any modules are enriched for up regulated genes
    if (any(df_up$adj.p < fdr.threshold)){
      
      # take the enriched modules
      enriched_modules <- df_up[df_up$adj.p < fdr.threshold,]
      
      # add enriched modules to dataframe
      enriched_modules_df <- cbind(TCT = tissue_cell_type, 
                                   Contrast = contrast, 
                                   enriched_modules,
                                   Status="Opening",
                                   Overlapping.Genes = NA) %>% rbind(enriched_modules_df) 
      
      for (k in 1:nrow(enriched_modules)){

        enriched_module_name   <- enriched_modules[k, "module.name"]
        module.genes <- modules[modules$Module.Name %in% enriched_module_name,]
        module.genes <- merge(module.genes, genes_up_tbl, 
                              by.x = "Gene.stable.ID", by.y = "Gene.Name", all.x = T)
        module.genes <- module.genes %>% arrange(FDR)
        
        enriched_modules_df$Overlapping.Genes[k] <- 
          module.genes[!is.na(module.genes$FDR),]$GeneName %>% 
          paste(collapse = ",")
        
        filename <- paste0("output/RNAseq_Enrichment/fm_seperate/",
                           tissue_cell_type,"-",geneset_name,
                           "-", contrast,"-", 
                           gsub("/", "_", enriched_module_name),"-UpRegulated.csv")
        cat (paste0("\tSaving ", filename, "\n"))
        write.csv(module.genes, file = filename)
      }
    }
    
    # check if any modules are enriched for down regulated genes
    # To-Do: I know I should make this a function...
    if (any(df_down$adj.p < fdr.threshold)){
      
      # take the enriched modules
      enriched_modules <- df_down[df_down$adj.p < fdr.threshold,]
      
     # add enriched modules to dataframe
      enriched_modules_df <- cbind(TCT = tissue_cell_type, 
                                   Contrast = contrast, 
                                   enriched_modules,
                                   Status="Closing",
                                   Overlapping.Genes = NA) %>% rbind(enriched_modules_df)  
      
      for (k in 1:nrow(enriched_modules)){
        
        enriched_module_name   <- enriched_modules[k, "module.name"]
        module.genes <- modules[modules$Module.Name %in% enriched_module_name,]
        module.genes <- merge(module.genes, genes_down_tbl, 
                              by.x = "Gene.stable.ID", by.y = "Gene.Name", all.x = T)
        module.genes <- module.genes %>% arrange(FDR)
        
        enriched_modules_df$Overlapping.Genes[k] <- 
          module.genes[!is.na(module.genes$FDR),]$GeneName %>% 
          paste(collapse = ",")
        
        filename <- paste0("output/RNAseq_Enrichment/fm_seperate/",
                           tissue_cell_type,"-",geneset_name,
                           "-", contrast,"-", 
                           gsub("/", "_", enriched_module_name),"-DownRegulated.csv")
        cat (paste0("\tSaving ", filename, "\n"))
        write.csv(module.genes, file = filename)
      }
    }
  }
  if (nrow(enriched_modules_df) == 0) return (NULL)
  cat ("\t\t\tEnriched module count:", nrow(enriched_modules_df), "\n")
  return (enriched_modules_df)
}
```

```{r er_plot}
er_plot <- function(){
  
  plot_and_save <- function(df_, strain_name){
    
    for (var in unique(df_$geneset.name)){
      p <- ggplot(df_[ df_$geneset.name == var,],
                  aes(x=module.name, y=Contrast, color = Status, size = -log10(p))) + 
      geom_point() + 
      coord_flip() + 
      scale_color_manual(values=c("blue", "red")) + ggtitle(var) + facet_wrap(~TCT) + theme_minimal()
      
      if (var == 'wp'){
        width_ = 20
        height_ = 20
      } else {
        width_ = 10
        height_ = 6
      }
      ggsave(paste0('output/F3/GSEA/RNAseq/fm_seperate/', var,"-",strain_name, '.pdf'),
          plot = p,
          units = "in",
          width = width_,
          height = height_,
          useDingbats = FALSE
          )
    }
  }
  
  
  file_list <- list.files(path)
  all_modules <- lapply(file_list, function(x){
    read.csv(paste0(path, x), stringsAsFactors = F) %>% data.frame
  })
  
  df <- do.call("rbind", all_modules)
  df$TCT[df$TCT == "naive"] <- "CD8+ Naive"
  df$TCT[df$TCT == "memory"] <- "CD8+ Memory"
  df$TCT[df$TCT == "bm"] <- "BM"
  df$TCT[df$TCT == "spleen"] <- "SPLEEN"
  df$TCT[df$TCT == "pbl"] <- "PBL"
  
  tissue_only <- T
  if (tissue_only){
    df <- df[df$TCT == "BM" | df$TCT == "SPLEEN" | df$TCT == "PBL",]
  }
  df_b6  <- df[df$Contrast == "Age18vs3_M.B6"  | df$Contrast == "Age18vs3_F.B6",]
  df_nzo <- df[df$Contrast == "Age18vs3_M.NZO" | df$Contrast == "Age18vs3_F.NZO",]
  
  plot_and_save(df_b6, "b6")
  plot_and_save(df_nzo, "nzo")
  
}
```


```{r}
#' @param doTMM if set TRUE, do all analyses with TMM which does not require limma package.
#' @param fit_cache given a parameter it will record trained models for each tissue.
#' @param tissue_cell_type pbl, spleen, naive, memory, bm
#' @return fit model matrix

run_tissue <- function(tissue_cell_type, fit_cache = NULL, doTMM = FALSE){
  
  data <- preprocess_data (tissue_cell_type)
  
  specs        <- data$specs
  gene.names   <- data$gene.names
  count.matrix <- data$count.matrix
  
  if (is.null(fit_cache)){
    if (doTMM) fit <- DE_fit(count.matrix, specs) 
    else       fit <- DE_fit_quantile(count.matrix, specs)
  } else {
    fit <- fit_cache[[tissue_cell_type]]
  }
  
  contrasts <- c("Gender_3.B6","Gender_3.NZO","Gender_18.B6","Gender_18.NZO",
                 "Age18vs3_M.B6","Age18vs3_M.NZO","Age18vs3_F.B6","Age18vs3_F.NZO")
  
  er_modules_list <- lapply(contrasts, function(contrast, tissue_cell_type){
    
      cat (paste0("For ", contrast,":\n"))
      tmp       <- DE_test_quantile(fit, contrast)
      top.table <- DE_toptags_quantile(tmp, p.value = 0.05)
      if (!is.null(top.table)){
        genes <- DE_annotate_genes(top.table, gene.names)
        cat (paste0("\tupreg_genes ----",   genes$opening_genes %>% nrow))
        cat (paste0("\tdownreg_genes ----", genes$closing_genes %>% nrow, "\n"))
        if (genes$opening_genes %>% nrow > 0 & genes$closing_genes %>% nrow > 0){
            er_modules <- check_genesets(genes, tissue_cell_type, contrast)
            return (er_modules)
        }
      }  
  }, tissue_cell_type = tissue_cell_type)

  # make the list a table!
  er_modules <- do.call("rbind", er_modules_list)
  filename_er_modules <- paste0("output/RNAseq_Enrichment/fm_seperate/er_modules/",
                                tissue_cell_type, "_er_summary.csv")
  write.csv(er_modules, file = filename_er_modules)
  
  return(fit)
}
```




```{r run_analyses}
list <- c("naive", "memory", "pbl", "spleen")
# This is an important parameter, if you want to train all fit matrices
# from strach you need to make this FALSE, so it won't use the cached models!
use_fit_cache = TRUE
load_geneset = FALSE
# PART 1
if (use_fit_cache){
  load("analysis/cache/rna_fit_matrices.RData")
  load("analysis/cache/rna_enrichment_analysis.Rdata")

  lapply(list, function(tissue_cell_type, fit_cache){
      run_tissue(tissue_cell_type, fit_cache)
    }, fit_cache) %>% invisible

} else {
  if (load_geneset)  load_genesets() # may take a while
  else load("analysis/cache/rna_enrichment_analysis.Rdata")
  
  fit_cache <- lapply(list, function(tissue_cell_type){
    run_tissue(tissue_cell_type)
  })
  names(fit_cache) <- list
  save(fit_cache, file = "analysis/cache/rna_fit_matrices.RData")
  save(human, mouse, union_size, selected_genesets, selected_genesets_mice, selected_genesets_labels,
       file = "analysis/cache/rna_enrichment_analysis.Rdata")
}

# PART 2
# Uncomment to compare TMM and Quantile Normalization
# lapply(list, DE_  correlation) %>% invisible
  
# PART 3
# Barplots 
# Before you can use it, you may need to run the PART 1 with use_fit_cache = F
# so that a fit_cache matrix will be generated and saved.
# if (use_fit_cache){
#   load("analysis/cache/fit_matrices.RData")
#   all_bar_plots <- lapply(list, function(tissue_cell_type, fit_cache){
#     create_bar_plots(tissue_cell_type, fit_cache)
#   }, fit_cache) %>% invisible
#   names(all_bar_plots) <- list
#   draw_bar_plots(all_bar_plots)
# }

# PART 4
# Barplots for enrichment analyses

```