METHYLDACKEL EXTRACT

Calculate per-base CpG metrics with MethylDackel

URL: https://github.com/dpryan79/MethylDackel

Example

This wrapper can be used in the following way:

rule test_methyldackel_extract_fraction:
    input:
        ref="chgchh.fa",
        aln="chgchh_aln.bam",
        aln_bai="chgchh_aln.bam.bai",
        bw="chgchh.bw",
    output:
        cpg="cpg.meth.bg",
        chh="chh.meth.bg",
        chg="chg.meth.bg",
    log:
        "logs/md_extract/meth.log",
    params:
        extra="--fraction",
    wrapper:
        "v3.9.0-14-g476823b/bio/methyldackel/extract"


rule test_methyldackel_extract_count:
    input:
        ref="chgchh.fa",
        aln="chgchh_aln.bam",
        aln_bai="chgchh_aln.bam.bai",
        bw="chgchh.bw",
    output:
        cpg="cpg.count.bg",
        chh="chh.count.bg",
        chg="chg.count.bg",
    log:
        "logs/md_extract/count.log",
    params:
        extra="--counts",
    wrapper:
        "v3.9.0-14-g476823b/bio/methyldackel/extract"


rule test_methyldackel_extract_logit:
    input:
        ref="chgchh.fa",
        aln="chgchh_aln.bam",
        aln_bai="chgchh_aln.bam.bai",
        bw="chgchh.bw",
    output:
        cpg="cpg.logit.bg",
        chh="chh.logit.bg",
        chg="chg.logit.bg",
    log:
        "logs/md_extract/logit.log",
    params:
        extra="--logit",
    wrapper:
        "v3.9.0-14-g476823b/bio/methyldackel/extract"


rule test_methyldackel_extract_report:
    input:
        ref="chgchh.fa",
        aln="chgchh_aln.bam",
        aln_bai="chgchh_aln.bam.bai",
        bw="chgchh.bw",
    output:
        cytosine_report="report.tsv",
    log:
        "logs/md_extract/report.log",
    params:
        extra="",
    wrapper:
        "v3.9.0-14-g476823b/bio/methyldackel/extract"

Note that input, output and log file paths can be chosen freely.

When running with

snakemake --use-conda

the software dependencies will be automatically deployed into an isolated environment before execution.

Software dependencies

• methyldackel=0.6.1

Input/Output

Input:

• ref: Path to reference genome (Fasta)

• aln: Path to aligned reads (BAM)

Output:

• cpg: Path to CpG metrics (BedGraph)

• chg: Optional path to CHG metrics (BedGraph)

• chh: Optional path to CHH metrics (BedGraph)

Params

• extra: Optional parameters provided to MethylDackel extract

Authors

• Thibault Dayris

Code

#!/bin/env python3
# coding: utf-8

from snakemake import shell
from tempfile import TemporaryDirectory


log = snakemake.log_fmt_shell(stdout=True, stderr=True, append=True)
extra = snakemake.params.get("extra", "")
extension = "bedGraph"
if "fraction" in extra:
    extension = "meth.bedGraph"
elif "counts" in extra:
    extension = "counts.bedGraph"
elif "logit" in extra:
    extension = "logit.bedGraph"

# Input files
bed = snakemake.input.get("bed", "")
if bed:
    extra += f" -l {bed} "


bw = snakemake.input.get("bw", "")
if bw:
    extra += f" --mappability {bw} "


bbm_in = snakemake.input.get("bbm", "")
if bbm_in:
    extra += f" --mappabilityBB {bbm_in}"


# Output files
chg = snakemake.output.get("chg")
if chg:
    extra += " --CHG "


chh = snakemake.output.get("chh")
if chh:
    extra += " --CHH "


methylkit = snakemake.output.get("methylkit")
if methylkit:
    extra += " --methylKit "
    extension = "methylKit"


report = snakemake.output.get("cytosine_report")
if report:
    extra += " --cytosine_report "


if not (snakemake.output.get("cpg") or report):
    extra += " --noCpG "

with TemporaryDirectory() as tempdir:
    shell(
        "MethylDackel extract "
        "{extra} "
        "-@ {snakemake.threads} "
        "--opref {tempdir}/methyldackel_results "
        "{snakemake.input.ref} "
        "{snakemake.input.aln} "
        "{log} "
    )

    if report:
        shell(
            "mv --verbose {tempdir}/methyldackel_results.cytosine_report.txt {report} {log}"
        )
    else:
        if snakemake.output.get("cpg"):
            shell(
                "mv --verbose {tempdir}/methyldackel_results_CpG.{extension} {snakemake.output.cpg} {log}"
            )

        if chg:
            shell(
                "mv --verbose {tempdir}/methyldackel_results_CHG.{extension} {chg} {log}"
            )

        if chh:
            shell(
                "mv --verbose {tempdir}/methyldackel_results_CHH.{extension} {chh} {log}"
            )