README.md

+## 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)
+
+---
+
 ### 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.
 
 ### 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 Python and R, which tracks and quantifies DNA damage patterns
+[mapDamage2](http://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` is developed at the [Centre for GeoGenetics](http://geogenetics.ku.dk/) by the [Orlando Group ](http://geogenetics.ku.dk/research/research_groups/palaeomix_group/).
-
-### Web interface
-
-Anna Kostikova from [insideDNA](https://insidedna.me) implemented a web interface for mapDamage.  
-Users can adjust the cloud resources in terms of RAM/CPU to perform their analysis. She wrote a [tutorial](https://insidedna.me/tutorials/view/Analysis-ancient-DNA-samples-using-mapDamage) explaining how to use this [web interface (sign up required)](https://insidedna.me/app#/tools/100648/)
+`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/).
 
 
 ### Citation
@@ -22,12 +23,21 @@ Jónsson H, Ginolhac A, Schubert M, Johnson P, Orlando L.
 _Bioinformatics_ 23rd April 2013. doi: 10.1093/bioinformatics/btt193](http://bioinformatics.oxfordjournals.org/content/early/2013/05/17/bioinformatics.btt193)
 
 
-The original mapDamage1 is described in the following article:  
+The original `mapDamage1` was published in the following article:  
 Ginolhac A, Rasmussen M, Gilbert MT, Willerslev E, Orlando L.
 [mapDamage: testing for damage patterns in ancient DNA sequences. _Bioinformatics_ 2011 **27**(15):2153-5
 http://bioinformatics.oxfordjournals.org/content/27/15/2153](http://bioinformatics.oxfordjournals.org/content/27/15/2153)
 
 
+### Test
+
+in the package, you can test `mapDamage` by running:
+
+```
+cd mapDamage/mapdamage/
+python3 mp_test.py
+```
+
 
 ### Contact
 Please report bugs and suggest possible improvements to Aurélien Ginolhac, Mikkel Schubert or Hákon Jónsson by email:

bin/mapDamage

@@ -139,6 +139,20 @@ def _check_mm_option():
             raise SystemExit("The --mapdamage-modules path (path=%s) must contain the mapdamage module" % path_to_mm)
     return path_to_mm
 
+def is_paired(alignfile):
+    """
+    check if an alignment file contains paired end data by
+    looking at the first 100 reads
+    """
+
+    n = 0
+    for read in alignfile:
+        if read.is_paired:
+            return True
+        n += 1
+        if n == 100:
+            return False
+
 
 def main():
     start_time = time.time()
@@ -227,6 +241,13 @@ def main():
     else:
         in_bam = pysam.Samfile(options.filename)
 
+    # exit SAM/BAM is paired (look at top 100 reads)
+    if is_paired(in_bam):
+        logger.error("alignment file must be single-ended")
+        sys.stderr.write("alignment file must be single-ended\n")
+        sys.exit(1)
+    # reset SAM/BAM to the begining
+    in_bam = pysam.Samfile(options.filename)
 
     reflengths = dict(list(zip(in_bam.references, in_bam.lengths)))
     # check if references in SAM/BAM are the same in the fasta reference file

debian/changelog

+mapdamage (2.1.1+dfsg-1) unstable; urgency=medium
+
+  * Team upload.
+  * New upstream version 
+    - eliminates two patches
+      + na.rm 
+      + import
+    - rescale_test.py no longer available - simpler build process
+  * Standards-Version: 4.4.1
+
+ -- Steffen Moeller <moeller@debian.org>  Sat, 30 Nov 2019 18:52:19 +0100
+
 mapdamage (2.1.0+dfsg-1) unstable; urgency=medium
 
   * New upstream version

debian/control

@@ -8,7 +8,7 @@ Build-Depends: debhelper-compat (= 12),
                dh-python,
                python3-setuptools,
                python3-pysam <!nocheck>
-Standards-Version: 4.4.0
+Standards-Version: 4.4.1
 Vcs-Browser: https://salsa.debian.org/med-team/mapdamage
 Vcs-Git: https://salsa.debian.org/med-team/mapdamage.git
 Homepage: https://ginolhac.github.io/mapDamage/

debian/install

-mapdamage/rescale_test	usr/share/doc/mapdamage/tests
+#No longer available?
+#mapdamage/rescale_test	usr/share/doc/mapdamage/tests

debian/patches/fix_quantile_na_rm.patch

-Author: Nadiya Sitdykova
-Description: change calcSampleStats such as any NA and Nan's are 
- removed before quantiles are computed
-Last-Update: 2017-03-28
-Bug: https://github.com/ginolhac/mapDamage/issues/18
-
---- a/mapdamage/Rscripts/stats/function.R
-+++ b/mapdamage/Rscripts/stats/function.R
-@@ -595,8 +595,8 @@ calcSampleStats <- function(da,X){
-                       pos=da[,"Pos"],
-                       mea=apply(X,1,mean),
-                       med=apply(X,1,median),
--                      loCI=apply(X,1,quantile,c(0.025)),
--                      hiCI=apply(X,1,quantile,c(0.975))
-+                      loCI=apply(X,1,quantile,c(0.025),na.rm=TRUE),
-+                      hiCI=apply(X,1,quantile,c(0.975),na.rm=TRUE)
-                       ))
- }
- 

debian/patches/fix_rescale_test.patch

-Author: Andreas Tille <tille@debian.org>
-Last-Update: Tue, 01 Oct 2019 17:13:12 +0200
-Description: Fix import statement in rescale_test.py
-
---- a/mapdamage/rescale_test.py
-+++ b/mapdamage/rescale_test.py
-@@ -1,6 +1,6 @@
- import unittest
- import optparse
--from . import rescale
-+from mapdamage import rescale
- import pysam
- import filecmp
- 

debian/patches/series

 use_debian_packaged_seqtk.patch
-fix_quantile_na_rm.patch
-fix_rescale_test.patch

debian/patches/use_debian_packaged_seqtk.patch

@@ -3,8 +3,10 @@ Last-Update: Thu, 28 Jul 2016 15:13:14 +0200
 Bug-Debian: https://bugs.debian.org/859090
 Description: Use Debian packaged seqtk
 
---- a/setup.py
-+++ b/setup.py
+Index: mapdamage/setup.py
+===================================================================
+--- mapdamage.orig/setup.py
++++ mapdamage/setup.py
 @@ -6,21 +6,6 @@ from distutils.command.install import in
  import os
  import subprocess
@@ -54,9 +56,11 @@ Description: Use Debian packaged seqtk
      scripts=['bin/mapDamage'],
      url='https://github.com/ginolhac/mapDamage',
      license='LICENSE.txt',
---- a/bin/mapDamage
-+++ b/bin/mapDamage
-@@ -149,7 +149,7 @@ def main():
+Index: mapdamage/bin/mapDamage
+===================================================================
+--- mapdamage.orig/bin/mapDamage
++++ mapdamage/bin/mapDamage
+@@ -163,7 +163,7 @@ def main():
          sys.path.insert(0,path_to_mm)
      import mapdamage
  
@@ -65,8 +69,10 @@ Description: Use Debian packaged seqtk
      if not (os.path.isfile(fpath_seqtk) and os.access(fpath_seqtk, os.X_OK)):
          sys.stderr.write("Seqtk executable not accessible; mapDamage has not\n"
                           "been intalled properly or current user does not\n"
---- a/mapdamage/composition.py
-+++ b/mapdamage/composition.py
+Index: mapdamage/mapdamage/composition.py
+===================================================================
+--- mapdamage.orig/mapdamage/composition.py
++++ mapdamage/mapdamage/composition.py
 @@ -33,7 +33,7 @@ def get_base_comp(filename,destination=F
      Gets the basecomposition of all the sequences in filename
      and returns the value to destination if given.

debian/rules

@@ -13,10 +13,15 @@ override_dh_fixperms:
 	find debian -name checkLibraries.R -exec chmod +x \{\} \;
 	find debian -name runGeneral.R     -exec chmod +x \{\} \;
 
-override_dh_install:
-	dh_install
-	mv `find debian -name rescale_test.py` $(TESTDIR)
-	sed -i 's/^import rescale/from mapdamage &/' $(TESTDIR)/rescale_test.py
+# rescale_test.py seems to have left the building
+#override_dh_install:
+#	dh_install
+#	mv `find debian -name rescale_test.py` $(TESTDIR)
+#	sed -i 's/^import rescale/from mapdamage &/' $(TESTDIR)/rescale_test.py
 
 override_dh_python3:
 	dh_python3 --no-ext-rename
+
+override_dh_auto_clean:
+	dh_auto_clean
+	rm -f mapdamage/_version.py

mapdamage/Rscripts/stats/function.R

@@ -595,8 +595,8 @@ calcSampleStats <- function(da,X){
                       pos=da[,"Pos"],
                       mea=apply(X,1,mean),
                       med=apply(X,1,median),
-                      loCI=apply(X,1,quantile,c(0.025)),
-                      hiCI=apply(X,1,quantile,c(0.975))
+                      loCI=apply(X,1,quantile,c(0.025), na.rm = TRUE),
+                      hiCI=apply(X,1,quantile,c(0.975), na.rm = TRUE)
                       ))
 }
 

mapdamage/mp_test.py

+import unittest
+import subprocess
+
+
+def read_nocomment(dnacomp):
+    with open(dnacomp, 'r') as f:
+        a = f.readlines()
+        return [x for x in a if not x.startswith('#')]
+
+class testCases(unittest.TestCase):
+
+    def test_no_stats(self):
+        """regular mapDamage"""
+        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"))
+
+if  __name__=='__main__':
+    unittest.main()

mapdamage/rescale_test.py

-import unittest
-import optparse
-from . import rescale
-import pysam
-import filecmp
-
-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,
-        "quiet":True
-        })
-
-
-class testPairedFile(unittest.TestCase):
-    """Tests if rescaling of a paired end file"""
-    def test_paired_end_file(self):
-        """Test, paired end SAM file"""
-        # here is the example
-        # >ref
-        # 5p CGA AAA CGA 3p
-        # 3p GCT TTT GCT 5p
-        # >r001/1
-        # CGA
-        # >r001/2
-        # GCA
-        # The sam file looks like this
-        #@HD VN:1.5 SO:coordinate
-        #@SQ SN:ref LN:9
-        #r001 163 ref 1 30 3M = 7 9 CGA III MR:f:0 #the normal ones
-        #r001 83 ref 7 30 3M = 1 -9 TCG III MR:f:0
-        #r002 163 ref 1 30 3M = 7 9 TGA III MR:f:0.9 #With one dam subs
-        #r002 83 ref 7 30 3M = 1 -9 CAA III MR:f:0.009
-        #r003 83 ref 1 30 3M = 7 9 TGA III #The reverse complement, should not rescale (thus no flags)
-        #r003 163 ref 7 30 3M = 1 -9 CAA III
-        # 
-        #hand calculated the correct rescaled sam file in pe_rescaled_correct.sam
-        options = mock_options("rescale_test/pe_test/pe.sam","rescale_test/pe_test/pe_rescaled.sam","rescale_test/pe_test/pe_output/")
-        ref = pysam.Fastafile("rescale_test/pe_test/ref.fa")
-        rescale.rescale_qual(ref,options,debug=True)
-        self.assertTrue(filecmp.cmp("rescale_test/pe_test/pe_rescaled.sam","rescale_test/pe_test/pe_rescaled_correct.sam"))
-
-class testCases(unittest.TestCase):
-    """Various cases that failed"""
-    def test_longalignshortread(self):
-        """Check if fails on an aligment longer than the read"""
-        prefix="rescale_test/long_align/"
-        options = mock_options(prefix+"pe.sam",prefix+"pe_out.sam",prefix+"pe_output/")
-        ref = pysam.Fastafile("rescale_test/pe_test/ref.fa")
-        rescale.rescale_qual(ref,options,debug=True)
-        #Should run without an error
-
-if  __name__=='__main__':
-    unittest.main()

mapdamage/rescale_test/long_align/pe.sam

-@HD	VN:1.5	SO:coordinate
-@SQ	SN:ref	LN:9
-r001	163	ref	1	30	3M2D	=	7	9	CGA	III

mapdamage/rescale_test/long_align/pe_output/Stats_out_MCMC_correct_prob.csv

-"","Position","C.T","G.A"
-"1",1,0.9,0
-"2",2,0.009,0
-"3",-2,0,0.09
-"4",-1,0,0.9

mapdamage/rescale_test/long_align/ref.fa

->ref
-CGAAAACGA

mapdamage/rescale_test/pe_test/pe.sam

-@HD	VN:1.5	SO:coordinate
-@SQ	SN:ref	LN:9
-r001	163	ref	1	30	3M	=	7	9	CGA	III
-r001	83	ref	7	30	3M	=	1	-9	CGA	III
-r002	163	ref	1	30	3M	=	7	9	TGA	III
-r002	83	ref	7	30	3M	=	1	-9	CAA	III
-r003	83	ref	1	30	3M	=	7	9	TGA	III
-r003	163	ref	7	30	3M	=	1	-9	CAA	III

mapdamage/rescale_test/pe_test/pe_output/Stats_out_MCMC_correct_prob.csv

-"","Position","C.T","G.A"
-"1",1,0.9,0
-"2",2,0.009,0
-"3",-2,0,0.09
-"4",-1,0,0.9

mapdamage/rescale_test/pe_test/pe_rescaled_correct.sam

-@HD	VN:1.5	SO:coordinate
-@SQ	SN:ref	LN:9
-r001	163	ref	1	30	3M	=	7	9	CGA	III	MR:f:0
-r001	83	ref	7	30	3M	=	1	-9	CGA	III	MR:f:0
-r002	163	ref	1	30	3M	=	7	9	TGA	!II	MR:f:0.9
-r002	83	ref	7	30	3M	=	1	-9	CAA	I5I	MR:f:0.009
-r003	83	ref	1	30	3M	=	7	9	TGA	III
-r003	163	ref	7	30	3M	=	1	-9	CAA	III

mapdamage/rescale_test/pe_test/ref.fa

->ref
-CGAAAACGA