diff --git a/workflow/Snakefile b/workflow/Snakefile index 29164fe..c954718 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -85,7 +85,6 @@ if assay == "cfchip": join(cfTool_dir, "Output", "H3K4me3", "Signatures", "{name}.Q5DD.csv"), name=chips )) - if has_inputs: rule_all_ins.append( join(qc_dir, "AllSamples-macsNarrow", "AllSamples-macsNarrow_DiffBindQC_TMMcounts.bed") diff --git a/workflow/rules/dba.smk b/workflow/rules/dba.smk index 7958029..26f244d 100644 --- a/workflow/rules/dba.smk +++ b/workflow/rules/dba.smk @@ -44,7 +44,7 @@ if reps == "yes": otherDirs.append(diffbind_dir) mk_dir_if_not_exist(PeakTools + otherDirs) -localrules: diffbind_csv_macsN, diffbind_csv_macsB +localrules: diffbind_csv_macsN, diffbind_csv_macsB, diffbind_csv_genrich # ~~ differential binding analysis ~~ # @@ -117,6 +117,40 @@ rule diffbind_csv_macsN: )) +rule diffbind_csv_genrich: + input: + bed = expand(join(genrich_dir, "{name}", "{name}.narrowPeak"), name=chips) + output: + csvfile = join( + diffbind_dir, + "{group1}_vs_{group2}-Genrich", + "{group1}_vs_{group2}-Genrich_Diffbind_prep.csv", + ), + params: + rname = "diffbind_csv_genrich", + this_peaktool = "Genrich", + peakcaller = "narrowPeak", + this_peakextension = "_peaks.narrowPeak", + pythonscript = join(bin_path, "prep_diffbind.py"), + bam_dir = bam_dir, + workpath = workpath, + log: join(local_log_dir, "diffbind_csv_genrich", "{group1}_vs_{group2}_diffbind_csv.log") + run: + for i, contrast in enumerate(group_combos): + shell(dedent( + """ + python {params.pythonscript} \\ + --con \"""" + contrast + """\" \\ + --wp {params.workpath} \\ + --pt {params.this_peaktool} \\ + --pe {params.this_peakextension} \\ + --bd {params.bam_dir} \\ + --pc {params.peakcaller} \\ + --csv {output.csvfile} + """ + )) + + rule diffbind_csv_macsB: input: beds = expand(join(macsB_dir, "{name}", "{name}_peaks.broadPeak"), name=chips), diff --git a/workflow/rules/qc.smk b/workflow/rules/qc.smk index ad003de..317c8af 100644 --- a/workflow/rules/qc.smk +++ b/workflow/rules/qc.smk @@ -426,6 +426,37 @@ rule FRiP_macsN: """ +rule FRiP_Genrich: + input: + bed = expand(join(genrich_dir, "{name}", "{name}.narrowPeak"), name=chips), + bam = join(bam_dir, "{name}.Q5DD.bam"), + output: + join(workpath, "PeakQC", "Genrich.{name}.Q5DD.FRiP_table.txt"), + params: + rname = "FRiP_macsN", + outroot = join(peakqc_dir, "macsNarrow"), + script = join(bin_path, "frip.py"), + genome = config['references'][genome]['REFLEN'], + tmpdir = tmpdir, + container: + config['images']['python'] + shell: + """ + # Setups temporary directory for + # intermediate files with built-in + # mechanism for deletion on exit + if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi + tmp=$(mktemp -d -p "{params.tmpdir}") + trap 'rm -rf "${{tmp}}"' EXIT + + python {params.script} \\ + -p {input.bed} \\ + -b {input.bam} \\ + -g {params.genome} \\ + -o {params.outroot} + """ + + rule FRiP_macsB: input: bed = expand(join(macsB_dir, "{name}", "{name}_peaks.broadPeak"), name=chips),