allow caching in gen_results. speed up examples

.github/workflows/rworkflows.yml

@@ -54,4 +54,3 @@ jobs:
         DOCKER_TOKEN: ${{ secrets.DOCKER_TOKEN }}
         runner_os: ${{ runner.os }}
         cache_version: cache-v1
-        docker_registry: ghcr.io

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: MultiEWCE
 Title: EWCE for Multiple Gene Lists
-Version: 0.1.7
+Version: 0.1.8
 Authors@R:
     c(person(given = "Robert",
            family = "Gordon-Smith",

NAMESPACE

@@ -14,6 +14,7 @@ export(get_unfinished_list_names)
 export(ggnetwork_plot_full)
 export(load_example_ctd)
 export(load_example_results)
+export(load_hpo_graph)
 export(map_tissues)
 export(merge_results)
 export(ontology_plot)

NEWS.md

@@ -1,3 +1,24 @@
+# MultiEWCE 0.1.8
+
+## New features
+
+* `MultiEWCE` finally gets a hex sticker!
+* `gen_results` / `gen_overlap`
+  - Check for existing results and import if already there.
+  - Name all results "gen_results.rds" or "gen_overlap.rds" to avoid 
+    rerunning duplicate analyses on HPC.
+* Update *rworkflows.yml*
+
+## Bug fixes
+
+* Fix unit tests and examples to use "hpo_id" instead of "hpo_name".
+* `load_hpo_graph`: 
+  - export
+  - Regenerate and update "hpo_graph.rds" file.
+* Drastically reduce time to run examples.
+* `ontology_plot`
+  - Fix function and add test.
+
 # MultiEWCE 0.1.7
 
 ## New features

R/correlation_heatmap.R

@@ -24,7 +24,7 @@
 #' @importFrom data.table := dcast.data.table setnafill
 #' @importFrom stats cor
 #' @examples
-#' top_targets <- example_targets$top_targets
+#' top_targets <- example_targets$top_targets[seq(100),]
 #' hm <- correlation_heatmap(top_targets = top_targets)
 correlation_heatmap <- function(top_targets,
                                 row_side_vars = c("ancestor_name",

R/frequency_barplot.R

@@ -8,7 +8,8 @@
 #' @export
 #' @importFrom HPOExplorer load_phenotype_to_genes
 #' @examples
-#' fp_res <- frequency_barplot()
+#' results <- load_example_results()[seq(5000),]
+#' fp_res <- frequency_barplot(results=results)
 frequency_barplot <- function(results = load_example_results(),
                               phenotype_to_genes = load_phenotype_to_genes(),
                               show_plot = TRUE,

R/frequency_histogram.R

@@ -10,7 +10,8 @@
 #' @importFrom HPOExplorer load_phenotype_to_genes
 #' @importFrom HPOExplorer add_gene_frequency add_ancestor
 #' @examples
-#' fp_res <- frequency_histogram()
+#' results <- load_example_results()[seq(5000),]
+#' fp_res <- frequency_histogram(results=results)
 frequency_histogram <- function(results = load_example_results(),
                                 phenotype_to_genes = load_phenotype_to_genes(),
                                 show_plot = TRUE,
@@ -35,7 +36,7 @@ frequency_histogram <- function(results = load_example_results(),
     measure.vars = measure.vars )
 
   g1 <- ggplot2::ggplot(d1, ggplot2::aes(x=value, fill=variable)) +
-    ggplot2::geom_histogram(stat = "count") +
+    ggplot2::geom_histogram(stat = "count", na.rm = TRUE) +
     ggplot2::scale_fill_manual(values = pals::viridis(4)) +
     ggplot2::facet_wrap(facets = "variable ~.") +
     ggplot2::theme_bw() +

R/gen_overlap.R

@@ -35,13 +35,26 @@ gen_overlap <- function(gene_data =
                         top_n = NULL,
                         long_format = FALSE,
                         save_dir = tempdir(),
+                        force_new = FALSE,
                         cores = 1,
                         verbose = TRUE){
   # o <- devoptera::args2vars(gen_overlap); list_name_column="hpo_id.disease_id";
   # gene_data[,hpo_id.disease_id:=paste(hpo_id,disease_id,sep=".")]
 
-  qval <- pval <-  gene_symbol <- NULL;
+  qval <- pval <- NULL;
 
+  #### Create save path ####
+  save_path <- gen_results_save_path(save_dir = save_dir,
+                                     prefix = "gen_overlap")
+  #### Check if results already exist ####
+  if(file.exists(save_path) &&
+     isFALSE(force_new)) {
+    messager("Results already exist at:",save_path,
+             "Use `force_new=TRUE` to overwrite.",v=verbose)
+    results_final <- readRDS(save_path)
+    return(results_final)
+  }
+  #### Run new analysis ####
   t1 <- Sys.time()
   ct_genes <- apply(ctd[[annotLevel]]$specificity_quantiles,
                     2,
@@ -64,10 +77,12 @@ gen_overlap <- function(gene_data =
   #                                        bg = bg),
   #                      by=list_name_column]
 
-  split.data.table <- utils::getFromNamespace("split.data.table","data.table")
   #### Remove all unnecessary columns to save memory ####
   gene_data <- gene_data[,c(list_name_column,gene_column), with=FALSE]
+  #### Subset data to only the list_names ####
+  gene_data <- gene_data[get(list_name_column) %in% list_names,]
   messager("Splitting data.",v=verbose)
+  split.data.table <- utils::getFromNamespace("split.data.table","data.table")
   gene_data_split <- split.data.table(x = gene_data,
                                       by = list_name_column,
                                       keep.by = FALSE)
@@ -97,8 +112,7 @@ gen_overlap <- function(gene_data =
   messager(difftime(Sys.time(),t1),v = TRUE)
   #### Save results ####
   save_path <- save_results(results = overlap,
-                            save_dir = save_dir,
-                            prefix = "gen_overlap_",
+                            save_path = save_path,
                             verbose = verbose)
   return(overlap)
 }

R/gen_results.R

@@ -33,7 +33,7 @@
 #' @importFrom stringr str_replace_all
 #' @examples
 #' gene_data <- HPOExplorer::load_phenotype_to_genes()
-#' list_names <- unique(gene_data$hpo_name)[seq(5)]
+#' list_names <- unique(gene_data$hpo_id)[seq(5)]
 #' ctd <- load_example_ctd()
 #' all_results <- gen_results(ctd = ctd,
 #'                            gene_data = gene_data,
@@ -45,7 +45,6 @@ gen_results <- function(ctd,
                         gene_column = "gene_symbol",
                         list_names = unique(gene_data[[list_name_column]]),
                         bg = unique(gene_data[[gene_column]]),
-                        force_new = FALSE,
                         reps = 100,
                         annotLevel = 1,
                         genelistSpecies = "human",
@@ -54,10 +53,22 @@ gen_results <- function(ctd,
                         parallel_boot = FALSE,
                         save_dir_tmp = NULL,
                         save_dir = tempdir(),
+                        force_new = FALSE,
                         verbose = 1) {
 
   # devoptera::args2vars(gen_results)
 
+  #### Create save path ####
+  save_path <- gen_results_save_path(save_dir = save_dir,
+                                     prefix = "gen_results")
+  #### Check if results already exist ####
+  if(file.exists(save_path) &&
+     isFALSE(force_new)) {
+    messager("Results already exist at:",save_path,
+             "Use `force_new=TRUE` to overwrite.",v=verbose)
+    results_final <- readRDS(save_path)
+    return(results_final)
+  }
   start <- Sys.time()
   #### Run analysis ####
   res_files <- ewce_para(ctd = ctd,
@@ -83,8 +94,7 @@ gen_results <- function(ctd,
            "seconds.",v=verbose)
   #### Save merged results ####
   save_path <- save_results(results = results_final,
-                            save_dir = save_dir,
-                            prefix = "gen_results",
+                            save_path = save_path,
                             verbose = verbose)
   return(results_final)
 }

R/gen_results_save_path.R

@@ -0,0 +1,10 @@
+gen_results_save_path <- function(save_dir,
+                                  prefix="gen_results",
+                                  suffix=NULL){
+  save_path <- file.path(
+    save_dir,
+    gsub(" ","_",paste0(prefix,suffix,".rds"))
+  )
+  dir.create(save_dir, showWarnings = FALSE, recursive = TRUE)
+  return(save_path)
+}

R/load_hpo_graph.R

@@ -11,11 +11,11 @@
 #' f <- file.path(tempdir(),"hpo_igraph.rds")
 #' saveRDS(g,file = f)
 #' piggyback::pb_upload(file = f,
-#'                      tag = "v0.0.1", repo = "neurogenomics/MultiEWCE")
+#'                      tag = "latest", repo = "neurogenomics/MultiEWCE")
 #' }
 #' @returns graph object
 #'
-#' @keywords internal
+#' @export
 #' @importFrom piggyback pb_download
 #' @importFrom tools R_user_dir
 #' @examples
@@ -37,5 +37,6 @@ load_hpo_graph <- function(file="hpo_igraph.rds",
                            overwrite = TRUE)
   }
   g <- readRDS(save_path)
+  g <- igraph::upgrade_graph(g)
   return(g)
 }

R/ontology_plot.R

@@ -26,7 +26,8 @@
 #' @importFrom ontologyPlot onto_plot
 #' @importFrom HPOExplorer get_hpo load_phenotype_to_genes
 #' @examples
-#' plt <- ontology_plot(cell_type="Amacrine cells")
+#' results = load_example_results()[seq(100000)]
+#' plt <- ontology_plot(cell_type="ENS_glia", results=results)
 ontology_plot <- function(cell_type,
                           results = load_example_results(),
                           color_var = c("fold_change","q","p"),
@@ -37,21 +38,21 @@ ontology_plot <- function(cell_type,
                             HPOExplorer::load_phenotype_to_genes(),
                           palette = "Spectral",
                           shape = "rect",
+                          verbose=TRUE,
                           ...
                           ){
-  # templateR:::source_all()
   # devoptera::args2vars(ontology_plot)
 
-  HPO_term_valid <- NULL;
-
-  message("Generating ontology plot.")
+  color <- hpo_id <- NULL;
+  messager("Generating ontology plot.",v=verbose)
+  #### Prepare data ####
   cells <- subset_results(cell_type = cell_type,
-                             results = results,
-                             q_threshold = q_threshold,
-                             fold_threshold = fold_threshold,
-                             phenotype_to_genes = phenotype_to_genes,
-                             hpo = hpo)
-  ### Check color_var ####
+                          results = results,
+                          q_threshold = q_threshold,
+                          fold_threshold = fold_threshold,
+                          phenotype_to_genes = phenotype_to_genes,
+                          hpo = hpo)
+  #### Prepare colors  ####
   color_var <- color_var[[1]]
   val_opts <- eval(formals(ontology_plot)$color_var)
   if(!color_var %in% val_opts){
@@ -63,10 +64,15 @@ ontology_plot <- function(cell_type,
                                       cells[[color_var]],
                                       reverse = color_var == "fold change")
   }
+  cols <- c("hpo_id",color_var)
+  color_dat <- unique(
+    cells[hpo_id %in% unique(cells$hpo_id),][,cols, with=FALSE]
+  )
+  color_dat[,color:=color_func(get(color_var))]
   #### Create plot ####
   plt <- ontologyPlot::onto_plot(ontology = hpo,
-                                 terms = cells$hpo_id,
-                                 fillcolor = color_func(cells[[color_var]]),
+                                 terms = color_dat$hpo_id,
+                                 fillcolor = color_dat$color,
                                  shape = shape,
                                  ...)
   return(plt)

R/prioritise_targets.R

@@ -122,7 +122,8 @@
 #' @import data.table
 #' @importFrom utils head
 #' @examples
-#' res <- prioritise_targets()
+#' results = load_example_results()[seq(5000),]
+#' res <- prioritise_targets(results=results)
 prioritise_targets <- function(#### Input data ####
                                results = load_example_results(),
                                ctd = load_example_ctd(),

R/save_results.R

@@ -1,16 +1,9 @@
 save_results <- function(results,
-                         save_dir,
-                         prefix="results_",
-                         suffix=NULL,
-                         # suffix=stringr::str_replace_all(Sys.time(),":","-"),
+                         save_path,
                          verbose=TRUE){
-  if (!is.null(save_dir) &&
+
+  if (!is.null(save_path) &&
       nrow(results)>0) {
-    save_path <- file.path(
-      save_dir,
-      gsub(" ","_",paste0(prefix,suffix,".rds"))
-    )
-    dir.create(save_dir, showWarnings = FALSE, recursive = TRUE)
     messager("\nSaving results ==>",save_path,v=verbose)
     saveRDS(results,save_path)
     return(save_path)