pval calculation corrected

enrichment/enrich_maps.R

@@ -11,7 +11,10 @@ library(httr)
 ask_GET <- function(furl, fask) {
   resp <- httr::GET(url = paste0(furl, fask),
                     httr::add_headers('Content-Type' = "application/x-www-form-urlencoded"))
-  return(httr::content(resp, as = "text"))
+  if(httr::status_code(resp) == 200) {
+    return(httr::content(resp, as = "text"))
+  }
+  return(NULL)
 }
 
 args = commandArgs(trailingOnly=TRUE)
@@ -34,14 +37,20 @@ cat(file = "enriched_disease_maps.txt",
 
 for(map in config[config$type == "map","resource"]) {
   message(paste0("Querying ", map))
-  cfg <- fromJSON(ask_GET(map, "configuration/"))
+  cfg <- ask_GET(map, "configuration/")
+  ### For erroneous response, skip to next map
+  if(is.null(cfg)) { next }
+  cfg <- fromJSON(cfg)
   project_id <- cfg$options[cfg$options$type == "DEFAULT_MAP","value"]
   
   mnv_base <- paste0(map,"projects/",project_id,"/")
   
   ### Ask for models
   message(paste0("Asking for models: ", mnv_base, "models/"))
-  models <- fromJSON(ask_GET(mnv_base, "models/"), flatten = F)
+  models <- ask_GET(mnv_base, "models/")
+  ### For erroneous response, skip to next map
+  if(is.null(models)) { next }
+  models <- fromJSON(models, flatten = F)
   
   model_elements <- lapply(models$idObject, 
                             function(x) fromJSON(ask_GET(paste0(mnv_base,"models/",x,"/"), 
@@ -51,6 +60,8 @@ for(map in config[config$type == "map","resource"]) {
   message(paste0("Setting up enrichment parameters..."))
   ### Prep run, go through all model elements and:
   for(me in model_elements) {
+    ### For erroneous response, skip to next model
+    if(is.null(me)) { next }
     ### Get all HGNCs (universe for enrichment)
     all_hgncs <- c(all_hgncs, sapply(me$references,
                                      function(x) x[x$type == "HGNC_SYMBOL", "resource"]))
@@ -58,14 +69,15 @@ for(map in config[config$type == "map","resource"]) {
                                   function(x) ifelse(is.character(x) & length(x) > 0, TRUE, FALSE))]
   }
   
+  this_input <- intersect(input, all_hgncs)
   N <- length(all_hgncs)
-  n <- length(input)
+  n <- length(this_input)
   
   message(paste0("Going through individual models..."))
   for(i in 1:length(model_elements)) {
     me <- model_elements[[i]]
     if(length(unique(me[me$type == "Pathway","name"])) == 0) { next }
-    path_pVals <- data.frame(pname = unique(me[me$type == "Pathway","name"]), pVal = 1)
+    path_pVals <- data.frame(pname = unique(me[me$type == "Pathway","name"]), pVal = 1, stringsAsFactors = F)
     ### Calculating the inclusion of the elements per bioentity set
     for(up in 1:nrow(path_pVals)) {
       res <- apply(me[me$name == path_pVals$pname[up],"bounds"], 1, 
@@ -74,16 +86,17 @@ for(map in config[config$type == "map","resource"]) {
 
       hgncs <- unlist(sapply(me$references[res & me$type %in% c("Protein", "RNA", "Gene")], 
                              function(x) x[x$type == "HGNC_SYMBOL", "resource"]))
+      
       m <- length(hgncs)
-      k <- length(intersect(input,hgncs))
+      k <- length(intersect(this_input,hgncs))
       
       if (k > 0) {
         path_pVals$pVal[up] <- (((m * (n - k)) / (k * (N - m))) ** k) * (((n * (N - m)) / (N * (n - k))) ** n)
       }
     }
     ### Bonferroni correction
-    # path_pVals$pVal <- path_pVals$pVal*nrow(path_pVals)
-
+    path_pVals$pVal <- path_pVals$pVal*nrow(path_pVals)
+    
     path_pVals <- path_pVals[path_pVals$pVal < 0.1, ]
 
     if(nrow(path_pVals) > 0) {