increase max count in frag file to 65535, give warning if count is over the limit while parsing the file

Author: taoliu

MACS3/IO/Parser.py

@@ -1,7 +1,7 @@
 # cython: language_level=3
 # cython: profile=True
 # cython: linetrace=True
-# Time-stamp: <2025-02-15 08:24:53 Tao Liu>
+# Time-stamp: <2025-04-11 13:41:46 Tao Liu>
 
 """Module for all MACS Parser classes for input. Please note that the
 parsers are for reading the alignment files ONLY.
@@ -34,7 +34,7 @@
 logger = logging.getLogger(__name__)
 debug = logger.debug
 info = logger.info
-
+warn = logger.warn
 # ------------------------------------
 # Other modules
 # ------------------------------------
@@ -1513,17 +1513,23 @@ def pe_parse_line(self, thisline: bytes):
 
         """
         thisfields: list
+        thiscount: cython.ushort
 
         thisline = thisline.rstrip()
 
         # still only support tabular as delimiter.
         thisfields = thisline.split(b'\t')
         try:
+            try:
+                thiscount = atoi(thisfields[4])
+            except OverflowError:
+                thiscount = 65535
+                warn(f"The count in this line is over 65535, and will be capped at 65535: {thisline}")
             return (thisfields[0],
                     atoi(thisfields[1]),
                     atoi(thisfields[2]),
                     thisfields[3],
-                    atoi(thisfields[4]))
+                    thiscount)
         except IndexError:
             raise Exception("Less than 5 columns found at this line: %s\n" %
                             thisline)
@@ -1537,7 +1543,7 @@ def build_petrack(self, max_count=0):
         left_pos: cython.int
         right_pos: cython.int
         barcode: bytes
-        count: cython.uchar
+        count: cython.ushort
         i: cython.long = 0          # number of fragments
         m: cython.long = 0          # sum of fragment lengths
         tmp: bytes = b""
@@ -1591,7 +1597,7 @@ def append_petrack(self, petrack, max_count=0):
         left_pos: cython.int
         right_pos: cython.int
         barcode: bytes
-        count: cython.uchar
+        count: cython.ushort
         i: cython.long = 0          # number of fragments
         m: cython.long = 0          # sum of fragment lengths
         tmp: bytes = b""

MACS3/Signal/PairedEndTrack.py

@@ -1,6 +1,6 @@
 # cython: language_level=3
 # cython: profile=True
-# Time-stamp: <2025-02-15 08:49:30 Tao Liu>
+# Time-stamp: <2025-04-11 13:24:36 Tao Liu>
 
 """Module for filter duplicate tags from paired-end data
 
@@ -779,7 +779,7 @@ class PETrackII:
 
     def __init__(self, anno: str = "", buffer_size: cython.long = 100000):
         # dictionary with chrname as key, nparray with
-        # [('l','i4'),('r','i4'),('c','u1')] as value
+        # [('l','i4'),('r','i4'),('c','u2')] as value
         self.locations = {}
         # dictionary with chrname as key, size of the above nparray as value
         # size is to remember the size of the fragments added to this chromosome
@@ -805,7 +805,7 @@ def add_loc(self,
                 start: cython.int,
                 end: cython.int,
                 barcode: bytes,
-                count: cython.uchar):
+                count: cython.ushort):
         """Add a location to the list according to the sequence name.
 
         chromosome: mostly the chromosome name
@@ -828,7 +828,7 @@ def add_loc(self,
             # note: ['l'] is the leftmost end, ['r'] is the rightmost end of fragment.
             # ['c'] is the count number of this fragment
             self.locations[chromosome] = np.zeros(shape=self.buffer_size,
-                                                  dtype=[('l', 'i4'), ('r', 'i4'), ('c', 'u1')])
+                                                  dtype=[('l', 'i4'), ('r', 'i4'), ('c', 'u2')])
             self.barcodes[chromosome] = np.zeros(shape=self.buffer_size,
                                                  dtype='i4')
             self.locations[chromosome][0] = (start, end, count)

MACS3/Signal/PileupV2.py

@@ -1,6 +1,6 @@
 # cython: language_level=3
 # cython: profile=True
-# Time-stamp: <2025-02-14 11:56:44 Tao Liu>
+# Time-stamp: <2025-04-11 13:26:28 Tao Liu>
 
 """Module Description:
 
@@ -112,7 +112,7 @@ def make_PV_from_LRC(LRC_array: cnp.ndarray,
     `mapping_func( L, R )` or simply 1 if mapping_func is the default.
 
     LRC array is an np.ndarray as with dtype
-    [('l','i4'),('r','i4'),('c','u1')] with length of N
+    [('l','i4'),('r','i4'),('c','u2')] with length of N
 
     PV array is an np.ndarray with
     dtype=[('p','i4'),('v','f4')] with length of 2N
@@ -122,7 +122,7 @@ def make_PV_from_LRC(LRC_array: cnp.ndarray,
     i: cython.ulong
     L: cython.int
     R: cython.int
-    C: cython.uchar
+    C: cython.ushort
     PV: cnp.ndarray
 
     l_LRC = LRC_array.shape[0]
@@ -292,7 +292,7 @@ def pileup_from_LRC(LRC_array: cnp.ndarray,
 
     User needs to provide a numpy array of left and right positions
     and the counts, with
-    dtype=[('l','i4'),('r','i4'),('c','u1')]. User also needs to
+    dtype=[('l','i4'),('r','i4'),('c','u2')]. User also needs to
     provide a mapping function to map the left and right position to
     certain weight.
 

docs/source/docs/bdgdiff.md

@@ -30,14 +30,16 @@ condition B while comparing treatment B over control B, and also the
 likelihood test while comparing treatment A and treatment B can't
 decide which condition is stronger.
 
+
+
 The likelihood function we used while comparing two conditions: ChIP
 (enrichment) or control (chromatin bias) is:
 
 $$ln(LR) = x*(ln(x)-ln(y)) + y - x$$
 
 Here $LR$ is the likelihood ratio, x is the signal (fragment pileup)
 we observed in condition 1, and y is the signal in condition
-2. And $ln$ is the natural logarithm.
+2. And $ln$ is the natural logarithm. 
 
 Note: All regions on the same chromosome in the bedGraph file should
 be continuous so only bedGraph files from MACS3 are acceptable.

docs/source/docs/callpeak.md

@@ -8,6 +8,24 @@ invoked by `macs3 callpeak` . If you type this command with `-h`, you
 will see a full description of command-line options. Here we only list
 the essentials.
 
+## Examples of Commandline Usage
+1. Peak calling for regular TF ChIP-seq:
+    $ macs3 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01
+2. Broad peak calling on Histone Mark ChIP-seq:
+    $ macs3 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1
+3. Peak calling on ATAC-seq (paired-end mode):
+    $ macs3 callpeak -f BAMPE -t ATAC.bam -g hs -n test -B -q 0.01
+4. Peak calling on ATAC-seq (focusing on insertion sites, and using single-end mode):
+    $ macs3 callpeak -f BAM -t ATAC.bam -g hs -n test -B -q 0.01 --shift -50 --extension 100
+5. Peak calling on scATAC-seq (paired-end mode):
+    $ macs3 callpeak -f BEDPE -t scATAC.fragments.tsv.gz -g hs -n test -B -q 0.01 -n test
+6. Peak calling on scATAC-seq (paired-end mode):
+    $ macs3 callpeak -f FRAG -t scATAC.fragments.tsv.gz -g hs -n test -B -q 0.01 -n test
+7. Peak calling on scATAC-seq (paired-end mode) and only for given barcodes:
+    $ macs3 callpeak -f FRAG -t scATAC.fragments.tsv.gz -g hs -n test -B -q 0.01 -n test --barcodes barcodes.txt
+
+
+
 ## Essential Commandline Options
 
 ### Input and Output

test/test_PairedEndTrack.py

@@ -1,5 +1,5 @@
 #!/usr/bin/env python
-# Time-stamp: <2025-02-12 13:44:51 Tao Liu>
+# Time-stamp: <2025-04-11 13:42:39 Tao Liu>
 
 import unittest
 from MACS3.Signal.PairedEndTrack import PETrackI, PETrackII
@@ -96,7 +96,6 @@ def setUp(self):
                               (b"chrY", 50, 160, b"0w#AAACGAACAAGTAAGA", 2),  # will be excluded
                               (b"chrY", 100, 170, b"0w#AAACGAACAAGTAAGA", 3)
                               ]
-
         self.exclude_regions = [(b"chrY", 10, 75)]
         self.x_regions = Regions()
         for r in self.exclude_regions:

test/test_Parser.py

@@ -1,5 +1,5 @@
 #!/usr/bin/env python
-# Time-stamp: <2025-02-15 08:23:55 Tao Liu>
+# Time-stamp: <2025-04-11 13:46:05 Tao Liu>
 
 import unittest
 
@@ -40,15 +40,9 @@ def setUp(self):
         self.bedpefile = "test/tiny.bedpe.gz"
         self.samfile = "test/tiny.sam.gz"
         self.bamfile = "test/tiny.bam"
-        self.fragfile = "test/test.fragments.tsv.gz"
-
+        self.fragfile = "test/tiny.frag.tsv.gz"
+        
     def test_fragment_file(self):
         p = FragParser(self.fragfile)
         petrack = p.build_petrack()
         petrack.finalize()
-        bdg = petrack.pileup_bdg()
-        bdg2 = petrack.pileup_bdg2()
-        peaks = bdg.call_peaks(cutoff=10, min_length=200, max_gap=100)
-        peaks2 = bdg2.call_peaks(cutoff=10, min_length=200, max_gap=100)
-        # print(peaks)
-        # print(peaks2)