Commit 66ce867a authored by Holger Brandl's avatar Holger Brandl

improved enrichement report

parent 58fc70b5
......@@ -20,12 +20,13 @@ Options:
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" )
coreCommons = "https://raw.githubusercontent.com/holgerbrandl/datautils/v1.39/R/core_commons.R"
coreCommons = "https://raw.githubusercontent.com/holgerbrandl/datautils/v1.40/R/core_commons.R"
devtools::source_url(coreCommons)
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.39/R/ggplot_commons.R")
# devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.39/R/bio/cp_utils.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.40/R/ggplot_commons.R")
# devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.40/R/bio/cp_utils.R")
devtools::source_url("https://www.dropbox.com/s/p2b8luxf7jteb63/cp_utils.R?dl=1")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.39/R/bio/diffex_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.40/R/bio/diffex_commons.R")
load_pack(knitr)
load_pack(DT)
......@@ -54,7 +55,8 @@ if (! is.null(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 <- if (str_length(opts$project) > 0) paste0(opts$project, ".") else basename(opts$gene_lists_tsv) %>% trim_ext("txt") #%>% paste0(".")
# resultsBaseName <- if (str_length(opts$project) > 0) paste0(opts$project, ".") else basename(opts$gene_lists_tsv) %>% trim_ext("txt") #%>% paste0(".")
resultsBaseName <- if (str_length(opts$project) > 0) paste0(opts$project, ".") else ""
#resultsBaseName=basename(opts$gene_lists_tsv) %>% trim_ext("txt") #%>% paste0(".")
q_cutoff <- as.numeric(opts$qcutoff)
......@@ -187,7 +189,7 @@ enrResults %<>% mutate(num_term_genes = str_split_fixed(BgRatio, fixed("/"), 2)[
erPlotData <- enrResults %>%
group_by_(.dots = c(group_col)) %>%
arrange(qvalue) %>%
slice(1 : 15) %>%
slice(1 : 20) %>%
## regroup because otherwise dplyr complains about corrupt df (which looks like a bug)
group_by_(.dots = c(group_col))
......@@ -310,11 +312,17 @@ term_barplot_files = erPlotData %>% do({
#+ results="asis"
walk(term_barplot_files$file, function(pngFile){ cat(paste0("<img src='", pngFile, "'><br>"))})
# install.packages("session")
session::save.session(".cp_enrichment.dat");
# session::restore.session(".cp_enrichment.dat");
########################################################################################################################
#' ## Enriched KEGG Pathways
#' To understand spatio-temporal changes in gene expression better we now overlay enriched kegg pathways with the -log10(q_value) of each contrast. The direction of the expression changes is encoded as color, whereby red indicates that sample_1 is overexpressed. Because we have multiple contrasts of interest, this defines a slice-color barcode for each gene. To relate the barcode to contrasts we define the following slice order:
#+ eval=(!is.null(opts$overlay_expr_data) & (nrow(enrResults %>% filter(ontology=="kegg")) >0))
hasEnrPathways=(!is.null(opts$overlay_expr_data) & (nrow(enrResults %>% filter(ontology=="kegg")) >0))
#+ eval=hasEnrPathways
## todo why is tidyr not processing an empty dataframe
keggPathways = enrResults %>% filter(ontology == "kegg") %$% ac(ID) %>% unique() # %>% head(3) ## todo remove debug head
......@@ -324,12 +332,12 @@ keggPathways = enrResults %>% filter(ontology == "kegg") %$% ac(ID) %>% unique()
# keggPathways = c("mmu04660", "mmu03060", "mmu04110", "mmu04260", "mmu04723", "mmu05322", "mmu04114")
keggPathways %<>% purrr::discard(~ . == "mmu04723") ## overloaded by scales packages
## #+ results='asis', echo=FALSE
if (! exists("keggPathways") || length(keggPathways) == 0) {
cat("No enriched pathways found")
return()
quit()
}
# ## #+ results='asis', echo=FALSE
# if (! exists("keggPathways") || length(keggPathways) == 0) {
# cat("No enriched pathways found")
# return()
# quit()
# }
......@@ -377,10 +385,6 @@ dir.create(pwPlotDir)
# set_names("entrez_gene_id", "ensembl_gene_id" ) %>%
# inner_join( sliceData %>% add_rownames("ensembl_gene_id") )
# install.packages("session")
session::save.session(".cp_enrichment.dat");
# session::restore.session(".cp_enrichment.debug.dat");
## from http://stackoverflow.com/questions/7505547/detach-all-packages-while-working-in-r
unload_packages()
......@@ -434,7 +438,7 @@ plot_pathway <- function(pathwayID, sliceData){
capture.output
#+ echo=FALSE, warning=FALSE
#+ echo=FALSE, warning=FALSE, eval=hasEnrPathways
pathwayPlots <- keggPathways %>%
# head(3) %>%
map(function(pathwayID){
......@@ -483,7 +487,7 @@ makeTooltip <- function(entrez_id){
#+ results="asis", echo=FALSE
#+ results="asis", echo=FALSE, eval=hasEnrPathways
## simple non-clickable plots
## http://stackoverflow.com/questions/12588323/r-how-to-extract-values-for-the-same-named-elements-across-multiple-objects-of
......@@ -527,7 +531,7 @@ createImgMap <- function(plotData){
mutate(x = x - 22, y = y - 8) %>% ## todo does not work for non-gene elements
mutate(link = paste0("http://www.kegg.jp/dbget-bin/www_bget?", keggOrCode, ":", kegg.names)) %>%
rowwise() %>%
do({ curNode = .; mutate(as.df(curNode), tooltip = makeTooltip(as.integer(curNode$kegg.names)))})
do({ curNode = .; mutate(as_df(curNode), tooltip = makeTooltip(as.integer(curNode$kegg.names)))})
## create tooltip by using mapping to
......@@ -558,11 +562,9 @@ pathwayPlots %>% plyr::l_ply(createImgMap)
# }
#' ## Notes
#' Other interesting tools to perform enrichment testing
#' * [piano](http://www.bioconductor.org/packages/release/bioc/html/piano.html)
#' System info:
#+ results="asis"
#+
devtools::session_info()
\ No newline at end of file
......@@ -10,6 +10,10 @@ Other enrichment tools
* http://amp.pharm.mssm.edu/Harmonizome/
Other interesting tools to perform enrichment testing
-----------------------------------------------------
* [piano](http://www.bioconductor.org/packages/release/bioc/html/piano.html)
IHW Eval
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment