From a7a5ae35731e22e1f9d0e039661ac1fef9b6398b Mon Sep 17 00:00:00 2001
From: Asya Khleborodova <asyakhl@users.noreply.github.com>
Date: Tue, 12 Dec 2023 13:29:06 -0500
Subject: [PATCH] adding DiffBind with blocking functionality (#27)

* adding DiffBind with blocking functionality

* DiffBind rule fix for blocking
---
 workflow/Snakefile                             |  8 ++++++++
 workflow/rules/dba.smk                         | 18 ++++++++++++++++--
 workflow/scripts/DiffBind_v2_ChIPseq_block.Rmd |  3 +++
 3 files changed, 27 insertions(+), 2 deletions(-)

diff --git a/workflow/Snakefile b/workflow/Snakefile
index 4a8966c..fc89e81 100644
--- a/workflow/Snakefile
+++ b/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"] 
diff --git a/workflow/rules/dba.smk b/workflow/rules/dba.smk
index 0cb5621..f934820 100644
--- a/workflow/rules/dba.smk
+++ b/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:
diff --git a/workflow/scripts/DiffBind_v2_ChIPseq_block.Rmd b/workflow/scripts/DiffBind_v2_ChIPseq_block.Rmd
index 1866add..255408a 100644
--- a/workflow/scripts/DiffBind_v2_ChIPseq_block.Rmd
+++ b/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