two pointer sweep logic fix

taoliu · taoliu · commit e777219ce4b7 · 2026-03-11T19:00:19.000-07:00
diff --git a/.github/workflows/publish-to-gh-with-sphinx.yml b/.github/workflows/publish-to-gh-with-sphinx.yml
@@ -2,7 +2,7 @@ name: Build and Publish Sphinx Documentation
 
 on:
   push:
-    branches: ["master"]
+    branches: ["main"]
   # Allows you to run this workflow manually from the Actions tab
   workflow_dispatch:
 
diff --git a/MACS3/Signal/PairedEndTrack.py b/MACS3/Signal/PairedEndTrack.py
@@ -1832,11 +1832,18 @@ def _anndata_ncls(self, regions):
         return adata
 
     def _two_pointer_sweep(self,regions):
+        barcode_items = list(self.barcode_dict.items())
+        barcodes = [b.decode() if isinstance(b, (bytes, bytearray)) else str(b) for b, _ in barcode_items]
+        barcode_ids = np.array([i for _, i in barcode_items], dtype=np.int64)
+
+        id_map = np.full(barcode_ids.max() + 1, -1, dtype=np.int64)
+        id_map[barcode_ids] = np.arange(len(barcode_ids), dtype=np.int64)
+        n_cells = len(barcodes)
 
         # # array of locataions and counts
-        # fragment_locs = self.locations[chrom][:self.size[chrom]]
+        # fragment_locs = petrack.locations[chrom][:petrack.size[chrom]]
         # # arrray of barcodes
-        # fragment_barcodes = self.barcodes[chrom][:self.size[chrom]]
+        # fragment_barcodes = petrack.barcodes[chrom][:petrack.size[chrom]]
 
         peak_names = []
         peak_data = []
@@ -1853,9 +1860,7 @@ def _two_pointer_sweep(self,regions):
         regions.sort()
         peak_counter = 0
         for chrom in regions.regions.keys():
-            size = self.size[chrom]
-            locs = self.locations[chrom]
-            #barcodes = self.barcodes
+            #barcodes = petrack.barcodes
             regions_c = regions.regions[chrom]
             local_peak = []
             ### regions empty skip
@@ -1865,61 +1870,67 @@ def _two_pointer_sweep(self,regions):
                 peak_data.append((chrom.decode(),start,end))
                 local_peak.append(peak_counter - 1)
 
-            fragment_locs = self.locations[chrom][:self.size[chrom]]
-            fragment_barcodes = self.barcodes[chrom][:self.size[chrom]]
+            fragment_locs = self.locations[chrom]
+            fragment_barcodes = self.barcodes[chrom]
 
             if len(fragment_locs) == 0:
                 continue
 
             frag_idx = 0
             local_peak_idx = 0
-            frag_iter = iter(fragment_locs).__next__ # fragment
-            peak_iter = iter(regions_c).__next__    # peaks
             frag_len = len(fragment_locs)
             peak_len = len(regions_c)
-            frag = frag_iter() # (l,r,c)
+            frag = fragment_locs[frag_idx]
             remaining_frag_len = frag_len - 1
-            peak = peak_iter() # (start, end)
+            peak = regions_c[local_peak_idx] # (start, end)
             remaining_peak_len = peak_len - 1
 
             # inside two_pointer_sweep, replace only the while-loop block with this
-
+            back_trace_frag = 0
             while True:
                 frag_start, frag_end = frag[0], frag[1]
                 peak_start, peak_end = peak[0], peak[1]
 
-                # peak inside fragment
-                if frag_start <= peak_start and peak_end <= frag_end:
+                # peak overlap fragment
+                if frag_start <= peak_end and peak_start <= frag_end:
                     bc_id = int(fragment_barcodes[frag_idx])
                     row_id = barcode_id_to_row.get(bc_id, -1)
                     if row_id >= 0:
                         rows.append(int(row_id))
                         columns.append(local_peak[local_peak_idx])
                         data.append(int(frag[2]))
 
-                    # move to next peak (same fragment may contain multiple peaks)
-                    if remaining_peak_len > 0:
-                        peak = peak_iter()
-                        remaining_peak_len -= 1
-                        local_peak_idx += 1
+                    if frag_end > peak_end:
+                        # overhang case, perhaps the frag can still contain next peak(s)
+                        back_trace_frag += 1
+                    remaining_frag_len -= 1
+                    if remaining_frag_len >= 0:
+                        frag_idx += 1
+                        frag = fragment_locs[frag_idx]
                         continue
                     else:
                         break
 
-                # if fragment ends before peak starts -> advance fragment
-                if frag_end < peak_start:
+                # if fragment ends before peak ends -> advance fragment
+                if frag_end < peak_end:
                     remaining_frag_len -= 1
                     if remaining_frag_len >= 0:
-                        frag = frag_iter()
                         frag_idx += 1
+                        frag = fragment_locs[frag_idx]
                     else:
                         break
                 else:
-                    # peak starts before fragment ends but not contained -> advance peak
-                    if remaining_peak_len:
-                        peak = peak_iter()
-                        remaining_peak_len -= 1
+                    # advance peak
+                    remaining_peak_len -= 1
+                    if remaining_peak_len >= 0:
                         local_peak_idx += 1
+                        peak = regions_c[local_peak_idx]
+                        # we will also check if backtrace > 0 or not, if so, we need to move back the fragment pointer to check those peaks that were skipped due to overhang
+                        if back_trace_frag > 0:
+                            frag_idx -= back_trace_frag
+                            remaining_frag_len += back_trace_frag
+                            frag = fragment_locs[frag_idx]
+                        back_trace_frag = 0
                     else:
                         break
 
@@ -1928,8 +1939,11 @@ def _two_pointer_sweep(self,regions):
         obs = pd.DataFrame(index=barcodes)
         #obs['n_fragments_in_peaks'] = np.asarray(x.sum(axis=1)).ravel()
         var = pd.DataFrame(peak_data, columns=['chrom', 'start', 'end'], index=peak_names)
-        adata = ad.AnnData(X=x, obs=obs, var=var)
-        return adata
+
+        adata_peaks_loop = ad.AnnData(X=x, obs=obs, var=var)
+        return adata_peaks_loop
+
+
     
     def return_anndata(self, regions, method = "ncls"):
         """
diff --git a/README.md b/README.md
@@ -52,3 +52,13 @@ MACS3 project is sponsored by [![CZI's Essential Open Source Software for Scienc
 2008: [Model-based Analysis of ChIP-Seq
 (MACS)](https://genomebiology.biomedcentral.com/articles/10.1186/gb-2008-9-9-r137)
 
+## Note
+
+Now the default branch of MACS3 has been renamed to 'main'. If you are still using old 'master' branch in your cloned Git repository for MACS3, please use the following command to rename it:
+
+```
+git branch -m master main
+git fetch origin
+git branch -u origin/main main
+git remote set-head origin -a
+```
diff --git a/notebooks/500-PBMC-test/anndata_implemetation_dev.ipynb b/notebooks/500-PBMC-test/anndata_implemetation_dev.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 3,
    "id": "5e322d83",
    "metadata": {},
    "outputs": [],
@@ -26,7 +26,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "id": "804f83bc",
    "metadata": {},
    "outputs": [],
@@ -35,20 +35,20 @@
     "\n",
     "petrack = PETrackII()\n",
     "petrack.add_loc(b\"chr1\", 0, 100, barcode=b\"A\", count=2)   # peak_1\n",
-    "petrack.add_loc(b\"chr1\", 70, 270, barcode=b\"A\", count=1)   # peak_2\n",
-    "petrack.add_loc(b\"chr1\", 0, 100, barcode=b\"B\", count=3)   # peak_2\n",
+    "petrack.add_loc(b\"chr1\", 70, 270, barcode=b\"A\", count=1)   # peak_1&2\n",
+    "petrack.add_loc(b\"chr1\", 0, 100, barcode=b\"B\", count=3)   # peak_1\n",
     "petrack.add_loc(b\"chr1\", 175, 325, barcode=b\"C\", count=4)   # peak_2\n",
     "petrack.finalize()\n",
     "\n",
     "regions = Regions()\n",
-    "regions.add_loc(b\"chr1\", 0, 100)    # peak_1\n",
-    "regions.add_loc(b\"chr1\", 200, 300)  # peak_2\n",
-    "regions.add_loc(b\"chr1\", 500, 600)  # peak_3"
+    "regions.add_loc(b\"chr1\", 0, 100)    # peak_1 A:3 B:3 C:0\n",
+    "regions.add_loc(b\"chr1\", 200, 300)  # peak_2 A:1 B:0 C:4\n",
+    "regions.add_loc(b\"chr1\", 500, 600)  # peak_3 A:0 B:0 C:0"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 48,
    "id": "74fff8ce",
    "metadata": {},
    "outputs": [],
@@ -118,7 +118,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 49,
    "id": "744ba5ee",
    "metadata": {},
    "outputs": [
@@ -142,7 +142,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 51,
    "id": "a83976eb",
    "metadata": {},
    "outputs": [],
@@ -226,7 +226,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 52,
    "id": "b7dcf943",
    "metadata": {},
    "outputs": [
@@ -306,55 +306,63 @@
     "\n",
     "        frag_idx = 0\n",
     "        local_peak_idx = 0\n",
-    "        frag_iter = iter(fragment_locs).__next__ # fragment\n",
-    "        peak_iter = iter(regions_c).__next__    # peaks\n",
     "        frag_len = len(fragment_locs)\n",
     "        peak_len = len(regions_c)\n",
-    "        frag = frag_iter() # (l,r,c)\n",
+    "        frag = fragment_locs[frag_idx]\n",
     "        remaining_frag_len = frag_len - 1\n",
-    "        peak = peak_iter() # (start, end)\n",
+    "        peak = regions_c[local_peak_idx] # (start, end)\n",
     "        remaining_peak_len = peak_len - 1\n",
     "\n",
     "        # inside two_pointer_sweep, replace only the while-loop block with this\n",
-    "\n",
+    "        back_trace_frag = 0\n",
     "        while True:\n",
     "            frag_start, frag_end = frag[0], frag[1]\n",
     "            peak_start, peak_end = peak[0], peak[1]\n",
     "\n",
-    "            # peak inside fragment\n",
-    "            if frag_start <= peak_start and peak_end <= frag_end:\n",
+    "            # peak overlap fragment\n",
+    "            if frag_start <= peak_end and peak_start <= frag_end:\n",
     "                bc_id = int(fragment_barcodes[frag_idx])\n",
     "                row_id = barcode_id_to_row.get(bc_id, -1)\n",
     "                if row_id >= 0:\n",
     "                    rows.append(int(row_id))\n",
     "                    columns.append(local_peak[local_peak_idx])\n",
     "                    data.append(int(frag[2]))\n",
+    "                    #print(f\"peak {local_peak_idx} overlaps with fragment {frag_idx}, row_id: {row_id}, count: {frag[2]}\")\n",
     "\n",
-    "                # move to next peak (same fragment may contain multiple peaks)\n",
-    "                if remaining_peak_len > 0:\n",
-    "                    peak = peak_iter()\n",
-    "                    remaining_peak_len -= 1\n",
-    "                    local_peak_idx += 1\n",
+    "                if frag_end > peak_end:\n",
+    "                    # overhang case, perhaps the frag can still contain next peak(s)\n",
+    "                    back_trace_frag += 1\n",
+    "                remaining_frag_len -= 1\n",
+    "                if remaining_frag_len >= 0:\n",
+    "                    frag_idx += 1\n",
+    "                    frag = fragment_locs[frag_idx]\n",
     "                    continue\n",
     "                else:\n",
     "                    break\n",
     "\n",
-    "            # if fragment ends before peak starts -> advance fragment\n",
-    "            if frag_end < peak_start:\n",
+    "            # if fragment ends before peak ends -> advance fragment\n",
+    "            if frag_end < peak_end:\n",
     "                remaining_frag_len -= 1\n",
     "                if remaining_frag_len >= 0:\n",
-    "                    frag = frag_iter()\n",
     "                    frag_idx += 1\n",
+    "                    frag = fragment_locs[frag_idx]\n",
     "                else:\n",
     "                    break\n",
     "            else:\n",
-    "                # peak starts before fragment ends but not contained -> advance peak\n",
-    "                if remaining_peak_len:\n",
-    "                    peak = peak_iter()\n",
-    "                    remaining_peak_len -= 1\n",
+    "                # advance peak\n",
+    "                remaining_peak_len -= 1\n",
+    "                if remaining_peak_len >= 0:\n",
     "                    local_peak_idx += 1\n",
+    "                    peak = regions_c[local_peak_idx]\n",
+    "                    # we will also check if backtrace > 0 or not, if so, we need to move back the fragment pointer to check those peaks that were skipped due to overhang\n",
+    "                    if back_trace_frag > 0:\n",
+    "                        frag_idx -= back_trace_frag\n",
+    "                        remaining_frag_len += back_trace_frag\n",
+    "                        frag = fragment_locs[frag_idx]\n",
+    "                    back_trace_frag = 0\n",
     "                else:\n",
     "                    break\n",
+    "        #print(f\"remaining fragments: {remaining_frag_len}, remaining peaks: {remaining_peak_len}\")\n",
     "\n",
     "    n_peaks = peak_counter\n",
     "    x = sparse.coo_matrix((data,(rows,columns)),shape=(n_cells,n_peaks)).tocsr()\n",
@@ -368,19 +376,25 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 37,
+   "execution_count": 54,
    "id": "2e4bd05d",
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
+      "peak 0 overlaps with fragment 0, row_id: 0, count: 2\n",
+      "peak 0 overlaps with fragment 1, row_id: 1, count: 3\n",
+      "peak 0 overlaps with fragment 2, row_id: 0, count: 1\n",
+      "peak 1 overlaps with fragment 2, row_id: 0, count: 1\n",
+      "peak 1 overlaps with fragment 3, row_id: 2, count: 4\n",
+      "remaining fragments: -1, remaining peaks: 1\n",
       "AnnData object with n_obs × n_vars = 3 × 3\n",
       "    var: 'chrom', 'start', 'end'\n",
-      "[[2 0 0]\n",
-      " [0 0 0]\n",
-      " [0 0 0]]\n"
+      "[[3 1 0]\n",
+      " [3 0 0]\n",
+      " [0 4 0]]\n"
      ]
     }
    ],
@@ -427,7 +441,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": ".venv (3.9.6)",
+   "display_name": ".venv (3.12.8)",
    "language": "python",
    "name": "python3"
   },
@@ -441,7 +455,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.6"
+   "version": "3.12.8"
   }
  },
  "nbformat": 4,
diff --git a/notebooks/500-PBMC-test/macs3-scATAC-pbmc-500.ipynb b/notebooks/500-PBMC-test/macs3-scATAC-pbmc-500.ipynb
