snakemake_macs.py

conf = ... # load the config file here


########### some helper functions to aggregate input files ##########


def prefix2samples(wildcards):
    # retreive the treatment or control samples for the given prefix here
    return expand("{sample}.bam", sample=...)


def prefix2merged(wildcards):
    types = "treatment control" if conf.control_targets else "treatment"
    return expand("{prefix}_{type}.bam", prefix=wildcards.prefix, type=wildcards.types)


def type2bdg:
    return expand("{{prefix}}_{type}.bdg.tmp", type="control_lambda" if wildcards.type == "control" else "treat_pileup")



################# Snakemake rules ###############


# collect results
rule all:
    input: expand("{prefix}_macs2.pdf", prefix=...) # collect all the results for all prefixes (retreived from the conf)


# merge bam files
rule merge:
    input:
        prefix2samples
    output:
        "{prefix}_{type}.bam"
    params:
        out="{prefix}_{type}"
    run:
        if len(input) == 1:
            shell("cp {input} {output}")
        else:
            shell("samtools merge {input} {params.out}")


# run macs
rule macs2:
    input:
        prefix2merged
    output:
        "{prefix}_peaks.bed",
        "{prefix}_control_lambda.bdg", ...
    params:
        keep_dup=1,
        shiftsize=73,
        description={prefix}
    run:
        input = " -c ".join(input)
        shell(
            "macs2 callpeak -B -q 0.01 --nomodel"
            "--keep-dup {params.keep_dup} "
            "--shiftsize={params.shiftsize} "
            "-n {params.description} "
            "-t {input}")


rule bdg_trim_control:
    input:
        bdg="{prefix}_control_lambda.bdg",
        chrom_bed=conf.get_path("lib", "chrom_bed")
    output:
        "{prefix}_control_lambda.bdg.tmp"
    shell:
        "bedtools intersect -a {input.bdg} -b {input.chrom_bed} -wa -f 1.00 > {output}"


rule bdg2bw:
    input:
        type2bdg,
        conf.get("lib", "chrom_len")
    output:
        "{prefix}_{type}.bw"
    shell:
        "bedGraphToBigWig {input} {output}"


rule stats:
    input:
        macs2_peaks_xls="{prefix}_peaks.xls",
        db=ChiLinQC_db,
        template=rlang_template
    output:
        json=expand("{json_prefix}_macs2.json", json_prefix=conf.json_prefix),
        R="{prefix}_macs2.R",
        pdf="{prefix}_macs2.pdf"
    run:
        stat_macs2(id=conf.id, **input, **output)