cont. dge analysis

dge_workflow/dge_analysis.R

-#+ include=FALSE
-## cd ~/mnt/meninges/meninges/plus_minus
-## source <(curl https://dl.dropboxusercontent.com/u/113630701/datautils/R/utils/spinr.sh 2>&1 2>/dev/null)
-## spinr ~/mnt/meninges/from_holger/scripts/plusminus/DiffexPlusMinus.R
-
-## cd ~/mnt/meninges/meninges/plus_minus_only
-
+#!/usr/bin/env Rscript
 
+#if(!file.exists(r_script)){
+#     stop(paste("file does not exist\n", doc))
+#}
 ########################################################################################################################
-#' # Mouse aRG +/- Differential gene Expression Analysis
+#' Differential gene Expression Analysis
+#+ include=FALSE
 
-#' Alignment done already in May 2014: `/Volumes/meninges/from_holger/mouse/mouse_aRGpm`
-#' 2 Cuffdiff configratutions possilbe (just plus and minus or complete mouse data from paper)
+#' TBD: Notes
 
 ########################################################################################################################
-#+ setup2, cache=FALSE, message=FALSE, echo=FALSE
+#+ setup
 
 
 #source("http://bioconductor.org/biocLite.R")
-#biocLite("RDAVIDWebService")
+#biocLite("RDAVIDWebService", suppressUpdates=T)
 
-require(biomaRt)
-require(RDAVIDWebService)
 
 
 devtools::source_url("https://dl.dropboxusercontent.com/u/113630701/datautils/R/core_commons.R")
@@ -30,29 +25,48 @@ require(cummeRbund)
 devtools::source_url("https://dl.dropboxusercontent.com/u/113630701/datautils/R/bio/cummerutils.R")
 
 
-options(width=300) ## always overflow boxes to keep tables in shape with scrollbars
-options(java.parameters = "-Xmx4g" ) # required by read.xlsx
+#options(width=300) ## always overflow boxes to keep tables in shape with scrollbars
+#options(java.parameters = "-Xmx4g" ) # required by read.xlsx
+#require.auto(xlsx)
 
-require.auto(xlsx)
 
+suppressMessages(require.auto(docopt))
 
-#dbDir="/home/brandl/mnt/meninges/from_holger/mouse/mouse_aRGpm/cuffdiff"
-dbDir=dirname(getwd())
+doc <- '
+Convert a cuffdiff into a dge-report
+Usage: dge_analysis.R <report_name>
 
-#setwd("~/mnt/meninges/meninges/plus_minus_only")
-#dbDir="/home/brandl/mnt/meninges/from_holger/mouse/mouse_aRGpm/cuffdiff_pm"
+Options:
+--cuffdir <DIR>        Cache results [default: .]
+'
 
-#mcdir(file.path(dirname(dbDir), "dge_analysis"))
-#mcdir(file.path(dbDir, "dge_analysis"))
-#knitr::opts_knit$set(root.dir = getwd())
+#opts <- docopt(doc)
+opts <- docopt(doc, "test_report")
+
+cuffdb_directory=normalizePath(opts$cuffdir)
+
+mcdir(opts$report_name)
+
+# print(opts)
 
 
+if(is.na(file.info(cuffdb_directory)$isdir)){
+    stop(paste("db directory does not exist", doc))
+}
+
+#knitr::opts_knit$set(root.dir = getwd())
+
 #######################################################################################################################
-#' # Cuffdiff Report
+#' # Cuffdiff DB Overview
 #+ load_db, cache=FALSE
 
+#' Used cuffdiff database in `r normalizePath(opts$cuffdir)`
 
-cuff <- readCufflinks(dir=dbDir)
+## workaround until cummerbund is fixed
+tmpdir <- tempfile("cuffdb_dir")
+system(paste("cp -r", cuffdb_directory, tmpdir))
+
+cuff <- readCufflinks(dir=tmpdir)
 
 cuff
 runInfo(cuff)
@@ -61,8 +75,7 @@ replicates(cuff)
 
 
 #######################################################################################################################
-#' ## Quality Control
-
+#' ## Score Distributions and correlation
 
 # create some global summary statistic plots
 dispersionPlot(genes(cuff)) #+ aes(alpha=0.01)
@@ -165,6 +178,8 @@ ggplot(degs, aes(paste(sample_1, "vs",  sample_2), fill=sample_1_overex)) + geom
 ## todo move up into table export section
 ## add gene description
 if(!file.exists("geneInfo.RData")){
+    require(biomaRt)
+
     mousemart = useDataset("mmusculus_gene_ensembl", mart=useMart("ensembl"))
     geneInfo <- getBM(attributes=c('ensembl_gene_id', 'external_gene_name', 'description', 'gene_biotype'), mart=mousemart);
     save(geneInfo, file="geneInfo.RData")
@@ -197,6 +212,8 @@ write.delim(degs, file="degs.txt")
 #` ## Term enrichment
 #+ echo=FALSE
 
+require(RDAVIDWebService) ## just works if installed on non-network-drive (e.g. /tmp/)
+
 
 geneLists <- degs %>%  transmute(gene_id, list_id=paste(sample_1, "vs", sample_2, "ovex", sample_1_overex, sep="_"))