diff --git a/workflow/rules/paired/peakcall.smk b/workflow/rules/paired/peakcall.smk index 3c1f0bc..8f9ba6a 100644 --- a/workflow/rules/paired/peakcall.smk +++ b/workflow/rules/paired/peakcall.smk @@ -12,7 +12,6 @@ tmpdir = config['options']['tmp_dir'] rule MACS2_broad: input: chip = join(bam_dir, "{name}.Q5DD.bam"), - txt = join(ppqt_dir, "{name}.Q5DD.ppqt.txt"), c_option = lambda w: join(bam_dir, f"{chip2input[w.name]}.Q5DD.bam") if chip2input[w.name] else [], output: @@ -37,6 +36,32 @@ rule MACS2_broad: """ +rule MACS2_narrow: + input: + chip = join(bam_dir, "{name}.Q5DD.bam"), + c_option = lambda w: join(bam_dir, f"{chip2input[w.name]}.Q5DD.bam") + if chip2input[w.name] else [], + output: + join(macsN_dir, "{name}", "{name}_peaks.narrowPeak"), + params: + rname = 'MACS2_narrow', + gsize = config['references'][genome]['EFFECTIVEGENOMESIZE'], + macsver = config['tools']['MACSVER'], + flag = lambda w: "-c" if chip2input[w.name] else "", + shell: + """ + module load {params.macsver}; + macs2 callpeak \\ + -t {input.chip} {params.flag} {input.c_option} \\ + -g {params.gsize} \\ + -n {wildcards.name} \\ + --outdir {macsN_dir}/{wildcards.name} \\ + -q 0.01 \\ + --keep-dup="all" \\ + -f "BAMPE" + """ + + rule SICER: input: chip = lambda w: join(bam_dir, w.name + ".Q5DD.bam"), @@ -111,31 +136,4 @@ rule SICER: mv ${{tmp}}/{params.name}.Q5DD-W300-G600.scoreisland {params.this_sicer_dir} fi - """ - - -rule MACS2_narrow: - input: - chip = join(bam_dir, "{name}.Q5DD.bam"), - txt = join(ppqt_dir, "{name}.Q5DD.ppqt.txt"), - c_option = lambda w: join(bam_dir, f"{chip2input[w.name]}.Q5DD.bam") - if chip2input[w.name] else [], - output: - join(macsN_dir, "{name}", "{name}_peaks.narrowPeak"), - params: - rname = 'MACS2_narrow', - gsize = config['references'][genome]['EFFECTIVEGENOMESIZE'], - macsver = config['tools']['MACSVER'], - flag = lambda w: "-c" if chip2input[w.name] else "", - shell: - """ - module load {params.macsver}; - macs2 callpeak \\ - -t {input.chip} {params.flag} {input.c_option} \\ - -g {params.gsize} \\ - -n {wildcards.name} \\ - --outdir {macsN_dir}/{wildcards.name} \\ - -q 0.01 \\ - --keep-dup="all" \\ - -f "BAMPE" """ \ No newline at end of file diff --git a/workflow/rules/single/peakcall.smk b/workflow/rules/single/peakcall.smk index 6f8825a..0657289 100644 --- a/workflow/rules/single/peakcall.smk +++ b/workflow/rules/single/peakcall.smk @@ -39,6 +39,35 @@ rule MACS2_broad: """ +rule MACS2_narrow: + input: + chip = join(bam_dir, "{name}.Q5DD.bam"), + txt = join(ppqt_dir, "{name}.Q5DD.ppqt.txt"), + c_option = lambda w: join(bam_dir, f"{chip2input[w.name]}.Q5DD_tagAlign.gz") + if w.name and chip2input[w.name] else [], + output: + join(macsN_dir, "{name}", "{name}_peaks.narrowPeak"), + params: + rname = 'MACS2_narrow', + gsize = config['references'][genome]['EFFECTIVEGENOMESIZE'], + macsver = config['tools']['MACSVER'], + flag = lambda w: "-c" if chip2input[w.name] else "", + shell: + """ + module load {params.macsver}; + ppqt_len=$(awk '{{print $1}}' {input.txt}) + macs2 callpeak \\ + -t {input.chip} {params.flag} {input.c_option} \\ + -g {params.gsize} \\ + -n {wildcards.name} \\ + --outdir {macsN_dir}/{wildcards.name} \\ + -q 0.01 \\ + --keep-dup="all" \\ + --nomodel \\ + --extsize $ppqt_len + """ + + rule SICER: input: chip = lambda w: join(bam_dir, w.name + ".Q5DD_tagAlign.gz"), @@ -111,30 +140,3 @@ rule SICER: """ -rule MACS2_narrow: - input: - chip = join(bam_dir, "{name}.Q5DD.bam"), - txt = join(ppqt_dir, "{name}.Q5DD.ppqt.txt"), - c_option = lambda w: join(bam_dir, f"{chip2input[w.name]}.Q5DD_tagAlign.gz") - if w.name and chip2input[w.name] else [], - output: - join(macsN_dir, "{name}", "{name}_peaks.narrowPeak"), - params: - rname = 'MACS2_narrow', - gsize = config['references'][genome]['EFFECTIVEGENOMESIZE'], - macsver = config['tools']['MACSVER'], - flag = lambda w: "-c" if chip2input[w.name] else "", - shell: - """ - module load {params.macsver}; - ppqt_len=$(awk '{{print $1}}' {input.txt}) - macs2 callpeak \\ - -t {input.chip} {params.flag} {input.c_option} \\ - -g {params.gsize} \\ - -n {wildcards.name} \\ - --outdir {macsN_dir}/{wildcards.name} \\ - -q 0.01 \\ - --keep-dup="all" \\ - --nomodel \\ - --extsize $ppqt_len - """ \ No newline at end of file diff --git a/workflow/rules/single/qc.smk b/workflow/rules/single/qc.smk index e0dc6fb..4e70175 100644 --- a/workflow/rules/single/qc.smk +++ b/workflow/rules/single/qc.smk @@ -109,7 +109,6 @@ rule kraken: --output {output.krakenout} \\ --gzip-compressed \\ {input.fq1} - fi # Generate Krona Report cut -f2,3 {output.krakenout} | \\ diff --git a/workflow/rules/single/trim_align_dedup.smk b/workflow/rules/single/trim_align_dedup.smk index e7d7005..cb9a1ca 100644 --- a/workflow/rules/single/trim_align_dedup.smk +++ b/workflow/rules/single/trim_align_dedup.smk @@ -171,7 +171,6 @@ rule dedup: out_bam = join(bam_dir, "{name}.Q5DD.bam"), flagstat = join(bam_dir, "{name}.Q5DD.bam.flagstat"), idxstat = join(bam_dir, "{name}.Q5DD.bam.idxstat"), - bwadups = join(bam_dir, "{name}.bwa.Q5.duplic"), tagalign = join(bam_dir, "{name}.Q5DD_tagAlign.gz") params: rname = 'dedup', @@ -281,7 +280,7 @@ rule bam2bw: """ input: bam = lambda w: join(bam_dir, f"{w.name}.{w.ext}.bam"), - ppqt = lambda w: join(ppqt_dir, f"{w.name}.{w.ext}.ppqt.txt") + ppqt = lambda w: join(ppqt_dir, f"{w.name}.{w.ext}.ppqt.txt"), output: outbw = join(bw_dir, "{name}.{ext}.RPGC.bw"), params: @@ -299,8 +298,7 @@ rule bam2bw: tmp=$(mktemp -d -p "{params.tmpdir}") trap 'rm -rf "${{tmp}}"' EXIT - bam_cov_option=--centerReads - echo "printing out value of bam-cov-option $bam_cov_option" + ppqt_len=$(awk '{{print $1}}' {input.ppqt}) bamCoverage \\ --bam {input.bam} \\ @@ -310,5 +308,5 @@ rule bam2bw: --numberOfProcessors {threads} \\ --normalizeUsing RPGC \\ --effectiveGenomeSize {params.effectivegenomesize} \\ - ${{bam_cov_option}}; + -e ${{ppqt_len}} """ \ No newline at end of file