add broad peak calling to macs3 (close #365)

snapatac2-core/Cargo.toml

@@ -15,7 +15,7 @@ anyhow = "1.0"
 bstr = "1.0"
 byteorder = "1.0"
 bigtools = { version = "0.5", features = ["read", "write"] }
-bed-utils = "0.6.2"
+bed-utils = "0.7.1"
 flate2 = "1.0"
 tokio = "1.34"
 hora = "0.1"

snapatac2-python/Cargo.toml

@@ -16,7 +16,7 @@ anndata = "0.4.2"
 anndata-hdf5 = "0.3"
 pyanndata = "0.4.1"
 anyhow = "1.0"
-bed-utils = "0.6.2"
+bed-utils = "0.7.1"
 flate2 = "1.0"
 itertools = "0.13"
 indicatif = "0.17"

snapatac2-python/pyproject.toml

@@ -31,7 +31,8 @@ dependencies = [
     'anndata >= 0.8.0, < 0.11.0',
     'kaleido',
     'multiprocess',
-    'MACS3 >= 3.0, < 3.1',
+    #'MACS3 >= 3.0, < 3.1',
+    "MACS3 @ git+https://github.com/kaizhang/MACS.git@1cceabc0420959e74ae610e34d4d71f88bf5f6a3",
     'natsort', 
     'numpy >= 1.16.0, < 2.0.0',
     'pandas >= 1.0, < 2.1.2',

snapatac2-python/python/snapatac2/tools/_call_peaks.py

@@ -11,6 +11,8 @@ def macs3(
     *,
     groupby: str | list[str] | None = None,
     qvalue: float = 0.05,
+    call_broad_peaks: bool = False,
+    broad_cutoff: float = 0.1,
     replicate: str | list[str] | None = None,
     replicate_qvalue: float | None = None,
     max_frag_size: int | None = None,
@@ -37,6 +39,17 @@ def macs3(
         `.obs[groupby]`. If None, peaks will be called for all cells.
     qvalue
         qvalue cutoff used in MACS3.
+    call_broad_peaks
+        If True, MACS3 will call broad peaks. The broad peak calling process
+        utilizes two distinct cutoffs to discern broader, weaker peaks (`broad_cutoff`)
+        and narrower, stronger peaks (`qvalue`), which are subsequently nested to
+        provide a detailed peak landscape. To conceptualize "nested" peaks, picture
+        a gene structure housing regions analogous to exons (strong peaks) and
+        introns coupled with UTRs (weak peaks). Please note that, if you only want to
+        call "broader" peak and not interested in the nested peak structure, please
+        simply use `qvalue` with weaker cutoff instead of using `call_broad_peaks` option.
+    broad_cutoff
+        qvalue cutoff used in MACS3 for calling broad peaks.
     replicate
         Replicate information. If provided, reproducible peaks will be called
         for each group.
@@ -119,8 +132,11 @@ def macs3(
     options.nolambda = nolambda
     options.smalllocal = 1000
     options.largelocal = 10000
-    options.call_summits = True
-    options.broad = False
+    options.call_summits = False if call_broad_peaks else True
+    options.broad = call_broad_peaks
+    if options.broad:
+        options.log_broadcutoff = log(broad_cutoff, 10) * -1
+
     options.fecutoff = 1.0
     options.d = extsize
     options.scanwindow = 2 * options.d
@@ -142,7 +158,7 @@ def _call_peaks(tags):
             tempfile.tempdir = tmpdirname  # Overwrite the default tempdir in MACS3
             merged, reps = _snapatac2.create_fwtrack_obj(tags)
             options.log_qvalue = log(qvalue, 10) * -1
-            logging.getLogger().setLevel(logging.CRITICAL + 1)
+            logging.getLogger().setLevel(logging.CRITICAL + 1) # temporarily disable logging
             peakdetect = PeakDetect(treat=merged, opt=options)
             peakdetect.call_peaks()
             peakdetect.peaks.filter_fc(fc_low = options.fecutoff)
@@ -157,7 +173,7 @@ def _call_peaks(tags):
                 peakdetect.peaks.filter_fc(fc_low = options.fecutoff)
                 others.append(peakdetect.peaks)
 
-            logging.getLogger().setLevel(logging.INFO)
+            logging.getLogger().setLevel(logging.INFO) # enable logging
             return _snapatac2.find_reproducible_peaks(merged, others, blacklist)
 
         logging.info("Calling peaks...")