This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository mapdamage.
commit d046a9d5a9deb7c2a41abd7caade686e8dd1e8cc
Author: Andreas Tille <ti...@debian.org>
Date: Mon Jun 19 09:36:00 2017 +0200
New upstream version 2.0.8+dfsg
---
README.md | 3 +++
mapdamage/parseoptions.py | 32 +++++++++++++++++++++++++-------
mapdamage/rescale.py | 23 ++++++++++++++++-------
mapdamage/version.py | 3 +--
setup.py | 2 +-
5 files changed, 46 insertions(+), 17 deletions(-)
diff --git a/README.md b/README.md
index f908dba..a0915b8 100644
--- a/README.md
+++ b/README.md
@@ -9,7 +9,10 @@ among ancient DNA sequencing reads generated by
Next-Generation Sequencing platf
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/)
### Citation
diff --git a/mapdamage/parseoptions.py b/mapdamage/parseoptions.py
index 7102cf2..bf5020f 100644
--- a/mapdamage/parseoptions.py
+++ b/mapdamage/parseoptions.py
@@ -128,20 +128,27 @@ def options():
help="How long sequence to use from each side [%default]", type =
int, default=12, action="store")
group3.add_option("--stats-only", dest="stats_only", help="Run only
statistical estimation from a valid result folder", \
default=False, action="store_true")
- group3.add_option("--rescale", dest="rescale", help="Rescale the quality
scores in the BAM file using the output from the statistical estimation", \
+ group3.add_option("--no-stats", help="Disabled statistical estimation,
active by default", default=False, action="store_true")
+ group3.add_option("--check-R-packages", help="Check if the R modules are
working", default=False, action="store_true")
+ parser.add_option_group(group3)
+
+ group4 = OptionGroup(parser,"Options for rescaling of BAM files")
+ group4.add_option("--rescale", dest="rescale", help="Rescale the quality
scores in the BAM file using the output from the statistical estimation",
default=False, action="store_true")
- group3.add_option("--rescale-only", dest="rescale_only", help="Run only
rescaling from a valid result folder", \
+ group4.add_option("--rescale-only", dest="rescale_only", help="Run only
rescaling from a valid result folder",
default=False, action="store_true")
- group3.add_option("--rescale-out", dest="rescale_out", help="Write the
rescaled BAM to this file", \
+ group4.add_option("--rescale-out", dest="rescale_out", help="Write the
rescaled BAM to this file",
default=None, action="store")
- group3.add_option("--no-stats", help="Disabled statistical estimation,
active by default", default=False, action="store_true")
- group3.add_option("--check-R-packages", help="Check if the R modules are
working", default=False, action="store_true")
+ group4.add_option("--rescale-length-5p", dest="rescale_length_5p",
+ help="How many bases to rescale at the 5' termini;
defaults to --seq-length.", type=int, action="store")
+ group4.add_option("--rescale-length-3p", dest="rescale_length_3p",
+ help="How many bases to rescale at the 5' termini;
defaults to --seq-length.", type=int, action="store")
+ parser.add_option_group(group4)
- parser.add_option_group(group3)
#Parse the arguments
(options, args) = parser.parse_args()
-
+
# check python version
if not check_py_version():
return None
@@ -238,6 +245,17 @@ def options():
sys.stderr.write("Error, %s does not exist while
plot/stats/rescale only was used\n" % options.folder)
return None
+ if options.rescale_length_3p is None:
+ options.rescale_length_3p = options.seq_length
+ elif not (0 <= options.rescale_length_3p <= options.seq_length):
+ parser.error("--rescale-length-3p must be less than or equal to "
+ "--seq-length and greater than zero")
+
+ if options.rescale_length_5p is None:
+ options.rescale_length_5p = options.seq_length
+ elif not (0 <= options.rescale_length_5p <= options.seq_length):
+ parser.error("--rescale-length-5p must be less than or equal to "
+ "--seq-length and greater than zero")
# check if the Rscript executable is present on the system
if not whereis('Rscript'):
diff --git a/mapdamage/rescale.py b/mapdamage/rescale.py
index 789d21e..ca71c3e 100644
--- a/mapdamage/rescale.py
+++ b/mapdamage/rescale.py
@@ -17,7 +17,7 @@ def phred_char_to_pval(ch):
"""Transforming ASCII character in the Phred scale to the error rate"""
return 10**(-(float(ord(ch))-float(33))/10)
-def get_corr_prob(folder):
+def get_corr_prob(folder, rescale_length_5p, rescale_length_3p):
"""
Reads the damage probability correction file, returns a
dictionary with this structure
@@ -36,6 +36,12 @@ def get_corr_prob(folder):
(folder, fi_handle.line_num, corr_prob[line["Position"]]))
else:
corr_prob[int(line["Position"])] = {'C.T':float(line["C.T"]),
'G.A':float(line["G.A"])}
+
+ # Exclude probabilities for positions outside of user-specified region
+ for key in 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"), \
@@ -269,11 +275,12 @@ def rescale_qual_read(bam, read, ref, corr_prob,subs,
debug = False,direction="b
read.qual = new_qual
# truncate this to 5 digits
number_of_rescaled_bases = float("%.5f" % number_of_rescaled_bases)
- # check if the read has a MR tag
- for tag in read.tags:
- if tag[0] == "MR":
- raise SystemExit("Read: %s already has a MR tag, can't rescale" %
(read))
- read.tags = read.tags + [("MR:f",number_of_rescaled_bases)]
+
+ if read.has_tag("MR"):
+ raise SystemExit("Read: %s already has a MR tag, can't rescale" % read)
+
+ read.set_tag("MR", number_of_rescaled_bases, 'f')
+
return read
@@ -300,7 +307,9 @@ def rescale_qual(ref, options,debug=False):
else:
write_mode = "wb"
bam_out = pysam.Samfile(options.rescale_out,write_mode, template = bam)
- corr_prob = get_corr_prob(options.folder)
+ corr_prob = get_corr_prob(options.folder,
+ rescale_length_5p=options.rescale_length_5p,
+ rescale_length_3p=options.rescale_length_3p)
subs = initialize_subs()
first_pair = True
number_of_non_proper_pairs = 0
diff --git a/mapdamage/version.py b/mapdamage/version.py
index 4e95d9f..5d2b8ff 100644
--- a/mapdamage/version.py
+++ b/mapdamage/version.py
@@ -2,5 +2,4 @@
try:
from mapdamage._version import __version__
except ImportError:
- __version__ = "2.0.6"
-
+ __version__ = "2.0.8"
diff --git a/setup.py b/setup.py
index b166396..d0c70ca 100644
--- a/setup.py
+++ b/setup.py
@@ -55,7 +55,7 @@ class compileInstall(DistutilsInstall):
setup(
cmdclass={'install': compileInstall},
name='mapdamage',
- version='2.0',
+ version='2.0.8',
author='Aurélien Ginolhac, Mikkel Schubert, Hákon Jónsson',
author_email='mschub...@snm.ku.dk, jonsson.ha...@gmail.com',
packages=['mapdamage'],
--
Alioth's /usr/local/bin/git-commit-notice on
/srv/git.debian.org/git/debian-med/mapdamage.git
_______________________________________________
debian-med-commit mailing list
debian-med-commit@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
