moved to common tools

chipseq_workflow/binding_profiles.R

+#!/usr/bin/env Rscript
+#+ echo=FALSE, message=FALSE
+
+suppressMessages(library(docopt))
+
+doc <- '
+Create a binding profile plot with ngs-plot
+Usage: binding_profiles.R [options] <comma_sep_bamfiles> <comma_sep_features>
+
+Options:
+--lists <comma_sep_genelists>    genelists to subset by
+--F     Propagated to ngsplot as -F (see their manual)
+--listslabel <listslabel>    Optional label for the supplied gene lists
+'
+
+#todo support control as defined in https://github.com/shenlab-sinai/ngsplot#examples
+
+
+opts <- docopt(doc, commandArgs(TRUE))
+
+## examples
+#opts <- docopt(doc, "--lists /projects/bioinfo/holger/projects/khan_chipseq_h1/resources/maternal_highexpr_giraldez.lst,/projects/bioinfo/holger/projects/khan_chipseq_h1/resources/maternal_lowexpr_giraldez.lst /projects/bioinfo/holger/projects/khan_chipseq_h1/alignments_mmfilt/H1M_Dome_mmf.bam 'tss,genebodies,/projects/bioinfo/holger/projects/khan_chipseq_h1/ngsplot/exons_all.bed'")
+#opts <- docopt(doc, "/projects/bioinfo/holger/projects/khan_chipseq_h1/alignments_mmfilt/H1M_Dome_mmf.bam 'tss,genebodies,/projects/bioinfo/holger/projects/khan_chipseq_h1/ngsplot/exons_all.bed'")
+
+devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.9/R/core_commons.R")
+
+require.auto(whisker) # http://stackoverflow.com/questions/10341114/alternative-function-to-paste
+
+
+features <- str_split(opts$comma_sep_features, ",") %>% unlist()
+bams <- str_split(opts$comma_sep_bamfiles, ",") %>% unlist()
+gene_lists <- str_split(opts$lists, ",") %>% unlist()
+
+stopifnot(gene_lists %>% sapply(file.exists) %>% all())
+stopifnot(bams %>% sapply(file.exists) %>% all())
+
+
+simplify_bam <- function(bamFile) basename(ac(bamFile)) %>% str_replace("_mmf.bam|.bam", "")
+simplify_group <- function(geneList) basename(ac(geneList)) %>% str_replace(".lst", "")
+simplify_bed <- function(bedFile) basename(ac(bedFile)) %>% str_replace(".bed", "")
+
+bamCombi=bams %>% laply(simplify_bam) %>% paste(collapse="__")
+
+listsLabel=as.character(opts$listslabel)
+
+## use concatenated lists as fallback label
+if(is.null(listsLabel)){
+    listsLabel=gene_lists %>% laply(simplify_group) %>% paste(collapse="__")
+}
+
+########################################################################################################################
+#' Handle default feature types
+
+
+#' Write a generic driver
+npConfig <- expand.grid(
+        bam_file=bams,
+        gene_list=if(is.null(gene_lists)) "." else c("-1", gene_lists)
+    ) %>%
+    ## define label
+    mutate(
+        label=if(length(bams)>1) paste0(simplify_bam(bam_file)) else "",
+        label=if(!is.null(gene_lists)) paste0(label, ifelse(label=="", "","."),  simplify_group(gene_list)),
+        label=ifelse(label=="-1", "all", label)
+    )
+
+
+driver_file=paste(bamCombi, listsLabel, "npd", sep=".")
+npConfig %>% write.delim(driver_file, header=F, quote=F)
+
+
+
+default_profile <- function(featureType){
+## DEBUG featureType="tss"
+
+stopifnot( featureType %in% c("tss", "tes", "genebody", "exon", "cgi", "enhancer", "dhs"))
+
+#system("bsub 'sleep 10; echo hello > test.txt' | joblist")
+#whisker.render ( "Hi {{name}}, ${other} Have a very nice {{noun}} ! " , list(name="Tom", noun="day") )
+#system("echo ${project}_ngpsplot")
+
+## submit the job
+'
+driver={{driver}}
+feature={{feature}}
+listsLabel={{listsLabel}}
+jobname=$(basename $driver .npd).${feature}
+mysub "${project}__ngsplot__${jobname}" "
+    ngs.plot.r -G Zv9 -R ${feature} -C ${driver} -O ${jobname} -L 2000
+" -q medium | joblist .ngsplot
+
+' %>% whisker.render(list(driver=driver_file, feature=featureType, listsLabel=listsLabel)) %T>% cat() %>% system()
+}
+
+
+features %>% subset(., !str_detect(., "[.]bed")) %>% ldply(default_profile)
+
+
+########################################################################################################################
+## Handle the more complicated bed files features
+
+
+bedFeatures <- features %>% subset(., str_detect(., "[.]bed"))
+
+stopifnot(bedFeatures %>% sapply(file.exists) %>% all())
+
+bed_profile <- function(bedFile){
+    ## DEBUG bedFile<-"/projects/bioinfo/holger/projects/khan_chipseq_h1/ngsplot/exons_all.bed"
+
+    echo("processing {{bedFile}}" %>% whisker.render)
+
+    listBeds <- llply(gene_lists, function(geneList){
+        bedFile <- bedFile
+        ## geneList <- gene_lists[1]
+        whisker.render('
+        geneList={{geneList}}
+        bedFile={{bedFile}}
+        listBed=$(basename $bedFile .bed).$(basename $geneList .lst).bed
+
+        grep -Ff $geneList $bedFile > $listBed
+        echo $listBed
+        ') %T>% cat %>% system(intern=T)
+    }) %>% unlist()
+
+    ## readd all as reference
+    listBeds <- c(listBeds, bedFile)
+
+    #' bed specific driver
+    npConfig <- expand.grid(
+            bam_file=bams,
+            filt_bed_file=listBeds
+        ) %>%
+        ## define label
+        mutate(
+                bed_grp=simplify_bed(filt_bed_file) %>% str_replace(paste0(simplify_bed(bedFile), "."), ""),
+                label=if(length(bams)>1) paste0(simplify_bam(bam_file)) else "",
+                label=if(!is.null(listBeds)) paste0(label, ifelse(label=="", "","."), bed_grp)
+        ) %>% select(-bed_grp)
+
+    driver_file=paste(bamCombi, simplify_bed(bedFile), listsLabel, "npd", sep=".")
+    npConfig %>% write.delim(driver_file, header=F, quote=F)
+
+    '
+    driver={{driver}}
+    listsLabel={{listsLabel}}
+    jobname=$(basename ${driver} .npd).${listsLabel}
+    mysub "${project}__ngsplot__${jobname}" "
+        ngs.plot.r -G Zv9 -R bed -C ${driver} -O ${jobname}   -L 2000
+    " -q medium | joblist .ngsplot
+    ' %>% whisker.render(list(driver=driver_file, listsLabel=listsLabel)) %T>% cat()  %>% system()
+}
+
+
+features %>% subset(., str_detect(., "[.]bed")) %>% ldply(bed_profile)
+