Skip to content
Snippets Groups Projects
Commit 57b3ece0 authored by Laura Cook's avatar Laura Cook
Browse files

pool peaks using bedtools merge instead

parent 28bc62ed
Branches
No related tags found
No related merge requests found
......@@ -421,57 +421,57 @@ rule frip:
# Based on ENCODE `overlap_peaks.py` - recommended for histone marks.
# Need to pool peaks first for each replicate
# rule pool_peaks:
# input:
# E10 = expand(["results/macs2/{case}_vs_{control}_E10.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
# E11 = expand(["results/macs2/{case}_vs_{control}_E11.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
# E12 = expand(["results/macs2/{case}_vs_{control}_E12.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
# E13 = expand(["results/macs2/{case}_vs_{control}_E13.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
# E14 = expand(["results/macs2/{case}_vs_{control}_E14.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
# E15 = expand(["results/macs2/{case}_vs_{control}_E15.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK)
# output:
# E10 = expand("results/macs2/E10.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
# E11 = expand("results/macs2/E11.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
# E12 = expand("results/macs2/E12.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
# E13 = expand("results/macs2/E13.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
# E14 = expand("results/macs2/E14.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
# E15 = expand("results/macs2/E15.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK)
# run:
# shell("cat {input.E10} > {output.E10}")
# shell("cat {input.E11} > {output.E11}")
# shell("cat {input.E12} > {output.E12}")
# shell("cat {input.E13} > {output.E13}")
# shell("cat {input.E14} > {output.E14}")
# shell("cat {input.E15} > {output.E15}")
rule pool_peaks:
input:
E10 = expand(["results/macs2/{case}_vs_{control}_E10.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
E11 = expand(["results/macs2/{case}_vs_{control}_E11.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
E12 = expand(["results/macs2/{case}_vs_{control}_E12.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
E13 = expand(["results/macs2/{case}_vs_{control}_E13.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
E14 = expand(["results/macs2/{case}_vs_{control}_E14.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
E15 = expand(["results/macs2/{case}_vs_{control}_E15.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK)
output:
E10 = expand("results/macs2/E10.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
E11 = expand("results/macs2/E11.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
E12 = expand("results/macs2/E12.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
E13 = expand("results/macs2/E13.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
E14 = expand("results/macs2/E14.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
E15 = expand("results/macs2/E15.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK)
run:
shell("cat {input.E10} | sort -k1,1 -k2,2n - | bedtools merge -i - > {output.E10}")
shell("cat {input.E11}| sort -k1,1 -k2,2n - | bedtools merge -i - > {output.E11}")
shell("cat {input.E12} | sort -k1,1 -k2,2n - | bedtools merge -i - > {output.E12}")
shell("cat {input.E13} | sort -k1,1 -k2,2n - | bedtools merge -i - > {output.E13}")
shell("cat {input.E14} | sort -k1,1 -k2,2n - | bedtools merge -i - > {output.E14}")
shell("cat {input.E15} | sort -k1,1 -k2,2n - | bedtools merge -i - > {output.E15}")
# rule overlap_peaks:
# input:
# E10= expand(["results/macs2/{case}_vs_{control}_E10.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
# E11= expand(["results/macs2/{case}_vs_{control}_E11.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
# E12= expand(["results/macs2/{case}_vs_{control}_E12.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
# E13= expand(["results/macs2/{case}_vs_{control}_E13.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
# E14= expand(["results/macs2/{case}_vs_{control}_E14.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
# E15= expand(["results/macs2/{case}_vs_{control}_E15.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
# E10_pooled=expand("results/macs2/E10.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
# E11_pooled=expand("results/macs2/E11.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
# E12_pooled=expand("results/macs2/E12.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
# E13_pooled=expand("results/macs2/E13.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
# E14_pooled=expand("results/macs2/E14.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
# E15_pooled=expand("results/macs2/E15.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
# output:
# E10= expand("results/macs2/{mark}_E10.5_overlap.narrowPeak", mark=MARK),
# E11= expand("results/macs2/{mark}_E11.5_overlap.narrowPeak", mark=MARK),
# E12= expand("results/macs2/{mark}_E12.5_overlap.narrowPeak", mark=MARK),
# E13= expand("results/macs2/{mark}_E13.5_overlap.narrowPeak", mark=MARK),
# E14= expand("results/macs2/{mark}_E14.5_overlap.narrowPeak", mark=MARK),
# E15= expand("results/macs2/{mark}_E15.5_overlap.narrowPeak", mark=MARK)
# run:
# shell("python2.7 scripts/overlap_peaks.py {input.E10} {input.E10_pooled} {output.E10}")
# shell("python2.7 scripts/overlap_peaks.py {input.E11} {input.E11_pooled} {output.E11}")
# shell("python2.7 scripts/overlap_peaks.py {input.E12} {input.E12_pooled} {output.E12}")
# shell("python2.7 scripts/overlap_peaks.py {input.E13} {input.E13_pooled} {output.E13}")
# shell("python2.7 scripts/overlap_peaks.py {input.E14} {input.E14_pooled} {output.E14}")
# shell("python2.7 scripts/overlap_peaks.py {input.E15} {input.E15_pooled} {output.E15}")
rule overlap_peaks:
input:
E10= expand(["results/macs2/{case}_vs_{control}_E10.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
E11= expand(["results/macs2/{case}_vs_{control}_E11.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
E12= expand(["results/macs2/{case}_vs_{control}_E12.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
E13= expand(["results/macs2/{case}_vs_{control}_E13.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
E14= expand(["results/macs2/{case}_vs_{control}_E14.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
E15= expand(["results/macs2/{case}_vs_{control}_E15.5_{mark}_macs2_peaks.narrowPeak"], zip, case=IPS, control=INPUTS, mark=MARK),
E10_pooled=expand("results/macs2/E10.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
E11_pooled=expand("results/macs2/E11.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
E12_pooled=expand("results/macs2/E12.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
E13_pooled=expand("results/macs2/E13.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
E14_pooled=expand("results/macs2/E14.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
E15_pooled=expand("results/macs2/E15.5_{mark}_macs2_pooled_peaks.narrowPeak", mark=MARK),
output:
E10= expand("results/macs2/{mark}_E10.5_overlap.narrowPeak", mark=MARK),
E11= expand("results/macs2/{mark}_E11.5_overlap.narrowPeak", mark=MARK),
E12= expand("results/macs2/{mark}_E12.5_overlap.narrowPeak", mark=MARK),
E13= expand("results/macs2/{mark}_E13.5_overlap.narrowPeak", mark=MARK),
E14= expand("results/macs2/{mark}_E14.5_overlap.narrowPeak", mark=MARK),
E15= expand("results/macs2/{mark}_E15.5_overlap.narrowPeak", mark=MARK)
run:
shell("python2.7 scripts/overlap_peaks.py {input.E10} {input.E10_pooled} {output.E10}")
shell("python2.7 scripts/overlap_peaks.py {input.E11} {input.E11_pooled} {output.E11}")
shell("python2.7 scripts/overlap_peaks.py {input.E12} {input.E12_pooled} {output.E12}")
shell("python2.7 scripts/overlap_peaks.py {input.E13} {input.E13_pooled} {output.E13}")
shell("python2.7 scripts/overlap_peaks.py {input.E14} {input.E14_pooled} {output.E14}")
shell("python2.7 scripts/overlap_peaks.py {input.E15} {input.E15_pooled} {output.E15}")
# ===============================================================================================
# 10. Combine all QC into multiqc output
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment