add hex sticker [skip ci]

.github/workflows/rworkflows.yml

@@ -4,14 +4,17 @@ name: rworkflows
     branches:
     - master
     - main
+    - devel
     - RELEASE_**
   pull_request:
     branches:
     - master
     - main
+    - devel
     - RELEASE_**
 jobs:
   rworkflows:
+    permissions: write-all
     runs-on: ${{ matrix.config.os }}
     name: ${{ matrix.config.os }} (${{ matrix.config.r }})
     container: ${{ matrix.config.cont }}
@@ -22,8 +25,8 @@ jobs:
         - os: ubuntu-latest
           bioc: devel
           r: auto
-          cont: bioconductor/bioconductor_docker:devel
-          rspm: https://packagemanager.rstudio.com/cran/__linux__/focal/release
+          cont: ghcr.io/bioconductor/bioconductor_docker:devel
+          rspm: ~
         - os: macOS-latest
           bioc: release
           r: auto
@@ -46,10 +49,9 @@ jobs:
         run_pkgdown: ${{ true }}
         has_runit: ${{ false }}
         has_latex: ${{ false }}
-        GITHUB_TOKEN: ${{ secrets.PAT_GITHUB }}
+        GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
         run_docker: ${{ true }}
-        docker_user: bschilder
-        docker_org: neurogenomicslab
         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.6
+Version: 0.1.7
 Authors@R:
     c(person(given = "Robert",
            family = "Gordon-Smith",
@@ -64,7 +64,11 @@ Suggests:
     grDevices,
     rvest,
     patchwork,
-    ggnetwork
+    ggnetwork,
+    ggpubr,
+    orthogene,
+    readxl,
+    tidyr
 Remotes: 
   github::neurogenomics/HPOExplorer,
   github::NathanSkene/EWCE

NAMESPACE

@@ -17,13 +17,15 @@ export(load_example_results)
 export(map_tissues)
 export(merge_results)
 export(ontology_plot)
+export(plot_ont_lvl)
 export(prioritise_targets)
 export(prioritise_targets_network)
 export(report_plot)
 export(subset_phenos)
 export(subset_results)
 export(summary_plot)
 export(terminal_celltypes)
+export(ttd_check)
 import(GeneOverlap)
 import(HPOExplorer)
 import(data.table)
@@ -69,3 +71,4 @@ importFrom(tools,R_user_dir)
 importFrom(utils,capture.output)
 importFrom(utils,getFromNamespace)
 importFrom(utils,head)
+importFrom(utils,tail)

NEWS.md

@@ -1,3 +1,14 @@
+# MultiEWCE 0.1.7
+
+## New features
+
+* Update to coordinate with `HPOExplorer` updates.
+* New funcs:
+  - `plot_ont_lvl`
+  - `ttd_check`
+  - `ttd_plot`
+  - `ttd_import`
+
 # MultiEWCE 0.1.6
 
 ## New features
@@ -11,6 +22,7 @@
   - "Modifier" --> "modifier"
   - "Aspect" --> "aspect"
   - "Gene" --> "gene_symbol"
+  - "DatabaseID" --> "disease_id" 
   - "LinkID" --> "disease_id"
 * Update all "data" objects.
 * `get_data`
@@ -23,6 +35,8 @@
   - Update colnames dynamically.
 * `load_example_ctd`
   - Change `tag` to "latest".
+* `load_hpo_graph`
+  - Change `tag` to "latest".
 * `map_tissues`
   - Fix docs.
 * `agg_results`

R/ewce_para.R

@@ -28,7 +28,7 @@
 #' @examples
 #' gene_data <- HPOExplorer::load_phenotype_to_genes()
 #' ctd <- MultiEWCE::load_example_ctd()
-#' list_names <- unique(gene_data$hpo_name)[seq_len(3)]
+#' list_names <- unique(gene_data$hpo_name)[seq(3)]
 #' res_files <- ewce_para(ctd = ctd,
 #'                        gene_data = gene_data,
 #'                        list_names = list_names,

R/gen_overlap.R

@@ -21,7 +21,7 @@
 #' @importFrom parallel mclapply
 #' @examples
 #' gene_data <- HPOExplorer::load_phenotype_to_genes()
-#' list_names <- unique(gene_data$disease_id)[seq_len(3)]
+#' list_names <- unique(gene_data$disease_id)[seq(3)]
 #' overlap <- gen_overlap(gene_data = gene_data,
 #'                        list_names = list_names)
 gen_overlap <- function(gene_data =

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_len(5)]
+#' list_names <- unique(gene_data$hpo_name)[seq(5)]
 #' ctd <- load_example_ctd()
 #' all_results <- gen_results(ctd = ctd,
 #'                            gene_data = gene_data,
@@ -58,6 +58,7 @@ gen_results <- function(ctd,
 
   # devoptera::args2vars(gen_results)
 
+  start <- Sys.time()
   #### Run analysis ####
   res_files <- ewce_para(ctd = ctd,
                          list_names = list_names,
@@ -77,10 +78,13 @@ gen_results <- function(ctd,
   #### Merge results into one dataframe ####
   results_final <- merge_results(res_files = res_files,
                                  list_name_column = list_name_column)
+  #### Report total time ####
+  messager("Done in:",round(difftime(Sys.time(),start,units = "s"), 1),
+           "seconds.",v=verbose)
   #### Save merged results ####
   save_path <- save_results(results = results_final,
-                           save_dir = save_dir,
-                           prefix = "gen_results_",
-                           verbose = verbose)
+                            save_dir = save_dir,
+                            prefix = "gen_results",
+                            verbose = verbose)
   return(results_final)
 }

R/get_unfinished_list_names.R

@@ -12,7 +12,7 @@
 #' @export
 #' @examples
 #' gene_data <- HPOExplorer::load_phenotype_to_genes()
-#' list_names <- unique(gene_data$hpo_name)[seq_len(3)]
+#' list_names <- unique(gene_data$hpo_name)[seq(3)]
 #' save_dir_tmp <- file.path(tempdir(),"results")
 #' ctd <- load_example_ctd()
 #' res_files <- ewce_para(ctd = ctd,

R/load_example_CTD.R

@@ -14,7 +14,8 @@
 #' CTD <- load_example_ctd()
 load_example_ctd <- function(file="CTD_Descartes_example.rds",
                              tag = "latest",
-                             save_dir=tools::R_user_dir(package = "MultiEWCE")
+                             save_dir=tools::R_user_dir(package = "MultiEWCE",
+                                                        which = "cache")
                              ) {
 
   dir.create(save_dir, showWarnings = FALSE, recursive = TRUE)

R/load_example_results.R

@@ -67,13 +67,15 @@
 #' @examples
 #' res <- load_example_results()
 load_example_results <- function(file=c(
+  # "results_DescartesHuman.csv.gz",
   "Descartes_All_Results_extras.symptoms.rds",
   "Descartes_All_Results_extras.symptoms.full_join.rds",
   "Descartes_All_Results_extras.rds",
   "gen_overlap.symptoms.filt.rds",
   "tabulamuris_merged.rds"),
   tag = "latest",
-  save_dir=tools::R_user_dir(package = "MultiEWCE"),
+  save_dir=tools::R_user_dir(package = "MultiEWCE",
+                             which = "cache"),
   force_new=FALSE
   ) {
 
@@ -92,7 +94,11 @@ load_example_results <- function(file=c(
                            dest = save_dir,
                            overwrite = TRUE)
   }
-  results <- readRDS(save_path)
+  if(grepl("\\.rds$",save_path, ignore.case = TRUE)){
+    results <- readRDS(save_path)
+  } else {
+    results <- data.table::fread(save_path)
+  }
   data.table::setnames(results,
                        c("HPO_ID","Phenotype","HPO_ID.disease_id","disease_id",
                          "HPO_ID.LinkID","LinkID"),

R/load_hpo_graph.R

@@ -21,8 +21,9 @@
 #' @examples
 #' g <- load_hpo_graph()
 load_hpo_graph <- function(file="hpo_igraph.rds",
-                           tag = "v0.0.1",
-                           save_dir=tools::R_user_dir(package = "MultiEWCE")
+                           tag = "latest",
+                           save_dir=tools::R_user_dir(package = "MultiEWCE",
+                                                      which = "cache")
                            ) {
 
   file <- file[[1]]

R/merge_results.R

@@ -15,7 +15,7 @@
 #' @examples
 #' gene_data <- HPOExplorer::load_phenotype_to_genes()
 #' ctd <- MultiEWCE::load_example_ctd()
-#' list_names <- unique(gene_data$hpo_name)[seq_len(3)]
+#' list_names <- unique(gene_data$hpo_name)[seq(3)]
 #' res_files <- ewce_para(ctd = ctd,
 #'                        gene_data = gene_data,
 #'                        list_names = list_names,
@@ -34,7 +34,7 @@ merge_results <- function(save_dir=NULL,
     names(res_files) <- gsub("_"," ",tolower(gsub(".rds$","",res_files)))
     messager(formatC(length(res_files),big.mark = ","),"results files found.")
   }
-  lapply(seq_len(length(res_files)),
+  lapply(seq(length(res_files)),
          function(i){
    if(is.null(res_files[[i]])){
      return(NULL)

R/plot_ont_lvl.R

@@ -0,0 +1,81 @@
+#' Plot ontology levels
+#'
+#' Generate plots comparing the ontology level of each HPO phenotype and
+#' several other metrics.
+#' @param p2g Phenotype to gene data.
+#' @param x_vars Variables to plot on the x-axis of each subplot.
+#' @returns A named list containing the data and the plot.
+#' @inheritParams prioritise_targets
+#'
+#' @export
+#' @import HPOExplorer
+#' @examples
+#' plts <- plot_ont_lvl()
+plot_ont_lvl <- function(results = load_example_results(),
+                         p2g = HPOExplorer::load_phenotype_to_genes(),
+                         x_vars = c("genes",
+                                    "celltypes",
+                                    "log(abs(fold_change))")){
+
+  requireNamespace("ggplot2")
+  requireNamespace("patchwork")
+  requireNamespace("ggpubr")
+
+  gene_symbol <- celltypes <- CellType <- ontLvl <- NULL;
+
+  results[,celltypes:=length(unique(CellType[q<0.05])),by="hpo_id"]
+  pcount <- p2g[,list(genes=length(unique(gene_symbol))),
+                by="hpo_id"]
+  pcount <- HPOExplorer::add_ont_lvl(pcount)
+  r2 <- merge(results[,c("hpo_id","CellType","celltypes",
+                         "p","q","fold_change")] |>
+                unique(),
+              pcount,
+              by="hpo_id")
+
+  plt <- function(x_var="log(genes)",
+                  y_var="log(abs(fold_change))",
+                  geom="hex",
+                  method = "loess",
+                  direction = 1,
+                  ...){
+    r2[,mean:=mean(get(gsub("[(]|[)]|log|abs","",y_var))),by="hpo_id"]
+    gp <- ggplot(r2,aes_string(x=x_var,y=y_var)) +
+      scale_fill_viridis_c(option = "plasma",
+                           direction = direction)
+    if(geom=="hex"){
+      gp <- gp + geom_hex(...)
+    } else if(geom=="violin"){
+      gp <- gp + geom_violin(orientation = "y",
+                             aes(fill=ontLvl),
+                             ...)
+    } else if(geom=="boxplot"){
+      gp <- gp + geom_boxplot(orientation = "y",
+                              aes(group=ontLvl,
+                                  fill=mean),
+                              ...)
+    }else {
+      gp <- gp + geom_jitter(width = 0,
+                             alpha=.25,
+                            ...)
+    }
+    gp <- gp +
+      geom_smooth(method = method) +
+      ggpubr::stat_cor(method = "pearson",
+                       label.x.npc = .5,
+                       label.y.npc = 1) +
+      theme_bw()
+    return(gp)
+  }
+  # plts1 <- lapply(c("log(genes)","log(celltypes)"), plt) |>
+  #   patchwork::wrap_plots(ncol = 1)
+  plts2 <- lapply(x_vars, plt,
+                  y="ontLvl",
+                  geom="boxplot",
+                  direction = -1,
+                  notch=TRUE) |>
+    patchwork::wrap_plots(ncol = 1)
+  return(list(data=r2,
+              plot=plts2))
+}
+

R/prioritise_targets.R

@@ -156,13 +156,13 @@ prioritise_targets <- function(#### Input data ####
                                  ),
                                #### Celltype level ####
                                q_threshold = 0.05,
-                               fold_threshold = 1,
+                               fold_threshold = 2,
                                symptom_p_threshold = NULL,
                                symptom_intersection_size_threshold = 1,
                                keep_celltypes = terminal_celltypes()$CellType,
                                #### Gene level ####
                                keep_evidence = seq(3,6),
-                               keep_seqnames = c(seq_len(22),"X","Y"),
+                               keep_seqnames = c(seq(22),"X","Y"),
                                gene_size = list("min"=0,
                                                 "max"=Inf),
                                gene_frequency_threshold = NULL,
@@ -181,7 +181,9 @@ prioritise_targets <- function(#### Input data ####
                                              "gene_freq_mean"=-1,
                                              "width"=1),
                                top_n = NULL,
-                               group_vars = c("hpo_id","CellType"),
+                               group_vars = c("disease_id",
+                                              "hpo_id",
+                                              "CellType"),
                                return_report = TRUE,
                                verbose = TRUE){
 

R/report_plot.R

@@ -71,7 +71,7 @@ report_plot <- function(rep_dt,
   dt1[,Tier:=gsub("Tier NA",NA,Tier)]
   dt1$step <- factor(dt1$step,
                     levels = unique(dt1$step),
-                    labels = paste0(seq_len(length(unique(dt1$step))),". ",
+                    labels = paste0(seq(length(unique(dt1$step))),". ",
                                     unique(dt1$step)),
                     ordered = TRUE)
   #### Make plot: tiers ####
@@ -101,7 +101,7 @@ report_plot <- function(rep_dt,
     data.table::melt(id.vars=names(filters))
   dt2$step <- factor(dt2$step,
                     levels = unique(dt2$step),
-                    labels = paste0(seq_len(length(unique(dt2$step))),". ",
+                    labels = paste0(seq(length(unique(dt2$step))),". ",
                                     unique(dt2$step)),
                     ordered = TRUE)
   messager("report_plot:: Preparing plot.",v=verbose)

R/save_results.R

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