cont. cp enrichment tool

common/cp_enrichment.R

@@ -13,14 +13,15 @@ Options:
 '
 
 opts <- docopt(doc, commandArgs(TRUE)) ## does not work when spining
-#opts <- docopt(doc, "--overlay_expr_data ../ctrl_fc_expr_filtered.txt geneClusters.RData" )
-#opts <- docopt(doc, "--overlay_expr_data ../ctrl_fc_expr_filtered.txt grpdGoiClusters.RData" )
+# opts <- docopt(doc, "--overlay_expr_data ctrl_fc_expr_filtered.txt geneClusters.RData" )
 
-devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.12/R/core_commons.R")
-devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.12/R/ggplot_commons.R")
-devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.12/R/bio/diffex_commons.R")
+devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.13/R/core_commons.R")
+devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.13/R/ggplot_commons.R")
+devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.13/R/bio/diffex_commons.R")
 
 require.auto(knitr)
+require.auto(clusterProfiler)
+require.auto(ReactomePA)
 
 ## to fix child support issue with knitr, see also
 ## http://stackoverflow.com/questions/20030523/knitr-nested-child-documents
@@ -38,6 +39,8 @@ if(!is.null(opts$overlay_expr_data)){
     overlayData <- read.delim(opts$overlay_expr_data)
 }
 
+## TODO expose options to run gesa instead (by assuming that gene lists are sorted)
+## see http://www.bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.pdf for details
 
 resultsBaseName=basename(opts$grouped_gene_lists_rdata) %>% trim_ext(".RData") %>% paste0(".")
 
@@ -57,8 +60,68 @@ geneLists %>% inner_join(listLabels) %>%
     ylab("")
 
 
-#enrResults <- geneLists %>% do(davidAnnotationChart(.$ensembl_gene_id))
+guess_cp_species <- function(ensIds){
+    an_id <-ensIds[1]
+
+    if(str_detect(an_id, "ENSG")){
+        return("human")
+    }else if(str_detect(an_id, "ENSMUSG")){
+        return("mouse")
+    }else if(str_detect(an_id, "ENSDARG")){
+        return("zebrafish")
+    }else if(str_detect(an_id, "FBgn")){
+        return("fly")
+    }else{
+        stop(paste("could not clusterProfiler species name from ", an_id))
+    }
+}
+
+guess_anno_db <- function(ensIds){
+    an_id <-ensIds[1]
+
+    if(str_detect(an_id, "ENSG")){
+        return("org.Hs.eg.db")
+    }else if(str_detect(an_id, "ENSMUSG")){
+        return("org.Mm.eg.db")
+    }else if(str_detect(an_id, "ENSDARG")){
+        return("org.Dr.eg.db")
+    }else if(str_detect(an_id, "FBgn")){
+        return("org.Dm.eg.db")
+    }else{
+        stop(paste("could not anno db mart from ", an_id))
+    }
+}
+
+#source("https://bioconductor.org/biocLite.R")
+#biocLite("org.Mm.eg.db")
+#biocLite("org.Hm.eg.db")
+#biocLite("org.Dr.eg.db")
+#biocLite("org.Dm.eg.db")
+
+
+cpSpecies <- guess_cp_species(geneLists$ensembl_gene_id)
+annoDb <- guess_anno_db(geneLists$ensembl_gene_id) # e.g. "org.Hs.eg.db"
+
+## supported ids
+#idType("org.Hs.eg.db")
+
+## convert to entrez gene ids
+geneLists %<>% mutate(entrez_gene_id=bitr(x, fromType="ENSEMBL", toType="ENTREZID", annoDb=annoDb))
+
+
+
+## TODO  expose as argument
+pCutoff <- 0.05
+
+
+geneIds <- geneLists %>% filter(cluster %in% c("cluster_2")) %$% geneLists
+enrichKEGG(gene = geneIds, organism = cpSpecies, pvalueCutoff = pCutoff, readable = TRUE)
+enrichPathway(gene = geneIds, organism = cpSpecies, pvalueCutoff = pCutoff, readable = TRUE)
+
+enrichGO(gene = gene, universe = names(geneList), organism = "human", ont = "CC", pAdjustMethod = "BH", pvalueCutoff = 0.01, qvalueCutoff = 0.05, readable = TRUE)
+
 enrResults <- quote(geneLists %>% do(davidAnnotationChart(.$ensembl_gene_id))) %>%
+
     cache_it(paste0("enrdata_", digest(geneLists)))