From 050f70a1da5457b100e08a9110d304c0197d8a0d Mon Sep 17 00:00:00 2001
From: Tobias Rausch <rauschtobi@gmail.com>
Date: Tue, 11 Apr 2023 13:11:13 +0200
Subject: [PATCH] fix bg noise

---
 src/profile.h | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/profile.h b/src/profile.h
index 816d749..e335e3e 100644
--- a/src/profile.h
+++ b/src/profile.h
@@ -29,7 +29,9 @@ namespace tracy
     p.resize(boost::extents[6][sz]);   // 'A', 'C', 'G', 'T', 'N', '-'
     for(int32_t j = trimleft; j < (trimleft + sz); ++j) {
       float totalsig = 0;
+      float allBaseSig = 0;
       for(uint32_t k = 0; k<4; ++k) {
+	allBaseSig += tr.traceACGT[k][bc.bcPos[j]];
 	if (_inBaseCalled(k, bc.primary[j], bc.secondary[j])) totalsig += tr.traceACGT[k][bc.bcPos[j]];
       }
       p[4][j-trimleft] = 0;
@@ -40,6 +42,11 @@ namespace tracy
 	for(uint32_t k = 0; k<4; ++k) {
 	  if (_inBaseCalled(k, bc.primary[j], bc.secondary[j])) p[k][j-trimleft] = ((float) (tr.traceACGT[k][bc.bcPos[j]]) / totalsig);
 	}
+	// Sometimes basecall signal is tiny fraction of all base signals (e.g., missing peaks in signal ramps)
+	float normfac = totalsig / allBaseSig;
+	for(uint32_t k = 0; k<4; ++k) {
+	  p[k][j-trimleft] = normfac * p[k][j-trimleft] + (1 - normfac) * 0.25;
+	}
       }
     }
   }