MAPDAMAGE2

tracking and quantifying damage patterns in ancient DNA sequences. For more information about MapDamage2 see MapDamage2 documentation.

Example

This wrapper can be used in the following way:

rule mapdamage2:
    input:
        ref="genome.fasta",
        bam="mapped/{sample}.bam",
    output:
        log="results/{sample}/Runtime_log.txt",  # output folder is infered from this file, so it needs to be the same folder for all output files
        GtoA3p="results/{sample}/3pGtoA_freq.txt",
        CtoT5p="results/{sample}/5pCtoT_freq.txt",
        dnacomp="results/{sample}/dnacomp.txt",
        frag_misincorp="results/{sample}/Fragmisincorporation_plot.pdf",
        len="results/{sample}/Length_plot.pdf",
        lg_dist="results/{sample}/lgdistribution.txt",
        misincorp="results/{sample}/misincorporation.txt",
#        rescaled_bam="results/{sample}.rescaled.bam", # uncomment if you want the rescaled BAM file
    params:
        extra="--no-stats"  # optional parameters for mapdamage2 (except -i, -r, -d, --rescale)
    log:
        "logs/{sample}/mapdamage2.log"
    threads: 1  # MapDamage2 is not threaded
    wrapper:
        "0.75.0-13-g0997adf/bio/mapdamage2"

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

  • mapdamage2=2.2

Authors

  • Filipe G. Vieira

Code

__author__ = "Filipe G. Vieira"
__copyright__ = "Copyright 2020, Filipe G. Vieira"
__license__ = "MIT"

import os.path
from snakemake.shell import shell

extra = snakemake.params.get("extra", "")
log = snakemake.log_fmt_shell(stdout=True, stderr=True)

in_bam = snakemake.input.get("bam", "")
if in_bam:
    in_bam = "--input " + in_bam

output_folder = os.path.dirname(snakemake.output.get("log", ""))
if not output_folder:
    raise ValueError("mapDamage2 rule needs output 'log'.")

rescaled_bam = snakemake.output.get("rescaled_bam", "")
if rescaled_bam:
    rescaled_bam = "--rescale-out " + rescaled_bam


shell(
    "mapDamage "
    "{in_bam} "
    "--reference {snakemake.input.ref} "
    "--folder {output_folder} "
    "{rescaled_bam} "
    "{extra} "
    "{log}"
)