Merge pull request #144 from gymreklab/develop

Fix bug #143

RELEASE_NOTES.rst

@@ -1,3 +1,11 @@
+4.0.1
+-----
+
+Bug fixes:
+
+* https://github.com/gymreklab/TRTools/issues/143 Fix HipstrMinSuppReads filter when
+  there are called samples but none have ALLREADS
+
 4.0.0
 -----
 

setup.py

@@ -15,7 +15,7 @@
 curdir = os.path.abspath(os.path.dirname(__file__))
 MAJ = 4
 MIN = 0
-REV = 0
+REV = 1
 VERSION = '%d.%d.%d' % (MAJ, MIN, REV)
 with open(os.path.join(curdir, 'trtools/version.py'), 'w') as fout:
         fout.write(

trtools/dumpSTR/filters.py

@@ -494,15 +494,27 @@ def __call__(self, record: trh.TRRecord):
             return np.full((record.GetNumSamples()), np.nan)
 
         if "ALLREADS" not in record.format:
-            sample_filter = np.zeros((record.GetNumSamples()), dtype=float)
-            return sample_filter
+            return np.zeros((record.GetNumSamples()), dtype=float)
         samples_to_check = (called_samples &
                           (record.format["ALLREADS"] != '') &
                           (record.format["ALLREADS"] != '.'))
 
-        delim = "|"
+        if not np.any(samples_to_check):
+            # all samples were either not called or were missing ALLREADS
+            # say that we filtered all the samples which were missing ALLREADS
+            sample_filter = np.full((record.GetNumSamples()), np.nan)
+            sample_filter[called_samples] = 0
+            return sample_filter
+
         # Going to assume that either all samples are phased or none are
-        if "/" in record.format["GB"][samples_to_check][0]: delim = "/"
+        first_gb = record.format["GB"][samples_to_check][0]
+        if "/" in first_gb:
+            delim = "/"
+        elif "|" in first_gb:
+            delim = '|'
+        else:
+            raise ValueError("Cant't identify phasing char ('|' or '/') in GB field")
+
         gb = np.char.split(record.format["GB"][samples_to_check], delim)
         gb = np.stack(gb).astype(int)
         # Format allreads like a python dict literal so we can interpret it

trtools/dumpSTR/tests/test_filters.py

@@ -311,12 +311,12 @@ class TestRec(DummyRecBase):
         def __init__(self):
             super().__init__()
             self.format['ALLREADS'] = np.array([
-                '0|23;1|123;2|5', '0|165;1|23;2|7',
+                '0|23;1|123;2|5', '0|15;1|23;2|7',
                 '0|23;1|444;2|12', '0|23;1|32;2|66',
                 '0|867;1|23;2|13', '0|848;1|92;2|483',
                 '', '', '.'])
-            self.format['GB'] = np.array(['1|1', '1|1', '0|2', '0|2', '0|2',
-                                          '0|2', '', '', '.'])
+            self.format['GB'] = np.array(['1|1', '1|1', '1|2', '2|1', '2|0',
+                                          '0|2', '1|1', '0|0', '1|0'])
         def GetNumSamples(self):
             return 9
         def GetCalledSamples(self):
@@ -330,13 +330,37 @@ def GetCalledSamples(self):
     assert np.isnan(out[0])
     assert out[1] == 23
     assert out[2] == 12
-    assert out[3] == 23
+    assert out[3] == 32
     assert out[4] == 13
     assert np.isnan(out[5])
     assert out[6] == 0 # If ALLREADS is missing, filter
     assert np.isnan(out[7]) # don't apply filter to nocalls
     assert np.isnan(out[8]) # don't apply filter to nocalls
 
+def test_HipstrMinSuppReads_no_called_samples_with_reads(tmpdir):
+    class TestRec(DummyRecBase):
+        def __init__(self):
+            super().__init__()
+            self.format['ALLREADS'] = np.array([
+                '0|23;1|123;2|5', '0|15;1|23;2|7',
+                '0|23;1|444;2|12', '0|23;1|32;2|66',
+                '0|867;1|23;2|13', '0|848;1|92;2|483',
+                '', '', '.'])
+            self.format['GB'] = np.array(['1|1', '1|1', '1|2', '2|1', '2|0',
+                                          '0|2', '1|1', '0|0', '1|0'])
+        def GetNumSamples(self):
+            return 9
+        def GetCalledSamples(self):
+            return np.array([False, False, False, False, False, False, True, False, True])
+
+    args = base_argparse(tmpdir)
+    args.hipstr_min_supp_reads = 50
+    call_filters = BuildCallFilters(args)
+    assert len(call_filters) == 1
+    out = call_filters[0](TestRec())
+    assert out.shape == (9,)
+    assert np.all(out[[6,8]] == 0)
+    assert np.all(np.isnan(out[[0,1,2,3,4,5,7]]))
 
 def test_HipstrDP(tmpdir):
     class TestRec(DummyRecBase):

trtools/version.py

@@ -1,4 +1,4 @@
 
 # THIS FILE IS GENERATED FROM SETUP.PY
-version = '4.0.0'
+version = '4.0.1'
 __version__ = version