added todos

dge_workflow/featcounts_deseq.R

@@ -2,6 +2,10 @@
 #+ include=FALSE
 
 
+## TODO Melanie: make sure reports can be spinned
+## Example: spin.R -c $NGS_TOOLS/dge_workflow/featcounts_deseq.R.R "\"count_matrix.txt""
+
+
 suppressMessages(require(docopt))
 
 doc <- '
@@ -15,9 +19,9 @@ Options:
 --min_count <min_count>   Minimal expression in any of the sample to be included in the final result list [default: 1]
 '
 
-opts <- docopt(doc, commandArgs(TRUE))
+#opts <- docopt(doc, commandArgs(TRUE))
 
-## opts <- docopt(doc, "--pcutoff 0.05 --contrasts ../time_contrasts.txt ../peak_clusters_tss5kb.count_matrix.txt")
+opts <- docopt(doc, "--pcutoff 0.05 --contrasts ../time_contrasts.txt ../peak_clusters_tss5kb.count_matrix.txt")
 ## opts <- docopt(doc, "countMatrix.txt")
 
 require(knitr)
@@ -66,9 +70,17 @@ geneInfo <- quote({
             biomaRt::getBM(mart=mart)
     }) %>% cache_it("geneInfo")
 
+########################################################################################################################
+#' # Differential Expresssion Analysis
+
+## TODO Melanie
+## * similar to example-report (see email) created with "dge_workflow/cuffdiff_report.R"
+## * sections to migrate are "Replicate Correlation and Clustering" and "Quality Control"
+
+
 
 ########################################################################################################################
-#' # Differential Binding Analysis
+#' # Differential Expresssion Analysis
 
 #' Working Directory: `r getwd()`
 
@@ -100,7 +112,7 @@ nrow(countMatrix)
 #' [Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2](http://genomebiology.com/2014/15/12/550/abstract)
 
 ## todo expose as argument
-get_sample_from_replicate <- function(repName) str_match(repName, "(.*)_[0-9]{1}")[,2]
+get_sample_from_replicate <- function(repName) str_match(repName, "(.*)_[0-9]{1}$")[,2]
 
 #' Define or load a contrasts matrix
 if(!is.null(contrasts_file)){
@@ -131,7 +143,7 @@ contrasts %>% kable()
 ## see "3.11 Sample-/gene-dependent normalization factors" in the DEseq2 manual for details
 colData <- data.frame(condition=colnames(countMatrix) %>% get_sample_from_replicate)
 dds <- DESeqDataSetFromMatrix(countMatrix, colData, formula(~ condition))
-dds <- estimateSizeFactors( dds )
+dds <- estimateSizeFactors(dds)
 
 #' #  Custom Lambda Normalization
 
@@ -317,39 +329,42 @@ degs %>%
     select(s1_over_s2_logfc, pvalue, ensembl_gene_id, sample_1, sample_2, external_gene_name, description, igv) %>%
     kable("html", table.attr = "class='dtable'", escape=F)
 
+## TODO Melanie compare with section "Differentially expressed genes" and complement missing bits
 
-########################################################################################################################
-#' ## Term enrichment Data Preparation
-#+ 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)
-}
 
-geneLists %<>% group_by(list_id)
-save(geneLists, file=".enrGeneLists.RData")
-# geneLists <- local(get(load("enrGeneLists.RData")))
-
-
-## redefine opts arguments and tweak knitr options
-opts_knit$set(root.dir = getwd())
-commandArgs <- function(x) c(paste("--overlay_expr_data ", count_matrix_file, " enrGeneLists.RData"))
-#source("/home/brandl/mnt/mack/bioinfo/scripts/ngs_tools/dev/common/david_enrichment.R")
-# child='/home/brandl/mnt/mack/bioinfo/scripts/ngs_tools/dev/common/david_enrichment.Rmd'
+## TODO later, reenable child-inclusion of enrichment analysis
+########################################################################################################################
+##' ## Term enrichment Data Preparation
+##+ 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)
+#}
+#
+#geneLists %<>% group_by(list_id)
+#save(geneLists, file=".enrGeneLists.RData")
+## geneLists <- local(get(load("enrGeneLists.RData")))
+#
+#
+### redefine opts arguments and tweak knitr options
+#opts_knit$set(root.dir = getwd())
+#commandArgs <- function(x) c(paste("--overlay_expr_data ", count_matrix_file, " enrGeneLists.RData"))
+##source("/home/brandl/mnt/mack/bioinfo/scripts/ngs_tools/dev/common/david_enrichment.R")
+## child='/home/brandl/mnt/mack/bioinfo/scripts/ngs_tools/dev/common/david_enrichment.Rmd'