adding DiffBind with blocking functionality (#27)

* adding DiffBind with blocking functionality

* DiffBind rule fix for blocking

workflow/Snakefile

@@ -26,6 +26,13 @@ workpath = config['project']['workpath']
 tmpdir = config['options']['tmp_dir']
 genome   = config['options']['genome']         # Reference genome of a set of samples
 assay   = config['options']['assay']
+blocks = config['project']['blocks']
+
+
+if None in list(blocks.values()):
+    blocking = False
+else:
+    blocking = True
 
 # Check for SE or PE FastQ files:
 convert = {1: False, 2: True}                     # 1 = SE, 2 = PE, -1 = Unknown
@@ -157,6 +164,7 @@ sicer_dir="sicer"
 
 uropa_dir = "UROPA_annotations"
 diffbind_dir = "DiffBind"
+diffbind_dir_block = "DiffBindBlock"
 
 if assay == "atac": 
     PeakTools = ["macsNarrow", "Genrich"] 

workflow/rules/dba.smk

@@ -163,7 +163,8 @@ rule diffbind:
         EdgeR_txt = join(workpath,diffbind_dir,"{group1}_vs_{group2}-{PeakTool}","{group1}_vs_{group2}-{PeakTool}_Diffbind_EdgeR.txt"),
         Deseq2_txt = join(workpath,diffbind_dir,"{group1}_vs_{group2}-{PeakTool}","{group1}_vs_{group2}-{PeakTool}_Diffbind_Deseq2.txt"),
         EdgeR_ftxt = join(workpath,diffbind_dir,"{group1}_vs_{group2}-{PeakTool}","{group1}_vs_{group2}-{PeakTool}_Diffbind_EdgeR_fullList.txt"),
-        Deseq2_ftxt = join(workpath,diffbind_dir,"{group1}_vs_{group2}-{PeakTool}","{group1}_vs_{group2}-{PeakTool}_Diffbind_Deseq2_fullList.txt"),  
+        Deseq2_ftxt = join(workpath,diffbind_dir,"{group1}_vs_{group2}-{PeakTool}","{group1}_vs_{group2}-{PeakTool}_Diffbind_Deseq2_fullList.txt"),
+        html_block = provided(join(workpath,diffbind_dir_block,"{group1}_vs_{group2}-{PeakTool}","{group1}_vs_{group2}-{PeakTool}_Diffbind_blocking.html"), blocking)
     params:
         rname="diffbind",
         rscript = join(workpath,"workflow","scripts","DiffBind_v2_ChIPseq.Rmd"),
@@ -175,7 +176,12 @@ rule diffbind:
         peakcaller= lambda w: FileTypesDiffBind[w.PeakTool],
         group1="{group1}",
         group2="{group2}",
-        PeakTool="{PeakTool}"
+        PeakTool="{PeakTool}",
+        blocking=blocking,
+        blocking_rscript = join(workpath,"workflow","scripts","DiffBind_v2_ChIPseq_block.Rmd"),
+        outdir_block= join(workpath,diffbind_dir_block,"{group1}_vs_{group2}-{PeakTool}"),
+        Deseq2_block = provided(join(workpath, diffbind_dir_block,"{group1}_vs_{group2}-{PeakTool}","{group1}_vs_{group2}-{PeakTool}_Diffbind_Deseq2_block.bed"), blocking),
+        EdgeR_block = provided(join(workpath, diffbind_dir_block,"{group1}_vs_{group2}-{PeakTool}","{group1}_vs_{group2}-{PeakTool}_Diffbind_EdgeR_block.bed"), blocking),
     container:
         config['images']['cfchip']
     shell: """
@@ -188,6 +194,14 @@ rule diffbind:
     params=list(csvfile= "{params.csvfile}", contrasts= "{params.contrast}", peakcaller= "{params.PeakTool}"))'
     if [ ! -f {output.Deseq2} ]; then touch {output.Deseq2}; fi
     if [ ! -f {output.EdgeR} ]; then touch {output.EdgeR}; fi
+
+    if [ '{params.blocking}' == True ]; then
+        echo "DiffBind with Blocking"
+        Rscript -e 'rmarkdown::render("{params.blocking_rscript}", output_file= "{output.html_block}", 
+        params=list(csvfile= "{params.csvfile}", contrasts= "{params.contrast}", peakcaller= "{params.PeakTool}", dir= "{params.outdir_block}"))'
+        if [ ! -f {params.Deseq2_block} ]; then touch {params.Deseq2_block}; fi
+        if [ ! -f {params.EdgeR_block} ]; then touch {params.EdgeR_block}; fi
+    fi
     """
 if assay == "cfchip":
     rule UROPA:

workflow/scripts/DiffBind_v2_ChIPseq_block.Rmd

@@ -13,6 +13,7 @@ params:
     csvfile: samplesheet.csv
     contrasts: "group1_vs_group2"
     peakcaller: "macs"
+    dir: "/path/to/DiffBindBlock/directory"
 ---
 
 <style type="text/css">
@@ -39,8 +40,10 @@ peakcaller <- params$peakcaller
 
 ```{r setup, echo=FALSE, warning=FALSE,message=FALSE}
 knitr::opts_chunk$set(echo = FALSE, include=TRUE, message=FALSE, warning=FALSE, error=FALSE)
+knitr::opts_knit$set(root.dir=params$dir)
 suppressMessages(library(DT))
 suppressMessages(library(DiffBind))
+suppressMessages(library(parallel))
 ```
 
 # Peak Data