fixed diffbind analysis by timepoint

misc/enr_simple.R

+
+
+########################################################################################################################
+#' ## Term enrichment
+
+#+ echo=FALSE
+
+#' This analysis was performed using [David](http://david.abcc.ncifcrf.gov/). The following ontologies were tested: `r paste(ontologies, collapse=', ')`
+
+geneLists <- degs %>%
+    transmute(ensembl_gene_id, list_id=paste(sample_1, "vs", sample_2))
+
+if(nrow(contrasts)<4){
+    geneLists <- rbind_list(
+        geneLists,
+        degs %>% filter(s1_overex) %>% transmute(ensembl_gene_id, list_id=paste(sample_1, ">", sample_2)),
+        degs %>% filter(!s1_overex) %>% transmute(ensembl_gene_id, list_id=paste(sample_1, "<", sample_2))
+    )
+}
+
+## additional overlaps before doing the intersection analysis
+geneLists %>% count(list_id) %>% kable()
+
+intersectLists <- function(geneLists, listIdA, listIdB, intersectListId) {
+    commonGenes <- setdiff(filter(geneLists, list_id==listIdA)$ensembl_gene_id, filter(geneLists, list_id==listIdB)$ensembl_gene_id)
+    data.frame(list_id=intersectListId, ensembl_gene_id=commonGenes)
+}
+
+## make this project specific setting
+#with(geneLists, as.data.frame(table(list_id))) %$% ac(unique(list_id))
+#if(any(str_detect(geneLists$list_id, "cyst"))){
+#    geneLists <- intersectLists(geneLists, "liver != unpolarised", "cyst != unpolarised", "unp_liv & unp_cys") %>% rbind_list(geneLists)
+#}
+
+geneLists <- distinct(geneLists) ## just precaution in case of multiple evals of statements above
+
+
+#+ eval=T
+enrResults <- quote(geneLists %>% group_by(list_id) %>% do(davidAnnotationChart(.$ensembl_gene_id))) %>%
+    cache_it(paste0("enrdata_", digest(geneLists)))
+
+write.delim(enrResults, file="enrResults.txt")
+# enrResults <- read.delim("enrResults.txt")
+#'  [Enrichment Results](enrResults.txt)
+
+#ac(annoChart$Genes[1]) %>% str_split(", ") %>% unlist() %>% length()
+
+#
+#' Because David is not too strigent by default we extract just those terms for which Bonferroni<0.01
+sigEnrResults <- filter(enrResults, Bonferroni <0.01)
+
+
+
+#' results='asis', error=TRUE
+if(nrow(sigEnrResults)>0){
+sigEnrResults %>% select(-Genes) %>% kable("html", table.attr = "class='dtable'", escape=F)
+}else{
+echo("no highly significant terms found")
+sigEnrResults = data.frame(Category=factor(c(), c()))
+}
+
+write.delim(sigEnrResults, file="sigEnrResults.txt")
+# sigEnrResults <- read.delim("sigEnrResults.txt")
+#'  [Very Significant Terms](sigEnrResults.txt)
+
+
+## plot the enrichment results
+#sigEnrResults %>% group_by(Category, add=T) %>% do({
+#    logPlot <- . %>% ggplot(aes(Term, PValue)) +
+#	    geom_bar(stat="identity")+coord_flip() +
+#	    xlab("Enriched Terms") +
+#	    ggtitle(.$Category[1]) +
+#	    scale_y_log10()
+#	    print(logPlot)
+#})
+
+#+ include=FALSE, eval=exists("sigEnrResults")
+## include=FALSE, eval=nrow(sigEnrResults)>1
+warning("dropping levels")
+sigEnrResults %<>% mutate(Category=ac(Category)) ## drop unsused level to get consistent color palette
+
+term_category_colors <- create_palette(unique(ac(sigEnrResults$Category)))
+
+term_barplot_files <- sigEnrResults %>%
+## chop and pad category names
+mutate(Term=str_sub(Term, 1, 70) %>% str_pad(70)) %>%
+do({
+#browser()
+    enrResultsGrp <- .
+    ## DEBUG enrResultsGrp <- sigEnrResults
+    label <- enrResultsGrp$list_id[1]
+
+    if(nrow(enrResultsGrp)==0)  return(data.frame())
+    browser()
+
+#    warning(paste0("processing terms", paste(ac(unique(enrResultsGrp$Category)), collapse=",")))
+    logPlot <- enrResultsGrp %>%
+        ## fix factor order
+        mutate(Term=reorder(Term, -PValue) %>% reorder(as.integer(as.factor(Category)))) %>%
+        ggplot(aes(Term, PValue, fill=Category)) +
+        geom_bar(stat="identity")+
+        scale_fill_manual(values = term_category_colors, drop=F, name="Ontology") +
+        coord_flip() +
+        xlab("Enriched Terms") +
+        ggtitle(label) +
+        scale_y_log10()
+
+    fileNameLabel <- label %>%
+        str_replace_all("!=", "ne") %>%
+        str_replace_all(">", "gt") %>%
+        str_replace_all("<", "lt") %>%
+        str_replace_all(fixed("&"), "AND") %>%
+        str_replace_all(fixed("|"), "OR") %>%
+        str_replace_all(" ", "_")
+
+#        ggsave(paste0("enrichmed_terms__", fileNameLabel, ".pdf"))
+#        print(logPlot)
+
+     tmpPng <- paste0("enrterms__", fileNameLabel, ".png")
+     ggsave(tmpPng, logPlot, width=10, height = 2+round(nrow(enrResultsGrp)/5), limitsize=FALSE)
+     data.frame(file=tmpPng)
+})
+
+#+ results="asis", error=TRUE
+l_ply(term_barplot_files$file, function(pngFile){ cat(paste0("<img src='", pngFile, "'><br>"))})
+