attempts to remove ladder peaks that resulf from oversaturated peaks …

…in other channels.

.vscode/settings.json

@@ -0,0 +1,15 @@
+{
+    "python.testing.unittestArgs": [
+        "-v",
+        "-s",
+        "./tests",
+        "-p",
+        "*test.py"
+    ],
+    "python.testing.pytestEnabled": false,
+    "python.testing.unittestEnabled": true,
+    "[python]": {
+        "editor.defaultFormatter": "ms-python.python"
+    },
+    "python.formatting.provider": "none"
+}

fatools/lib/fautil/algo2.py

@@ -24,6 +24,10 @@
 
 import attr
 
+import logging, os
+LOGLEVEL = os.environ.get("LOGLEVEL", "WARNING").upper()
+logging.basicConfig(level=LOGLEVEL)
+
 class LadderMismatchException(Exception):
     """Raised when number of peaks in ladder channel not equal to number of ladder steps."""
     pass
@@ -237,15 +241,19 @@ def call_peaks( channel, params, func, min_rtime, max_rtime ):
         if is_verbosity(4):
             print(allele)
 
-def align_peaks(channel, params, ladder, anchor_pairs=None):
+def align_peaks(channel, params, ladder, anchor_pairs=None, saturated_peak_rtimes=[]):
     """
     returns (score, rss, dp, aligned_peak_number)
     """
 
     alleles = channel.get_alleles()
 
     if (len(alleles) != len(ladder['sizes'])):
-        raise LadderMismatchException( ("alleles not same length as ladder for file: %s!") % channel.fsa.filename)
+        # attempt to filter out cross talk peaks from saturated signal in other channels
+        logging.debug(f"observed rtimes:{[a.rtime for a in alleles]}, saturated_rtimes:{saturated_peak_rtimes}")
+        alleles = [ a for a in alleles if not any(np.isclose(a.rtime,s, atol= 2) for s in saturated_peak_rtimes)] 
+        if (len(alleles) != len(ladder['sizes'])):  
+            raise LadderMismatchException( ("alleles not same length as ladder for file: %s!") % channel.fsa.filename)
 
     if anchor_pairs:
         return align_pm( alleles, ladder, anchor_pairs)

fatools/lib/fautil/mixin2.py

@@ -1,4 +1,5 @@
 
+import itertools
 from fatools.lib.fautil import algo2 as algo
 from fatools.lib.utils import cout, cerr, cexit, is_verbosity
 from fatools.lib import const
@@ -87,8 +88,7 @@ def preannotate(self, parameters):
         algo.preannotate_peaks(self, params)
 
     # ChannelMixIn align method
-    def align(self, parameters, ladder=None, anchor_pairs=None):
-
+    def align(self, parameters, ladder=None, anchor_pairs=None, saturated_peak_rtimes=[]):
         # sanity checks
         if self.marker.code != 'ladder':
             raise RuntimeError('E: align() must be performed on ladder channel!')
@@ -101,7 +101,7 @@ def align(self, parameters, ladder=None, anchor_pairs=None):
                                             ladder['strict'], ladder['relax'] )
 
         start_time = time.process_time()
-        result = algo.align_peaks(self, parameters, ladder, anchor_pairs)
+        result = algo.align_peaks(self, parameters, ladder, anchor_pairs, saturated_peak_rtimes)
         dpresult = result.dpresult
         fsa = self.fsa
         fsa.z = dpresult.z
@@ -295,7 +295,7 @@ def align(self, parameters=None):
         self.scan(parameters)
 
         ladder = self.get_ladder_channel()
-        ladder.align(parameters)
+        ladder.align(parameters, None, None, self.get_saturated_rtimes())
 
     # FSAMixIn call method
     def call(self, parameters):
@@ -337,6 +337,13 @@ def get_ladder_channel(self):
             if c.marker.code == 'ladder':
                 return c
         raise RuntimeError('E: ladder channel not found')
+
+    def get_saturated_rtimes(self):
+        saturation_threshold = 29000
+        #makes a flat list of alleles from all channels
+        alleles = itertools.chain(*[c.get_alleles() for c in self.channels])
+        saturated_rtimes = [a.rtime for a in alleles if a.rfu >= saturation_threshold]
+        return saturated_rtimes
 
 
 

tests/channel_test.py

@@ -0,0 +1,42 @@
+import unittest
+
+from fatools.lib.fautil.mixin2 import FSAMixIn
+
+
+class GetSaturatedRtimesTest(unittest.TestCase):
+
+    #mocks up a channels and alleles
+    class Channel:
+        def __init__(self, alleles):
+            self.alleles = alleles
+        def get_alleles(self):
+            return self.alleles
+
+    class Allele:
+        def __init__(self, rfu, rtime):
+            self.rfu = rfu
+            self.rtime = rtime
+
+    high_alleles = [Allele(29000,100),Allele(29001,200), Allele(35000,300), Allele(40000,400)]
+    low_alleles = [Allele(28999,100),Allele(3000,200), Allele(1000,300), Allele(500,400)]  
+
+    def setUp(self):
+        self.fsa = FSAMixIn()
+
+    def test_get_saturated_rtimes_empty_channels(self):
+        self.fsa.channels = []
+        saturated_rtimes = self.fsa.get_saturated_rtimes()
+        self.assertEqual(saturated_rtimes, [])
+
+    def test_get_saturated_rtimes_no_saturated_alleles(self):
+        self.fsa.channels = [self.Channel(self.low_alleles)]
+        saturated_rtimes = self.fsa.get_saturated_rtimes()
+        self.assertEqual(saturated_rtimes, [])
+
+    def test_get_saturated_rtimes_some_saturated_allele(self):
+        self.fsa.channels = [self.Channel(self.low_alleles),self.Channel(self.high_alleles)]
+        saturated_rtimes = self.fsa.get_saturated_rtimes()
+        self.assertEqual(saturated_rtimes, [100,200,300,400])
+
+if __name__ == "__main__":
+    unittest.main()