remove pybigtools backend & bug fixes

Author: yztxwd

pyproject.toml

@@ -19,7 +19,6 @@ dependencies = [
     'pybigwig>=0.3.0',
     'torch>=1.10.0',
     'webdataset>=0.2.0',
-    'pybigtools>=0.2.5, <0.3'
 ]
 license = "MIT AND (Apache-2.0 OR BSD-2-Clause)"
 keywords = ["interpretation",  "attribution", "concept", "genomics", "deep learning"]

seqchromloader/__init__.py

@@ -8,6 +8,6 @@
 
 import logging
 def mute_warning():
-    loader_logger.setLevel(logging.ERROR)
-    writer_logger.setLevel(logging.ERROR)
-    utils_logger.setLevel(logging.ERROR)
+    loader_logger.setLevel(logging.INFO)
+    writer_logger.setLevel(logging.INFO)
+    utils_logger.setLevel(logging.INFO)

seqchromloader/config.py

@@ -3,6 +3,6 @@
 def set_bigwig_backend(backend):
     global BIGWIG_BACKEND
 
-    assert backend in ['pyBigWig', 'pybigtools', 'memmap'], "Backend options are pyBigWig, memmap or pybigtools!"
+    assert backend in ['pyBigWig', 'memmap'], "Backend options are pyBigWig or memmap!"
 
     BIGWIG_BACKEND = backend

seqchromloader/loader.py

@@ -7,16 +7,13 @@
 import math
 import logging
 import torch
-import random
 import pysam
 import pyfaidx
-import pybigtools
 import numpy as np
 import pandas as pd
 import webdataset as wds
-from math import sqrt, ceil
 from itertools import islice
-from torch.utils.data import Dataset, IterableDataset, DataLoader
+from torch.utils.data import IterableDataset, DataLoader
 from pybedtools import BedTool
 
 from seqchromloader import utils, config

seqchromloader/utils.py

@@ -9,7 +9,6 @@
 import logging
 import pysam
 import pyBigWig
-import pybigtools
 from Bio import motifs
 from pyfaidx import Fasta
 from multiprocessing import Pool
@@ -25,9 +24,6 @@ def __init__(self, bw_path, memmap_dir=None, backend='pyBigWig'):
         if backend == 'pyBigWig':
             self.bw = pyBigWig.open(bw_path)
             self._chroms = self.bw.chroms()
-        elif backend == 'pybigtools':
-            self.bw = pybigtools.open(bw_path)
-            self._chroms = self.bw.chroms()
         elif backend == 'memmap':
             # load all chromosomes using np.memmap
             if memmap_dir is None:
@@ -36,8 +32,8 @@ def __init__(self, bw_path, memmap_dir=None, backend='pyBigWig'):
             self.memmap_dir = memmap_dir
             os.makedirs(self.memmap_dir, exist_ok=True)
 
-            # open bigwig via pybigtools
-            self.bw = pybigtools.open(bw_path)
+            # open bigwig via pyBigWig
+            self.bw = pyBigWig.open(bw_path)
 
             self.arrays = {}
             self._chroms = {}
@@ -50,7 +46,7 @@ def __init__(self, bw_path, memmap_dir=None, backend='pyBigWig'):
                     logger.info(f"[memmap] dumping {chrom}...")
 
                     # load full chromosome
-                    data = self.bw.values(chrom, 0, length, missing=0.0)
+                    data = np.nan_to_num(self.bw.values(chrom, 0, length, numpy=True), nan=0.0).astype(np.float32)
 
                     # save
                     np.save(npy_path, data)
@@ -65,34 +61,18 @@ def __init__(self, bw_path, memmap_dir=None, backend='pyBigWig'):
         self.backend = backend
 
     def intervals(self, chrom):
-        if self.backend == 'pyBigWig':
-            return self.bw.intervals(chrom)
-        else:
-            return self.bw.records(chrom)
+        return self.bw.intervals(chrom)
 
     def stats(self, chrom, type='mean', exact=True):
         if self.backend == 'pyBigWig':
-            return self.bw.stats(chrom, type=type, exact=exact)[0]
-        elif self.backend == 'pybigtools':
-            return self.bw.values(chrom, missing=np.nan, bins=1, exact=exact, summary=type)[0].item()
+            return np.float32(self.bw.stats(chrom, type=type, exact=exact)[0])
         elif self.backend == 'memmap':
-            arr = self.arrays[chrom]
-            if type == 'mean':
-                return float(np.mean(arr))
-            elif type == 'max':
-                return float(np.max(arr))
-            elif type == 'min':
-                return float(np.min(arr))
-            elif type == 'std':
-                return float(np.std(arr))
-            else:
-                raise NotImplementedError(f"stat {type} not implemented for memmap backend")
+            # use pyBigWig to obtain stats to avoid load full chromosome arrays into memory
+            return np.float32(self.bw.stats(chrom, type=type, exact=exact)[0])
 
     def values(self, chrom, start, end, missing=0):
         if self.backend == 'pyBigWig':
-            return np.nan_to_num(self.bw.values(chrom, start, end), nan=missing).astype(np.float32)
-        elif self.backend == 'pybigtools':
-            return self.bw.values(chrom, start, end, missing=missing).astype(np.float32)
+            return np.nan_to_num(self.bw.values(chrom, start, end, numpy=True), nan=missing).astype(np.float32)
         elif self.backend == 'memmap':
             arr = self.arrays[chrom]
 

tests/test_backend_memmap.py

@@ -0,0 +1,8 @@
+from test_writer_loader import Test
+import unittest, seqchromloader
+
+class TestBackendMemmap(Test, unittest.TestCase):
+    seqchromloader.config.set_bigwig_backend("memmap")
+
+if __name__ == "__main__":
+    unittest.main(verbosity=2)

tests/test_backend_pybigwig.py

@@ -0,0 +1,9 @@
+from test_writer_loader import Test
+import unittest, seqchromloader
+
+class TestBackendPyBigWig(Test, unittest.TestCase):
+    seqchromloader.config.set_bigwig_backend("pyBigWig")
+
+if __name__ == "__main__":
+    unittest.main(verbosity=2)
+

tests/test_writer_loader.py

@@ -20,9 +20,9 @@
 
 import seqchromloader
 
-seqchromloader.config.set_bigwig_backend('pybigtools')
+#seqchromloader.config.set_bigwig_backend('memmap')
 
-class Test(unittest.TestCase):
+class Test:
     def setUp(self) -> None:
         pass
 
@@ -40,6 +40,10 @@ def tearDownClass(cls) -> None:
     def assertIsFile(self, path):
         if not pl.Path(path).resolve().is_file():
             raise AssertionError("File does not exist: %s" % str(path))
+
+    def test_compute_mean_std_bw(self):
+        bw_f = "data/sample.bw"
+        mean, std = seqchromloader.utils.compute_mean_std_bigwig(bw_f)
 
     def test_get_genome_sizes(self):
         genome_sizes_nochr10=get_genome_sizes(genome="mm10", to_filter=["chr10"])
@@ -389,6 +393,3 @@ def test_chrom_transform(chrom):
 
 def test_target_transform(target):
     return target * 3
-
-if __name__ == "__main__":
-    unittest.main(verbosity=2)