added ensembl wrapper utiltiy for cluster profiler

R/bio/cp_utils.R

@@ -40,12 +40,23 @@ guess_anno_db <- function(ensIds){
 #biocLite("org.Dr.eg.db")
 #biocLite("org.Dm.eg.db")
 #biocLite("KEGG.db")
+#biocLite("ReactomePA")
 
 
 
-load_pack(ReactomePA)
+#load_pack(ReactomePA)
 
-cp_test <- function(geneIds){
+## clusterProfiler convenience wrapper for ensembl ids
+find_enr_terms = function(ensemblIds, ...){
+    #browser()
+    annoDb = guess_anno_db(ensemblIds)
+    clusterProfiler::bitr(ensemblIds, fromType="ENSEMBL", toType="ENTREZID", OrgDb=annoDb) %>%
+        with(cp_test(ENTREZID, annoDb=annoDb, cp_species=guess_cp_species(ensemblIds), ...))
+}
+
+
+
+cp_test = function(geneIds, annoDb, cp_species, q_cutoff=0.05){
     # DEBUG geneIds <- glMapped %>% filter(cluster %in% c("cluster_9")) %$% entrez_gene_id %>% as.integer
     # DEBUG geneIds <- head(glMapped,30)$entrez_gene_id
     #    geneIds=.$entrez_gene_id
@@ -57,14 +68,13 @@ cp_test <- function(geneIds){
     echo("testing", length(geneIds), " genes for enrichment")
 
     #    PANTHER10_ontology <- read.delim("http://data.pantherdb.org/PANTHER10.0/ontology/Protein_Class_7.0")
+    #    pantherResults <-     enricher(gene = geneIds, organism = cp_species, qvalueCutoff = q_cutoff, readable = TRUE, TERM2GENE = PANTHER10_ontology) %>% summary()
 
-    #    browser()
-    #    pantherResults <-     enricher(gene = geneIds, organism = cpSpecies, qvalueCutoff = qCutoff, readable = TRUE, TERM2GENE = PANTHER10_ontology) %>% summary()
-    keggResults <-        clusterProfiler::enrichKEGG(gene = geneIds, organism = cpSpecies, qvalueCutoff = qCutoff, use_internal_data=T) %>% summary()
-    reactomeResults <-      enrichPathway(gene = geneIds, organism = cpSpecies, qvalueCutoff = qCutoff) %>% summary()
-    goResultsCC <-          clusterProfiler::enrichGO(gene = geneIds, OrgDb = annoDb, qvalueCutoff = qCutoff, ont = "CC") %>% summary()
-    goResultsMF <-          clusterProfiler::enrichGO(gene = geneIds, OrgDb = annoDb, qvalueCutoff = qCutoff, ont = "MF") %>% summary()
-    goResultsBP <-          clusterProfiler::enrichGO(gene = geneIds, OrgDb = annoDb, qvalueCutoff = qCutoff, ont = "BP") %>% summary()
+    keggResults <-          clusterProfiler::enrichKEGG(gene = geneIds, organism = cp_species, qvalueCutoff = q_cutoff, use_internal_data=T) %>% as.data.frame()
+    reactomeResults <-      ReactomePA::enrichPathway(gene = geneIds, organism = cp_species, qvalueCutoff = q_cutoff) %>% as.data.frame()
+    goResultsCC <-          clusterProfiler::enrichGO(gene = geneIds, OrgDb = annoDb, qvalueCutoff = q_cutoff, ont = "CC") %>% as.data.frame()
+    goResultsMF <-          clusterProfiler::enrichGO(gene = geneIds, OrgDb = annoDb, qvalueCutoff = q_cutoff, ont = "MF") %>% as.data.frame()
+    goResultsBP <-          clusterProfiler::enrichGO(gene = geneIds, OrgDb = annoDb, qvalueCutoff = q_cutoff, ont = "BP") %>% as.data.frame()
 
     #cp-bug: if no pathways are enriched odd strucuture is retured ##todo file issue
     if(!("data.frame" %in% class(keggResults))) keggResults <- filter(goResultsBP, Description="foobar")
@@ -81,16 +91,24 @@ cp_test <- function(geneIds){
 }
 
 
-## example
+## example usage
 if(F){
 
-someGenes = c("") ## todo continue
+someGenes = c("ENSMUSG00000002014", "ENSMUSG00000002015", "ENSMUSG00000008892",
+"ENSMUSG00000015733", "ENSMUSG00000020720", "ENSMUSG00000020869",
+"ENSMUSG00000022884", "ENSMUSG00000026202", "ENSMUSG00000026276",
+"ENSMUSG00000028648", "ENSMUSG00000030432", "ENSMUSG00000030894",
+"ENSMUSG00000032249", "ENSMUSG00000038286", "ENSMUSG00000038965",
+"ENSMUSG00000038991")
+
+## using ensembl conveience wrapper
+enrResults = find_enr_terms(someGenes)
 
+## or the old-style way:
 cpSpecies <- guess_cp_species(someGenes)
 annoDb <- guess_anno_db(someGenes) # e.g. "org.Hs.eg.db"
-
-
 enrResults <- cp_test(someGenes)
+
 ## or for grouped data
 #enrResults <-  quote(glMapped %>% do(cp_test(.$entrez_gene_id))) %>% cache_it(paste0("enrdata_", digest(glMapped)))
 
@@ -98,4 +116,6 @@ enrResults <- cp_test(someGenes)
 #enrichrFile_TFchip = "/projects/bioinfo/holger/data/enrichr_datasets/ENCODE_TF_ChIP-seq_2015.txt"
 #tfChipOnt <- convert_enrichr_cp(enrichrFile_TFchip)
 
+
+
 }
\ No newline at end of file