ettore gsea/ora

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: Ulisse
 Type: Package
 Title: Pathway, pathway cross-talk and cell-cell communication 
-Version: 1.2.5
-Date: 2022-03-28
+Version: 1.3.0
+Date: 2023-11-27
 Authors@R: 
     c(person(given = "Alice",
            family = "Chiodi",

NAMESPACE

@@ -18,14 +18,12 @@ export(gsea)
 export(ora)
 export(ora2enrich)
 export(ora_clusters)
-export(ora_pipeline)
 export(pathway_data)
 export(pathway_sim_comm)
 export(pathway_similarity)
 export(perm_link)
 export(plot_functional_relevance)
 export(plot_network_CT)
-export(plot_ora_comparison)
 export(preparing_DEG_list)
 export(preparing_expr_list)
 export(preparing_gs_list)
@@ -41,12 +39,15 @@ import(grid)
 import(igraph)
 import(kit)
 import(msigdbr)
+import(openxlsx)
 import(pals)
 import(parallel)
 import(plotrix)
 import(stats)
 import(stringr)
 importClassesFrom(DOSE,enrichResult)
+importClassesFrom(DOSE,gseaResult)
+importFrom(ComplexHeatmap,Heatmap)
 importFrom(RColorBrewer,brewer.pal)
 importFrom(circlize,colorRamp2)
 importFrom(ggforce,geom_mark_rect)
@@ -55,18 +56,19 @@ importFrom(grDevices,adjustcolor)
 importFrom(grDevices,dev.off)
 importFrom(grDevices,jpeg)
 importFrom(grDevices,rainbow)
+importFrom(grid,gpar)
+importFrom(grid,grid.text)
 importFrom(gtools,permutations)
 importFrom(methods,as)
 importFrom(methods,is)
 importFrom(methods,new)
-importFrom(parallel,mclapply)
+importFrom(pals,brewer.purples)
 importFrom(qvalue,qvalue)
 importFrom(reshape2,acast)
 importFrom(scales,rescale)
 importFrom(stats,dhyper)
+importFrom(stats,p.adjust)
 importFrom(stats,phyper)
 importFrom(stats,setNames)
 importFrom(stringi,stri_c)
 importFrom(utils,write.table)
-importFrom(viridis,turbo)
-importFrom(viridis,viridis)

R/calc_gs_perm.R

@@ -1,11 +1,10 @@
 #' Internal function of pathtools
-#' @param rll numeric matrix of genes-by-ranking criteria; 
-#' each column contains numeric values; rownames are mandatory
+#' @param rll list of named ranked vectors
 #' @param perm vector of permuted names
 #' @param gs gene set
-calc_gs_perm <- function(rll, perm, gs){
-  out <- unlist(lapply(rll, function(x) es(which(perm %in% gs), 
-                                           array(x, dimnames=list(perm)))))
+calc_gs_perm <- function(rll=NULL, perm=NULL, gs=NULL){
+
+  out <- unlist(lapply(rll, function(x) es(idx = which(perm %in% gs), array(x, dimnames=list(perm)))[, 1]))
 
   return(out)
 

R/es.R

@@ -2,19 +2,18 @@
 #'
 #' @param idx vector of indices a subset of elements of x
 #' @param x named vector, ranked list
-#' @param le logical
-#' @return enrichment score
+#' @return data.frame with as, enrichmet score; tags, leading edge size; tags_perc, leading edge size percent over gene set; list_top, rank of the ES; list_top_perc, rank of the ES percent over full ranked list; lead_edge, gene names of the leading edge.
 #' @export
 #'
-es <- function(idx, x, le=F){
-
+es <- function(idx=NULL, x=NULL){
+  
   #idx: indexes of a subset of elements of x
   #x: array of elements
   ## ES score
   #Phit(S,i) <- SUM_{i in idx, j=1..i}( r_j / Nr)
   #r: score
   #Nr: total score in x[idx]
-
+  
   N <- length(x)
   Nh <- length(idx)
   hits <-  rep(0, length(x))
@@ -23,60 +22,52 @@ es <- function(idx, x, le=F){
   misses[idx] <- 0
   hits.cumsum <- cumsum(abs(hits))
   Nr <- sum(abs(hits)) #equal to sum(abs(x[idx]))
-
+  
   if(Nr==0){   #nothing to do
-
+    
     es <- 0
-    if(le){
-      deviation <- 0
-      tags <- 0
-      tags_perc <- 0
-      list_top <- 0
-      list_top_perc <- 0
-      lead_edge <- 0
-      lead_edge_subset <- 0
-    }
-
+
+    deviation <- 0
+    tags <- 0
+    tags_perc <- 0
+    list_top <- 0
+    list_top_perc <- 0
+    lead_edge <- 0
+    lead_edge_subset <- 0
+
   }else{
-
+    
     misses.cumsum <- cumsum(misses)
     Phits <- hits.cumsum / Nr
     Pmiss <- misses.cumsum / (N - Nh)
     deviation <- Phits - Pmiss
     wm <- which.max(abs(deviation))
-
+    
     #es <- deviation[which.max(abs(deviation))]
     es <- deviation[wm]
-
-    if(le){
-      if(es >=0){
-
-        tags <- sum(idx <= wm)
-        list_top <- wm
-        lead_edge_subset <- paste0(names(x)[intersect(1:wm, idx)], collapse = ";")
-
-      }else{
-
-        tags <- sum(idx > wm)
-        list_top <- N - wm
-        lead_edge_subset <- paste0(names(x)[intersect(wm:N, idx)], collapse = ";")
-
-      }
-
-      tags_perc <- tags / length(idx)
-      list_top_perc <- list_top / length(x)
+
+
+    if(es >=0){
+
+      tags <- sum(idx <= wm)
+      list_top <- wm
+      lead_edge_subset <- paste0(names(x)[intersect(1:wm, idx)], collapse = ";")
+
+    }else{
+
+      tags <- sum(idx >= wm)
+      list_top <- N - (wm-1)
+      lead_edge_subset <- paste0(names(x)[intersect(wm:N, idx)], collapse = ";")
+
     }
+
+    tags_perc <- tags / length(idx)
+    list_top_perc <- list_top / length(x)
+
   }
-
-  if(le){
-
-    lead_edge <- tags_perc * (1-list_top_perc) * N / (N - Nh)
-
-    return(list(es=es, deviation=deviation, lea=data.frame(tags=tags, tags_perc=tags_perc, list_top=list_top, list_top_perc=list_top_perc, lead_edge=lead_edge, lead_edge_subset=lead_edge_subset, stringsAsFactors = F)))
-
-  }else{
-
-    return(es)
-  }
-
+
+  lead_edge <- tags_perc * (1-list_top_perc) * N / (N - Nh)
+
+  return(data.frame(es=es, tags=tags, tags_perc=tags_perc, list_top=list_top, list_top_perc=list_top_perc, lead_edge=lead_edge, lead_edge_subset=lead_edge_subset, stringsAsFactors = F))
+
 }

R/filter_gsl.R

@@ -6,11 +6,14 @@
 #' @export
 
 
-filter_gsl <- function(gsl, universe, min_size=5, max_size=500){
+filter_gsl <- function(gsl=NULL, universe=NULL, min_size=5, max_size=500){
 
   ans <- lapply(gsl, function(x) x[x %in% universe])
   ans <- lapply(ans, unique)
-  idx_keep <- unlist(lapply(ans, function(x) length(x) >= min_size & length(x) <= max_size))
+  #idx_keep <- unlist(lapply(ans, function(x) length(x) >= min_size & length(x) <= max_size))
+
+  gs_size <- lengths(ans)
+  idx_keep <- gs_size >= min_size & gs_size <= max_size
 
   ans <- ans[idx_keep]
 

R/get_genes.R

@@ -0,0 +1,13 @@
+#' Internal function
+
+get_genes <- function(gene_set=NULL, wb=NULL){
+
+  ans <- gene_set[gene_set %in% wb]
+  if(!is.null(eg2sym)){
+    ans <- sort(eg2sym$symbol[eg2sym$gene_id %in% ans])
+  }
+
+  ans <- paste0(ans, collapse = ";")
+  return(ans)
+
+}