.. _`bio/methyldackel/extract`:

METHYLDACKEL EXTRACT
====================


.. image:: https://img.shields.io/github/issues-pr/snakemake/snakemake-wrappers/bio/methyldackel/extract?label=version%20update%20pull%20requests
   :target: https://github.com/snakemake/snakemake-wrappers/pulls?q=is%3Apr+is%3Aopen+label%3Abio/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:

.. code-block:: python

    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.0.1/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.0.1/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.0.1/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.0.1/bio/methyldackel/extract"

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

When running with

.. code-block:: bash

    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
----

.. code-block:: python

    #!/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}"
                )


.. |nl| raw:: html

   <br>