cont. dge analysis for stefania

dge_workflow/dge_analysis.R

@@ -19,7 +19,7 @@ Options:
 # setwd("/Volumes/project-ste/rnaseq_isnm1/stefania_isnm1_2nd_batch/cuffdiff_PC/dge_report")
 
 #opts <- docopt(doc, commandArgs(TRUE)) ## does not work when spining
-# Stefania: opts <- docopt(doc, "--undirected --qcutoff 0.05 --minfpkm 2  ..")
+# Stefania: opts <- docopt(doc, "--undirected --qcutoff 0.05 --minfpkm 2 ..")
 # opts <- docopt(doc, "--undirected  ..")
 # opts <- docopt(doc, paste(Sys.getenv("DGE_PARAMS"), ".."))
 ## todo use commandArgs here!!
@@ -44,9 +44,9 @@ if(is.na(file.info(cuffdb_directory)$isdir)){
 
 devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.4/R/core_commons.R")
 devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.4/R/ggplot_commons.R")
+devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.4/R/bio/diffex_commons.R")
 
 require(cummeRbund)
-devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.4/R/bio/diffex_commons.R")
 
 require(knitr)
 require(tidyr)
@@ -56,8 +56,8 @@ require(tidyr)
 ########################################################################################################################
 #' # Differential Gene Expression Analysis
 
-print("configuration")
-vec2df(unlist(opts)) %>% kable()
+print("Configuration")
+vec2df(unlist(opts)) %>% filter(!str_detect(name, "^[<-]")) %>% kable()
 
 #' [CuffDiff](http://cufflinks.cbcb.umd.edu/manual.html) and [cummeRbund](http://compbio.mit.edu/cummeRbund/) were used to perform this analysis. For details about how the test for differntial expression works, refer to [Differential analysis of gene regulation at transcript resolution with RNA-seq](http://www.nature.com/nbt/journal/v31/n1/full/nbt.2450.html).
 
@@ -67,11 +67,18 @@ vec2df(unlist(opts)) %>% kable()
 
 #+ load_db, cache=FALSE
 
-## workaround until sqlite  is fixed on lustre
-tmpdir <- tempfile("cuffdb_dir")
-system(paste("cp -r", cuffdb_directory, tmpdir))
+isCluster=Sys.getenv("LSF_SERVERDIR")!=""
+dbTempCopyRequired=isCluster | str_detect(normalizePath(cuffdb_directory), fixed("/Volumes/projects"))
+
+if(dbTempCopyRequired){
+    ## workaround until sqlite  is fixed on lustre
+    tmpdir <- tempfile("cuffdb_dir")
+    system(paste("cp -r", cuffdb_directory, tmpdir))
+    cuff <- readCufflinks(dir=tmpdir)
+}else{
+    cuff <- readCufflinks(dir="..")
+}
 
-cuff <- readCufflinks(dir=tmpdir)
 
 cuff
 runInfo(cuff)
@@ -141,7 +148,6 @@ repCountMatrix(genes(cuff)) %>%
 
 #######################################################################################################################
 #' ## Score Distributions
-#+ include=FALSE
 
 
 #if(!as.logical(opts$S)){
@@ -166,6 +172,8 @@ fpkmSCVPlot(genes(cuff))
 
 #csVolcano(genes(cuff),"aRGm", "N") + ggtitle("odd line-shaped fan-out in volcano plot")
 
+## todo consider to add more pretty volcano plots (see stefania's extra plots
+
 
 #csVolcano(genes(cuff),"a", "b")
 ## todo replace or complement with marta plot
@@ -254,9 +262,6 @@ allDiff <- diffData(genes(cuff)) %>%
 
 #getExpressedGenes(cuff, minFPKM=minFPKM) %>% writeLines(paste0("genes_expressed_gt_", minFPKM[1], "_in_at_least_one_sample.txt"))
 
-## save as reference
-merge(allDiff, geneInfo, by="ensembl_gene_id", all.x=T) %>% write.delim(file="diffData.txt")
-
 # diffData <- read.delim("diffData.txt")
 #'  [diffData](diffData.txt)
 
@@ -271,6 +276,9 @@ if(!is.null(qcutoff)){
     allDiff <- transform(allDiff, isHit=p_value<=pcutoff)
 }
 
+## save as reference
+merge(allDiff, geneInfo, by="ensembl_gene_id", all.x=T) %>% write.delim(file="diffData.txt")
+
 
 #+
 # #' Used cutoff criterion was: q_value<0.01
@@ -281,6 +289,8 @@ degs <- subset(allDiff, isHit)
 #degs <- subset(allDiff, significant=="yes")
 
 
+## todo embedd as paginated DataTable
+
 with(degs, as.data.frame(table(sample_1, sample_2, sample_1_overex))) %>% filter(Freq >0)
 
 ## add gene info
@@ -463,13 +473,16 @@ save(pathwayPlots, file=".pathwayPlots.RData")
 #unlist(lapply(pathwayPlots, "[[", "plotfile"))
 # unlist(lapply(pathwayPlots, "[[", "plotfile")) %>% l_ply(function(pngFile){ cat(paste0("<img src='", pngFile, "'><br>"))})
 
-cat("
-<style type='text/css'>
- img {
-   max-width: 100%;
- }
-</style>
-")
+#cat("
+#<style type='text/css'>
+# img {
+#   max-width: 100%;
+# }
+#</style>
+#")
+
+## todo add tooltips with additinal info
+## http://stackoverflow.com/questions/5478940/how-to-create-a-html-tooltip-on-top-of-image-map
 
 ## extended version with clickable links
 pathwayPlots %>% l_ply(function(plotData){
@@ -486,7 +499,8 @@ pathwayPlots %>% l_ply(function(plotData){
 #    "http://www.kegg.jp/dbget-bin/www_bget?mmu:16412
 
     ## first the image itself
-    cat(paste0("<img usemap='#", pathway_id,"' src='", pngFile, "'><br>"))
+#    cat(paste0("<img usemap='#", pathway_id,"' src='", pngFile, "'><br>"))
+    cat(paste0("<div style='width: 2000px'><img style='float: left' usemap='#", pathway_id,"' src='", pngFile, "'></div><br>"))
     cat(paste0("<map name='", pathway_id,"'>"))
     #keggNodes %>% rowwise() %>% {curNode=.; cat(curNode$name)}