started gsea enrichment tool

common/gsea_enrichment.R

+#!/usr/bin/env Rscript
+#+ echo=FALSE, error=F
+
+
+# cd
+
+suppressMessages(require(docopt))
+
+## todo use textual input here for ease of use
+doc <- '
+Perform an enrichment analysis for a set of genes
+Usage:gsea_enrichment.R[options] <sorted_gene_lists_tsv> <group_col>
+
+Options:
+--list_id_col Column containing the grouping variable to speparate different lists [default: ]
+--project <project_prefix>   Name to prefix all generated result files [default: ]
+--qcutoff <qcutoff>             Use a q-value cutoff of 0.01 instead of a q-value cutoff [default: 0.01]
+
+'
+
+opts <- docopt(doc, commandArgs(TRUE)) ## does not work when spining
+# opts <- docopt(doc, "--overlay_expr_data ../plot_score_matrix.txt ../degs_by_contrast.txt contrast" )
+
+devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.27/R/core_commons.R")
+devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.27/R/ggplot_commons.R")
+devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.27/R/bio/diffex_commons.R")
+
+loadpack(knitr)
+loadpack(DT)
+
+#loadpack(clusterProfiler)
+
+#devtools::session_info() # nice!
+
+## to fix child support issue with knitr, see also
+## http://stackoverflow.com/questions/20030523/knitr-nested-child-documents
+## https://github.com/yihui/knitr/issues/38
+# todo disabled because root.dir in parent document seems the only working solution
+#if(exists('project_dir')) setwd(project_dir)
+#print(getwd())
+
+## load the data
+geneLists <- read_tsv(opts$sorted_gene_lists_tsv)
+group_col=opts$group_col
+geneLists %<>% group_by_(.dots=group_col)
+
+#geneLists %<>% filter(cluster %in% c("cluster_1", "cluster_2"))
+
+## 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 <- if(str_length(opts$project)>0) paste0(opts$project, ".") else basename(opts$sorted_gene_lists_tsv) %>% trim_ext("txt") #%>% paste0(".")
+#resultsBaseName=basename(opts$gene_lists_tsv) %>% trim_ext("txt") #%>% paste0(".")
+
+qCutoff <- as.numeric(opts$qcutoff)
+
+
+#reload_dplyr()
+
+########################################################################################################################
+#' ## Enrichment Analysis
+
+#' Run configuration was
+vec2df(unlist(opts)) %>% filter(!str_detect(name, "^[<-]")) %>% kable()
+
+#' This analysis was performed using [clusterProfiler](http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html). The following ontologies were tested: Kegg, Go, Reactome, Dose,
+
+listLabels <- geneLists %>% select(-ensembl_gene_id) %>% distinct
+listLabels %<>% transform(list_label=do.call(paste, c(listLabels, sep="__")))
+
+geneLists %>% inner_join(listLabels) %>%
+    ggplot(aes(list_label)) +
+    geom_bar() +
+    coord_flip() +
+    ggtitle("gene list sizes to be tested for term enrichment") +
+    ylab("")
+
+
+## todo move to diffex commons
+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("http://bioconductor.org/biocLite.R")
+#biocLite("org.Mm.eg.db")
+#biocLite("org.Hs.eg.db")
+#biocLite("org.Dr.eg.db")
+#biocLite("org.Dm.eg.db")
+#biocLite("KEGG.db")
+
+
+#data(gcSample)
+#yy = enrichKEGG(gcSample[[5]], pvalueCutoff=0.01)
+#head(summary(yy))
+#plot(yy)
+## seems broken
+
+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
+
+glMappedRaw <-  clusterProfiler::bitr(geneLists$ensembl_gene_id, fromType="ENSEMBL", toType="ENTREZID", OrgDb=annoDb) %>%
+    set_names("ensembl_gene_id", "entrez_gene_id") %>% right_join(geneLists)
+
+
+
+#' Check how many failed to map
+count(glMappedRaw, is.na(entrez_gene_id))
+
+unloadNamespace('clusterProfiler')
+#loadpack(clusterProfiler)
+
+reload_dplyr()
+
+glMapped <- glMappedRaw %>% filter(!is.na(entrez_gene_id)) %>% select(-ensembl_gene_id)  %>%
+    ## regroup the data
+    #    group_by_(cluster)
+    group_by_(.dots=group_col %>% ac)
+
+
+## retain just 1500 genes at max per group
+#glMappedSub <- glMapped %>% do({shuffle(.) %>% head(1500)})
+#count(glMappedSub, cluster)
+
+## sync reactome pacakge to node
+#if(Sys.getenv("LSF_SERVERDIR")!="" && dir.exists("/tmp/local_r_packages")){
+#    system("if [ ! -d '/tmp/local_r_packages/reactome.db/' ]; then scp -r falcon:/tmp/local_r_packages /tmp ; fi")
+#}
+
+loadpack(ReactomePA)
+
+cp_test <- function(geneIds){
+    # 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
+
+    if(length(geneIds)>1500){
+        geneIds <- sample(geneIds) %>% head(1500)
+    }
+
+    echo("testing", length(geneIds), " genes for enrichment")
+
+#    PANTHER10_ontology <- read.delim("http://data.pantherdb.org/PANTHER10.0/ontology/Protein_Class_7.0")
+
+#    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()
+
+    #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")
+
+    enrResults <- rbind_list(
+        mutate(keggResults, ontology="kegg"),
+        mutate(reactomeResults, ontology="reactome"),
+        mutate(goResultsBP, ontology="go_bp"),
+        mutate(goResultsMF, ontology="go_mf"),
+        mutate(goResultsCC, ontology="go_cc")
+    )
+    enrResults
+#    echo("numResults", nrow(enrResults))
+}
+
+## prety slow
+enrResults <-  quote(glMapped %>% do(cp_test(.$entrez_gene_id))) %>% cache_it(paste0("enrdata_", digest(glMapped)))
+
+## run the actual enrichment test for all clusters and all ontologies
+#library(foreach); library(doMC); registerDoMC(cores=20)
+## https://support.bioconductor.org/p/38541/
+## library(GOstats) inside your %dopar% loop. And then start with a fresh session, so GOstats is not loaded outside the loop. It worked for me.
+#enrResults <- dlply(glMapped, groups(glMapped) %>% ac, function(curGroup){
+#  cp_test(curGroup$entrez_gene_id)
+#}, .progress="text", .parallel=F) ##%>% cache_it(paste0("enrdata_", digest(glMapped)))
+#rbind_all(enrResults)
+
+#reload_dplyr()
+
+#enrResults %<>% rename(Category=Description)
+#enrResults %<>% rename(ontology=Category)
+
+## remove to clumsy gene_id columns
+enrResults %<>% select(-geneID)
+
+write_tsv(enrResults, path=paste0(resultsBaseName, "enrResults.txt"))
+# enrResults <- read.delim(paste0(resultsBaseName, "enrResults.txt"))
+#'  [Enrichment Results](`r paste0(resultsBaseName, "enrResults.txt")`)
+
+loadpack(DT)
+datatable(enrResults)
+
+#enrResults %>% ggplot(aes(pvalue)) + geom_histogram() + scale_x_log10()
+#enrResults %>% ggplot(aes(ontology)) + facet_wrap(~cluster) + geom_bar() + rot_x_lab()
+#facetSpecs <- paste("~", groups(geneLists) %>% ac %>% paste(collapse=" + "))
+facetSpecs <- paste("~", group_col %>% ac %>% paste(collapse=" + "))
+
+#' Visualize term-pvalues per list
+
+#http://stackoverflow.com/questions/11028353/passing-string-variable-facet-wrap-in-ggplot-using-r
+enrResults %>% ggplot(aes(ontology)) + facet_wrap(as.formula(facetSpecs), ncol=3) + geom_bar() + rot_x_lab() + ggtitle("enriched term counts by cluster")
+
+enrResults %<>% mutate(num_term_genes=str_split_fixed(BgRatio, fixed("/"),2)[,1] %>% as.numeric)
+
+#' Keep at max 100 terms for visualzation per group
+#if(table(enrResults$cluster))
+## restablish the grouping and limit results per group to 100
+#enrResults %<>% group_by_(.dots=groups(geneLists) %>% ac)
+#enrResults %<>% group_by_(.dots=group_col)
+#enrResults %<>% sample_frac(1) %>% filter((row_number()<100)) %>% group_by_(.dots=group_col)
+
+#tt <- enrResults %>% group_by_(.dots=c(group_col,"ontology")) %>% arrange(-qvalue) %>% dplyr::slice(1:10)
+#enrResults %>% filter(cluster=="cluster_1")
+#enrResults %>% group_by_(.dots=c(group_col)) %>% ggplot(aes(pvalue, qvalue, color=ontology))+ geom_point()
+#enrResults %>% group_by_(.dots=c(group_col)) %>% ggplot(aes(pvalue, p.adjust, color=cluster))+ geom_point()
+erPlotData <- enrResults %>% group_by_(.dots=c(group_col)) %>% arrange(qvalue) %>% slice(1:15) %>%
+      ## regroup because otherwise dplyr complains about corrupt df (which looks like a bug)
+      group_by_(.dots=c(group_col))
+
+#count(erPlotData, cluster)
+
+
+
+#+ include=FALSE, eval=FALSE
+## 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)
+#})
+
+#sigEnrResults %>%
+##    select(-Genes) %>%
+#    do({
+#    enrResultsGrp <- .
+#    ## DEBUG enrResultsGrp <- sigEnrResults
+##    geneLists %>% first(1) %>% select(-ensembl_gene_id) %>% paste0(., collapse="_")
+##browser()
+#    label <- geneLists %>% semi_join(enrResultsGrp) %>% first(1) %>% dplyr::select(-ensembl_gene_id) %>% paste0(., collapse="_")
+#
+#    echo("processing", label)
+#
+#    logPlot <- enrResultsGrp %>%
+#        ## fix factor order
+#        mutate(Term=reorder(Term, -PValue) %>% reorder(as.integer(Category))) %>%
+#        ggplot(aes(Term, PValue, fill=Category)) +
+#	    geom_bar(stat="identity")+
+#	    coord_flip() +
+#	    xlab("Enriched Terms") +
+#	    ggtitle(label) +
+#	    scale_y_log10()
+#
+#	    ggsave(paste0(resultsBaseName, label, ".enrichmed_terms.pdf"))
+#	    print(logPlot)
+#})
+##ggsave2()
+
+
+#+ error=TRUE, echo=FALSE
+warning("dropping levels")
+erPlotData %<>% mutate(ontology=ac(ontology)) ## drop unsused level to get consistent color palette
+
+erPlotData %<>% rename(Term=Description)
+term_category_colors <- create_palette(unique(ac(erPlotData$ontology)))
+
+figDir <- "enr_charts"
+dir.create(figDir)
+
+term_barplot_files <- erPlotData %>%
+    ## chop and pad category names
+    mutate(Term=str_sub(Term, 1, 70) %>% str_pad(70)) %>%
+do({
+    enrResultsGrp <- .
+
+    ## DEBUG   enrResultsGrp <- erPlotData %>% ungroup() %>% filter(contrast==contrast[1])
+#    with(enrResultsGrp, as.data.frame(table(contrast)))
+
+#    browser()
+    ## DEBUG enrResultsGrp <- sigEnrResults
+#    label <- geneLists %>% semi_join(enrResultsGrp) %>% first(1) %>% dplyr::select(-ensembl_gene_id) %>% paste0(., collapse="_")
+    label <- subset(enrResultsGrp, select=group_col)[1,1] %>% as.matrix %>% ac
+
+
+#stopifnot(all(!is.na(enrResultsGrp$Term)))
+#stopifnot(all(!is.na(enrResultsGrp$ontology)))
+#stopifnot(all(!is.na(enrResultsGrp$num_term_genes)))
+
+    #browser()
+#    warning(paste0("processing terms", paste(ac(unique(enrResultsGrp$ontology)), collapse=",")))
+    logPlot <- enrResultsGrp %>%
+        ## fix factor order
+#        mutate(Term=reorder(Term, -qvalue) %>% reorder(as.integer(as.factor(ontology)))) %>%
+        mutate(Term=reorder(Term, -qvalue)) %>%
+        ggplot(aes(Term, num_term_genes, fill=-log10(qvalue), color=ontology)) +
+        geom_bar(stat="identity") +
+        scale_color_manual(values = term_category_colors, drop=F, name="Ontology") +
+        coord_flip() +
+        xlab("Enriched Terms") +
+        ggtitle(label) +
+        scale_y_log10()
+
+        # print(logPlot)
+
+    ## todo  use builtin method to create filesystem-compatible name
+    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)
+
+     stopifnot(str_length(fileNameLabel)>0)
+     tmpPng <- paste0(figDir, "/enrterms__", fileNameLabel, ".png")
+     ggsave(tmpPng, logPlot, width=10, height = 2+round(nrow(enrResultsGrp)/5), limitsize=FALSE)
+     data.frame(file=tmpPng)
+})
+
+#+ results="asis"
+l_ply(term_barplot_files$file, function(pngFile){ cat(paste0("<img src='", pngFile, "'><br>"))})