New upstream version 2.1.0+dfsg

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

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

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/).

`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 ressources 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/)

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

### Citation

### Contact

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

aginolhac at snm.ku.dk, MSchubert at snm.ku.dk or jonsson.hakon at gmail.com.

ginolhac at gmail.com, mikkelsch at gmail.com or jonsson.hakon at gmail.com.

#!/usr/bin/python

#!/usr/bin/env python3

# -*- coding: utf-8 -*-

from __future__ import print_function

import logging

import random

import time

URL = "http://code.google.com/p/pysam/"

try:

    __import__(MODULE)

except ImportError, e:

except ImportError as e:

    sys.stderr.write("Error: Could not import required module '%s':\n\t- %s\n" % (MODULE, e))

    sys.stderr.write("       If module is not installed, please download from '%s'.\n" % URL)

    sys.stderr.write("       A local install may be performed using the following command:\n")

            if index >= downsample_to:

                continue

        sample[index] = record

    return filter(None, sample)

    return [_f for _f in sample if _f]

def _read_bamfile(bamfile, options):

    in mapdamage.parseoptions

    """

    path_to_mm = None

    for nr,arg in zip(xrange(len(sys.argv)),sys.argv):

    for nr,arg in zip(list(range(len(sys.argv))),sys.argv):

        if arg.startswith("--mapdamage-modules"):

            try:

                if "=" in arg:

        in_bam = pysam.Samfile(options.filename)

    reflengths = dict(zip(in_bam.references, in_bam.lengths))

    reflengths = dict(list(zip(in_bam.references, in_bam.lengths)))

    # check if references in SAM/BAM are the same in the fasta reference file

    fai_lengths = mapdamage.seq.read_fasta_index(options.ref + ".fai")

    if not fai_lengths:

write.mutation.table <- function(tbl, end, mismatch, filename)

  {

    columns <- c("pos", sprintf("5p%s", mismatch))

    columns <- c("pos", sprintf("%s%s", end, mismatch))

    tbl[is.na(tbl)] <- 0 # Python doesn't like NAs

    write.table(tbl[tbl$End == end, c("Pos", mismatch)],

                row.names = FALSE, col.names = columns,

    REsStats <- calcSampleStats(da,REs)

    #Plotting the posterior predictive intervals

    p <- ggplot()+

         geom_point(aes(x,mea,colour="C->T",aes_string="subs"),data=CTsStats)+

         geom_point(aes(x,mea,colour="C->T"),data=CTsStats)+

         geom_point(aes(x,mea,colour="G->A"),data=GAsStats)+

         geom_point(aes(x,mea,colour="Others"),data=REsStats)+

         geom_errorbar(aes(x=x,y=med,ymin=loCI,ymax=hiCI,color="C->T"),data=CTsStats)+

         geom_errorbar(aes(x=x,y=med,ymin=loCI,ymax=hiCI,color="G->A"),data=GAsStats)+

         geom_errorbar(aes(x=x,y=med,ymin=loCI,ymax=hiCI,color="Others"),data=REsStats)+

         geom_errorbar(aes(x=x,ymin=loCI,ymax=hiCI,color="C->T"),data=CTsStats)+

         geom_errorbar(aes(x=x,ymin=loCI,ymax=hiCI,color="G->A"),data=GAsStats)+

         geom_errorbar(aes(x=x,ymin=loCI,ymax=hiCI,color="Others"),data=REsStats)+

         geom_line(aes(oneBased,C.T/C),color="red",data=data.frame(da))+

         geom_line(aes(oneBased,G.A/G),color="green",data=data.frame(da))+

         geom_line(aes(oneBased,((A.C+A.G+A.T)/A+(C.A+C.G)/C+(G.C+G.T)/G+(T.A+T.C+T.G)/T)/10),color="blue",data=data.frame(da))+

         ylab("Substitution rate")+

         xlab("Relative position")+

         scale_x_continuous(breaks=bres,labels=labs)+

         labs(colour = "Subs. type")+

         ggtitle("Posterior prediction intervals")

         labs(y = "Substitution rate",

              x = "Relative position",

              colour = "Subs. type",

              title = "Posterior prediction intervals")+

         scale_x_continuous(breaks=bres,labels=labs)

    if (output$cu_pa$use_bw_theme){

        p <- p+theme_bw()

    }

    rownames(summStat)[3] <- "Acceptance ratio"

    write.csv(summStat,paste(filename,"_summ_stat.csv",sep=""))

}

#The main work flow of the package

#Load the libraries

suppressMessages(library(inline))  #Already checked the libraries 

suppressMessages(library(ggplot2)) #thus ignoring any messages from them

suppressMessages(library(Rcpp))

suppressMessages(library(gam)) 

suppressMessages(library(RcppGSL))

suppressPackageStartupMessages(library(inline))  #Already checked the libraries

suppressPackageStartupMessages(library(ggplot2)) #thus ignoring any messages from them

suppressPackageStartupMessages(library(Rcpp))

suppressPackageStartupMessages(library(gam))

suppressPackageStartupMessages(library(RcppGSL))

#Miscellaneous functions

source(paste(path_to_mapDamage_stats,"function.R",sep=""))