From b3a88395ec191d99d2b5db30cb1c81bf8c828db7 Mon Sep 17 00:00:00 2001
From: fransua <serra.francois@gmail.com>
Date: Wed, 23 Nov 2016 23:28:59 +0100
Subject: [PATCH] BAM: decay correction adjusted.

---
 _pytadbit/_version.py         |  2 +-
 scripts/matrix_from_BAM.py    |  2 ++
 scripts/normalize_from_BAM.py | 25 ++++++++++++++++++-------
 3 files changed, 21 insertions(+), 8 deletions(-)

diff --git a/_pytadbit/_version.py b/_pytadbit/_version.py
index 6aa4bcd0..5bdf0da2 100644
--- a/_pytadbit/_version.py
+++ b/_pytadbit/_version.py
@@ -1 +1 @@
-__version__ = "3DAROC_2016.14"
\ No newline at end of file
+__version__ = "3DAROC_2016.16"
\ No newline at end of file
diff --git a/scripts/matrix_from_BAM.py b/scripts/matrix_from_BAM.py
index 74bdafd6..6f243b5a 100644
--- a/scripts/matrix_from_BAM.py
+++ b/scripts/matrix_from_BAM.py
@@ -57,6 +57,8 @@ def read_bam_frag(inbam, filter_exclude, sections,
             pos1 = r.reference_start + 1
             crm2 = refs[r.mrnm]
             pos2 = r.mpos + 1
+            print r
+            print crm1, pos1, crm2, pos2
             try:
                 pos1 = sections[(crm1, pos1 / resolution)]
                 pos2 = sections[(crm2, pos2 / resolution)]
diff --git a/scripts/normalize_from_BAM.py b/scripts/normalize_from_BAM.py
index 9e7b45b4..ddf4f809 100644
--- a/scripts/normalize_from_BAM.py
+++ b/scripts/normalize_from_BAM.py
@@ -220,9 +220,16 @@ def read_bam(inbam, filter_exclude, resolution, min_count=2500,
     pool.join()
 
     # to correct biases
-    sumdec = sum(p.get() for p in procs)
-    target = sumdec / float(size * size * factor)
-    diasum = dict([(d, diasum[d] * target) for d in diasum])
+    sumdec = {}
+    for proc in procs:
+        for k, v in proc.get().iteritems():
+            try:
+                sumdec[k] += v
+            except KeyError:
+                sumdec[k]  = v
+    target = dict([(d, sumdec[d] / float((total - d) * factor * 2))
+                   for d in sumdec])
+    diasum = dict([(d, diasum[d] * target[d]) for d in diasum])
 
     return biases, diasum, badcol
 
@@ -236,12 +243,16 @@ def sum_nrm_matrix(fname, biases):
 
 def sum_dec_matrix(fname, biases, decay):
     dico = load(open(fname))
-    sumdec = 0
-    for (i, j), v in dico.iteritems(): 
+    sumdec = {}
+    for (i, j), v in dico.iteritems():
+        k = abs(i-j)
         try:
-            sumdec += v / (biases[i] * biases[j] * decay[abs(i-j)])
+            sumdec[k] += v / (biases[i] * biases[j] * decay[k])
         except KeyError:
-            pass
+            try:
+                sumdec[k]  = v / (biases[i] * biases[j] * decay[k])
+            except KeyError:
+                pass
     os.system('rm -f %s' % (fname))
     return sumdec