diff --git a/mouse/Snakefile b/mouse/Snakefile index 86fd5c51a10487cb4d18c5e242e81cbd079e5bea..bbb72d022c4350795265277a202ccdbcbb90fd33 100644 --- a/mouse/Snakefile +++ b/mouse/Snakefile @@ -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