New upstream version 2.2.0+dfsg

## mapDamage

[![install with bioconda](https://img.shields.io/badge/install%20with-bioconda-brightgreen.svg?style=flat)](http://bioconda.github.io/recipes/mapdamage2/README.html) [![Conda](https://img.shields.io/conda/dn/bioconda/mapdamage2.svg)](https://anaconda.org/bioconda/mapdamage2/files)

[![install with bioconda](https://img.shields.io/badge/install%20with-bioconda-brightgreen.svg?style=flat)](http://bioconda.github.io/recipes/mapdamage2/README.html) [![Conda](https://img.shields.io/conda/dn/bioconda/mapdamage2.svg)](https://anaconda.org/bioconda/mapdamage2/files) ![Conda](https://anaconda.org/bioconda/mapdamage2/badges/latest_release_date.svg) ![Conda](https://anaconda.org/bioconda/mapdamage2/badges/version.svg) [![Project Status: Inactive – The project has reached a stable, usable state but is no longer being actively developed; support/maintenance will be provided as time allows.](https://www.repostatus.org/badges/latest/inactive.svg)](https://www.repostatus.org/#inactive)

#### `bioconda` installation

- python3 version **2.2.0**

```

conda install -c bioconda mapdamage2=2.2.0

```

- python3 version **2.2.0** **with** R and 4 mandatory packages for the Bayesian inference:

```

conda install -c bioconda mapdamage2=2.2.0=pyr36_1

```

---

### Important

Users with versions dating prior to June 12 2013 please update. A nasty bug that caused the statistical part of `mapDamage` to use half of the data for estimation of the damage parameters, sorry for the inconvenience.

- From version `2.2.0` the `master` branch is requiring **python3** as `python2` is not supported from 2020-01-01.

- Users with versions dating prior to June 12 2013 please update. A nasty bug that caused the statistical part of `mapDamage` to use half of the data for estimation of the damage parameters, sorry for the inconvenience.

### Introduction

Complete documentation, instructions, examples, screenshots and FAQ are available at [this address](http://ginolhac.github.io/mapDamage/).

[mapDamage2](http://geogenetics.ku.dk/publications/mapdamage2.0/) is a computational framework written in **Python3** and **R**, which tracks and quantifies DNA damage patterns

[mapDamage2](https://geogenetics.ku.dk/publications/mapdamage2.0/) is a computational framework written in **Python3** and **R**, which tracks and quantifies DNA damage patterns

among ancient DNA sequencing reads generated by Next-Generation Sequencing platforms.

`mapDamage` was developed at the [Centre for GeoGenetics](http://geogenetics.ku.dk/) by the [Orlando Group ](http://geogenetics.ku.dk/research/research_groups/palaeomix_group/).

`mapDamage` was developed at the [Centre for GeoGenetics](https://geogenetics.ku.dk/) by the [Orlando Group ](https://geogenetics.ku.dk/research_groups/palaeomix_group/).

### Citation

http://bioinformatics.oxfordjournals.org/content/27/15/2153](http://bioinformatics.oxfordjournals.org/content/27/15/2153)

### Test

### Test the no-stats part and rescaling

in the package, you can test `mapDamage` by running:

you can test `mapDamage` by running:

```

cd mapDamage/mapdamage/

python3 mp_test.py

```

should return

```

Started with the command: /usr/local/bin/mapDamage -i tests/test.bam -r tests/fake1.fasta -d tests/results --no-stats

	Reading from 'tests/test.bam'

	Writing results to 'tests/results/'

pdf tests/results/Fragmisincorporation_plot.pdf generated

additional tests/results/Length_plot.pdf generated

Successful run

.

----------------------------------------------------------------------

Ran 2 tests in 3.357s

OK

```

### Contact

Please report bugs and suggest possible improvements to Aurélien Ginolhac, Mikkel Schubert or Hákon Jónsson by email:

import unittest

import subprocess

import pysam

import mapdamage

import optparse

import filecmp

def read_nocomment(dnacomp):

    with open(dnacomp, 'r') as f:

        a = f.readlines()

        return [x for x in a if not x.startswith('#')]

def mock_options(filename,rescale_out,folder):

    """Make the options object with nice values for testing"""

    return optparse.Values({

        "filename":filename,

        "rescale_out":rescale_out,

        "verbose":True,

        "folder":folder,

        "rescale_length_5p":12, # default values as in --seq-length

        "rescale_length_3p":12, # default values as in --seq-length

        "quiet":True

        })

class testCases(unittest.TestCase):

    def test_no_stats(self):

        subprocess.run(["mapDamage",  "-i",  "tests/test.bam", "-r", "tests/fake1.fasta", "-d", "tests/results", "--no-stats"], check = True)

        self.assertTrue(read_nocomment("tests/dnacomp.txt") == read_nocomment("tests/results/dnacomp.txt"))

class testRescaling(unittest.TestCase):

    def test_single_end_file(self):

        """Test, rescaling BAM file"""

        #

	    # The expected substition frequencies before and after scaling using the scaled qualities as probalities:

	    # CT	0.06226411977920493		0.04163524443356556

	    # TC	0.020395286584806528		0.020395286584806528

	    # GA	0.04400459948304954		0.03794905109091021

	    # AG	0.05355350777642983		0.05355350777642983

	    # Quality metrics before and after scaling

	    # CT-Q0 	5		5

	    # CT-Q10 	5		2

	    # CT-Q20 	3		2

	    # CT-Q30 	3		2

	    # CT-Q40 	0		0

	    # GA-Q0 	5		5

	    # GA-Q10 	5		4

	    # GA-Q20 	1		0

	    # GA-Q30 	1		0

	    # GA-Q40 	0		0

        options = mock_options("tests/test.bam","tests/test.rescaled.sam","tests/probs/")

        ref = pysam.Fastafile("tests/fake1.fasta")

        mapdamage.rescale.rescale_qual(ref,options,debug=True)

        self.assertTrue(filecmp.cmp("tests/test.rescaled.sam","tests/test.rescaled.correct.sam"))

if  __name__=='__main__':

    unittest.main()

    if not os.path.isfile(full_path):

        sys.exit("Missing file, the file \n\tStats_out_MCMC_correct_prob.csv\nshould be in the folder\n\t"+folder+"\nDid you run the MCMC estimates of the parameters?")

    try:

        fi_handle = csv.DictReader(open(full_path))

        with open(full_path) as fi:

            fi_handle = csv.DictReader(fi)

            corr_prob = {}

            for line in fi_handle:

                if (line["Position"] in corr_prob):

            for key in list(corr_prob.keys()):

                if key < -rescale_length_3p or key > rescale_length_5p:

                    corr_prob.pop(key)

            return corr_prob

    except csv.Error as e:

        sys.exit('File %s, line %d: %s' % (os.path.join(folder,"Stats_out_MCMC_correct_prob.csv"), \

"","Position","C.T","G.A"

"1",1,0.653616777018436,0

"2",2,0.574821404505089,0

"3",3,0.524660717195785,0

"4",4,0.491676371805506,0

"5",5,0.46919312975022,0

"6",6,0.453361082263795,0

"7",7,0.441882878430246,0

"8",8,0.433345603590333,0

"9",9,0.426853067214069,0

"10",10,0.4218187429273,0

"11",11,0.417847465539491,0

"12",12,0.414666179579865,0

"13",-12,0,0.427459648772127

"14",-11,0,0.430678304009036

"15",-10,0,0.434694598982698

"16",-9,0,0.439783285536968

"17",-8,0,0.446341431282045

"18",-7,0,0.454957458647015

"19",-6,0,0.466528745659899

"20",-5,0,0.482466321045583

"21",-4,0,0.505055143009322

"22",-3,0,0.538101706832429

"23",-2,0,0.588153916474448

"24",-1,0,0.666277684522864

@SQ	SN:fake1	LN:201

@PG	ID:bwa	PN:bwa	VN:0.7.17-r1188	CL:bwa samse -f input.sam fake1.fasta - input.fastq

fake1:50-100	0	fake1	69	37	31M	*	0	0	CCATGTCGGGCAGGCTGGTCTCGAACTCCTG	////////E6/EE/E/A////E//AAAE/EE	XT:A:U	NM:i:0	X0:i:1	X1:i:0	XM:i:0	XO:i:0	XG:i:0	MD:Z:31	MR:f:0

fake1:50-100b	0	fake1	69	37	31M	*	0	0	CTATGTCGGGCAGGCTGGTCTCGAACTCCTG	/#//////E6/EE/E/A////E//AAAE/EE	XT:A:U	NM:i:1	X0:i:1	X1:i:0	XM:i:1	XO:i:0	XG:i:0	MD:Z:1C29	MR:f:0.57482

fake1:50-100c	4	*	0	0	*	*	0	0	CGGTAGAGATGGAGTTTCACCATGTCGGGCAAG	//////////////A////////////E6/EE/

fake1:50-100d	4	*	0	0	*	*	0	0	TAATAGAGATGGAGTTTCACCATGTCGTGCAAG	//////////////A////////////E6/EE/

fake1:70-120	0	fake1	70	37	38M	*	0	0	TATGTCGGGCAGGCTGGTCTCGAACTCCTGACCTCAGA	#////6EE//AA6E//E66E/EAEEEEEAEEEAAEEA#	XT:A:U	NM:i:2	X0:i:1	X1:i:0	XM:i:2	XO:i:0	XG:i:0	MD:Z:0C36G0	MR:f:1.31989

fake1:80-130	16	fake1	80	25	41M	*	0	0	GAGCTGGTCTCGAACTCCTGACCTCAGGCGATCTGCCTGTC	AAABBCCDDFEEEEEAEEEEA/EE/E/AAE//A//E/A///	XT:A:U	NM:i:3	X0:i:1	X1:i:0	XM:i:3	XO:i:0	XG:i:0	MD:Z:0A0G37C1	MR:f:0

fake1:80-130b	16	fake1	80	25	51M	*	0	0	AGGCCGGTCTCGAACTCCTGACCTCAGGCGATCTGCCTGCCTTAACCTCCC	AAABBCCDDFEEEEEAEEEEA/EE/E/AAE//A//E/A///E/A$/EE66/	XT:A:U	NM:i:3	X0:i:1	X1:i:0	XM:i:3	XO:i:0	XG:i:0	MD:Z:4T37C1G6	MR:f:0.44188

fake1:80-130c	16	fake1	80	37	51M	*	0	0	AGGCTGGTCTCGAACTCCTGACCTCAGGCGATCTGCCTGCCTTAGCCCCCC	AAABBCCDDFEEEEEAEEEEA/EE/E/AAE//A//E/A///E/A//EE66/	XT:A:U	NM:i:2	X0:i:1	X1:i:0	XM:i:2	XO:i:0	XG:i:0	MD:Z:42C4T3	MR:f:0

fake1:80-130d	16	fake1	80	25	51M	*	0	0	AGGCTGGTCTCGAACTCCTGACCTCAGACGATCTGCCTGCCTTAGCCCCCC	AAABBCCDDFEEEEEAEEEEA/EE/E/AAE//A//E/A///E/A//EE66/	XT:A:U	NM:i:3	X0:i:1	X1:i:0	XM:i:3	XO:i:0	XG:i:0	MD:Z:27G14C4T3	MR:f:0