Problem class: python-20261012-pangenome-minimizer-chaining-with-structural-variants
All 10 tests pass and SOLUTION.md is regenerated with the cleaned imports. Here is the complete, verified solution.
Problem class: python-20261012-pangenome-minimizer-chaining-with-structural-variants
Environment: Python 3.11+ (verified on 3.14) · Dependencies: standard library only
Artifacts written to ~/: pangenome_chain.py, test_pangenome_chain.py, SOLUTION.md.
Given millions of noisy long reads, a colored variation graph, and a memory cap:
A "solve the problem" answer is therefore a small, verifiable system, not a one-line patch.
| # | Naive approach | Root cause of failure | Symptom |
|---|---|---|---|
| RC1 | Compare every anchor pair (O(n^2)) |
For millions of anchors ~10^12 pairs, and the DP predecessor matrix is O(n^2) memory |
Timeout / OOM |
| RC2 | Greedily extend over read coordinate only | A repetitive read minimizer maps to several graph positions; a read base can be reused | Chimeric / inflated alignments |
| RC3 | Linearise by node id / DFS and 1-D chain | Parallel bubble branches are alternatives; nothing forbids A then B | Path no sample carries |
| RC4 | Ignore orientation | Half of all mappings vanish | Reads lost / wrong loci |
| RC5 | Python set of sample ids per anchor/chain |
Millions of sets blow the cap; provenance lost on merge | OOM / provenance dropped |
| RC6 | Buffer all anchors for all reads | Memory grows with corpus, not graph | Cap violated |
The minimizer-chain objective is
score(chain) = sum_i w_i - lambda * sum over j->i of ( (r_i-r_j) + (g_i-g_j) )
where r is read coordinate and g graph coordinate. For current anchor i:
max_j [ f[j] - lambda*((r_i-r_j)+(g_i-g_j)) ]
= max_j [ f[j] + lambda*(r_j+g_j) ] - lambda*(r_i+g_i)
= max_j [ V_j ] - lambda*(r_i+g_i)
The term -lambda*(r_i+g_i) is identical for every predecessor, so the DP becomes range-max queries of V_j = f[j]+lambda*(r_j+g_j) over graph coordinates. A segment tree gives O(n log n) time and O(n) memory.
Consecutive anchors are feasible only if r_i - r_j >= k (k = minimizer length). This is exact, enforced with delayed two-pointer insertion: j enters the tree only once r_i >= r_j + k. Every read base is consumed at most once.
Bubble branches get (bubble_id, branch_id) tags; a transition inside the same bubble that changes branch is forbidden. Bubbles are bounded, so the conflicting set is tiny; per-branch buckets give exact same-branch predecessors, and the segment-tree query is split around the bubble span to exclude other branches.
Minimizers are canonical min(kmer, revcomp(kmer)) with a strand bit. Read maps forward iff s_r*s_g == +1; otherwise the read coordinate is mirrored r = L-k-p, restoring monotonicity in both r and g. Anchors chain per orientation, best wins.
Colour sets are big-int bitmasks (one bit per sample); a chain's provenance is the AND of its anchors' masks — exactly the samples containing the whole path.
stream_reads consumes one read at a time; only the graph index + current read anchors are resident. Minimizers occurring more than max_occ times are dropped.
reads --minimizers--> canonical k-mers --graph index--> anchors (r, g, strand, colors)
|
split by strand -------+
v
graph linearisation: gidx = topo prefix sum
bubble detection: (bubble, branch) tags
|
v
DP over read order: f[i] = w_i + max(0, max_j V_j - lambda*(r_i+g_i))
V_j = f[j] + lambda*(r_j+g_j) [segment tree, range max]
delayed insertion enforces r_i - r_j >= k
branch exclusion + per-branch buckets for bubbles
max_gap bounds an SV jump
|
v
chain reconstruction + provenance = AND of masks
pangenome_chain.py#!/usr/bin/env python3
"""
pangenome_chain.py
==================
Reference implementation for:
Maximum-scoring chaining of minimizer anchors over a *colored variation
graph*, under a memory cap, tolerating bounded structural-variant (SV)
bubbles, preserving sample provenance, and handling reverse-complement
read orientations.
Design summary
--------------
* Minimizers are canonical k-mers (min(kmer, revcomp(kmer))); the strand bit
records orientation. A monotonic deque computes them in O(len(seq)).
* A colored variation graph is linearised with a topological prefix sum so
that every directed path has strictly increasing graph coordinates
``gidx``. Parallel bubble branches therefore occupy disjoint coordinate
ranges.
* Anchors are points (r, g) with r = read coordinate and g = graph coordinate.
A chain is an increasing subsequence in both coordinates.
* Score model: chain score = sum(anchor weights) - lambda * Manhattan gap,
where the gap between consecutive anchors j -> i is
(r_i - r_j) + (g_i - g_j).
* The DP is O(n log n) because the predecessor term decomposes:
f[i] = w_i - lambda*(r_i+g_i) + max_j (f[j] + lambda*(r_j+g_j))
and the max is a prefix/range maximum over graph coordinates, maintained in
a segment tree.
* Pruning (the sub-quadratic guarantee): predecessor j is *dominated* by k if
r_j <= r_k, g_j <= g_k and V_j <= V_k, where V = f + lambda*(r+g). Only the
per-coordinate maximum V is retained, so a segment tree with n leaves answers
every query and the total work is O(n log n). See `proof_domination()`.
* Double counting is prevented by the read-overlap constraint
``r_i - r_j >= k``: a read base is used by at most one anchor. It is
implemented with delayed insertion (a two-pointer sweep), so it stays exact
and does not change the asymptotic cost.
* Bubbles: anchors carry (bubble_id, branch_id). A transition that stays
inside the same bubble but changes branch is forbidden. The segment-tree
query is split around the bubble's coordinate range and a per-branch maximum
supplies same-branch predecessors.
* Provenance: node colour sets are represented as Python big-int bitmasks.
The chain provenance is the AND (intersection) of its anchors' masks, which
is exactly the set of samples whose graph contains the whole chain path.
Author: reference solution
"""
from __future__ import annotations
import zlib
from bisect import bisect_left, bisect_right
from collections import defaultdict, deque
from dataclasses import dataclass
from typing import Dict, Iterable, List, Optional, Sequence, Tuple
NEG = float("-inf")
# ---------------------------------------------------------------------------
# DNA helpers
# ---------------------------------------------------------------------------
_COMP = str.maketrans("ACGTNacgtn", "TGCANtgcan")
def revcomp(s: str) -> str:
"""Reverse complement of a DNA string."""
return s.translate(_COMP)[::-1]
def canonical_kmer(kmer: str) -> Tuple[str, int]:
"""Return (canonical kmer, strand).
strand == +1 iff kmer is equal to its canonical form.
strand == -1 iff the reverse complement is the canonical form.
Canonical form is invariant under reverse complementation.
"""
rc = revcomp(kmer)
if kmer <= rc:
return kmer, 1
return rc, -1
def hash_kmer(kmer: str) -> int:
"""Stable, fast, process-independent hash (CRC32 of the ASCII bytes)."""
return zlib.crc32(kmer.encode("ascii"))
# ---------------------------------------------------------------------------
# Minimizers (O(n) sliding-window minimum)
# ---------------------------------------------------------------------------
def minimizers(seq: str, k: int, w: int, hash_fn=hash_kmer):
"""Return minimizers as a list of (pos, canonical_kmer, strand, h).
The window contains the hashes of ``w`` consecutive k-mers. Ties are
broken by the left-most position, matching the usual minimizer definition
min by (hash, position).
"""
n = len(seq)
if k <= 0 or w <= 0 or n < k:
return []
m = n - k + 1
hs = [0] * m
strands = [0] * m
cans: List[str] = [""] * m
for i in range(m):
c, s = canonical_kmer(seq[i : i + k])
cans[i] = c
strands[i] = s
hs[i] = hash_fn(c)
out: List[Tuple[int, str, int, int]] = []
dq: deque[int] = deque() # indices, hashes increasing from front
for i in range(m):
# Keep the minimum by (hash, position): pop strictly larger hashes only,
# so an earlier equal-hash element wins.
while dq and hs[dq[-1]] > hs[i]:
dq.pop()
dq.append(i)
while dq[0] <= i - w:
dq.popleft()
if i >= w - 1:
j = dq[0]
if not out or out[-1][0] != j:
out.append((j, cans[j], strands[j], hs[j]))
return out
def minimizers_naive(seq: str, k: int, w: int, hash_fn=hash_kmer):
"""Reference O(n*w) minimizer implementation, used only for tests."""
n = len(seq)
if n < k:
return []
m = n - k + 1
recs = []
for i in range(m):
c, s = canonical_kmer(seq[i : i + k])
recs.append((hash_fn(c), i, c, s))
out = []
for start in range(m - w + 1):
best = min(recs[start : start + w], key=lambda x: (x[0], x[1]))
if not out or out[-1][0] != best[1]:
out.append((best[1], best[2], best[3], best[0]))
return out
# ---------------------------------------------------------------------------
# Colored variation graph
# ---------------------------------------------------------------------------
@dataclass
class Node:
seq: str
colors: int = 0 # bitmask of sample ids
class ColoredVariationGraph:
"""A small, explicit colored variation graph.
Parameters
----------
nodes : dict[int, Node]
edges : iterable of (u, v) directed edges
"""
def __init__(self, nodes: Dict[int, Node], edges: Iterable[Tuple[int, int]]):
self.nodes = dict(nodes)
self.edges = list(dict.fromkeys(edges)) # dedup, keep order
self.adj: Dict[int, List[int]] = defaultdict(list)
for u, v in self.edges:
self.adj[u].append(v)
self.gidx: Dict[int, int] = {}
self.glen: Dict[int, int] = {}
# -- linearisation -----------------------------------------------------
def topological_order(self) -> List[int]:
indeg = {v: 0 for v in self.nodes}
for u, v in self.edges:
indeg.setdefault(v, 0)
indeg[v] += 1
dq = deque(sorted(v for v in self.nodes if indeg[v] == 0))
order: List[int] = []
while dq:
u = dq.popleft()
order.append(u)
for v in sorted(self.adj[u]):
indeg[v] -= 1
if indeg[v] == 0:
dq.append(v)
if len(order) != len(self.nodes):
raise ValueError("variation graph is not a DAG; unfold cycles first")
return order
def linearize(self) -> List[int]:
"""Assign strictly increasing coordinates along every path.
``gidx[v]`` is the cumulative sequence length of preceding nodes in
topological order. Because parallel branches occupy disjoint ranges,
no directed path can revisit a coordinate, and any path is strictly
increasing in ``gidx``.
"""
order = self.topological_order()
pos = 0
for v in order:
self.gidx[v] = pos
self.glen[v] = len(self.nodes[v].seq)
pos += len(self.nodes[v].seq)
return order
# -- minimizer index ---------------------------------------------------
def index_minimizers(self, k: int, w: int, max_occ: int = 1000, hash_fn=hash_kmer):
"""Build canonical-kmer -> [(node, offset, strand)] index.
Keys occurring in more than ``max_occ`` positions are dropped: this is
the explicit memory cap for repetitive sequence and prevents a single
repeat from exploding into millions of anchors.
"""
counts: Dict[str, int] = defaultdict(int)
per_node: Dict[int, list] = {}
for nid, node in self.nodes.items():
ms = minimizers(node.seq, k, w, hash_fn)
per_node[nid] = ms
for _, can, _, _ in ms:
counts[can] += 1
idx: Dict[str, List[Tuple[int, int, int]]] = defaultdict(list)
dropped = 0
for nid, ms in per_node.items():
for pos, can, strand, _ in ms:
if counts[can] <= max_occ:
idx[can].append((nid, pos, strand))
else:
dropped += 1
return idx, dropped
# -- bubble detection (2-branch diamonds) ------------------------------
def detect_bubbles(self, max_nodes: int = 10000) -> Tuple[dict, dict]:
"""Detect simple 2-branch diamond bubbles.
A diamond is u -> v -> w and u -> x -> w. Node v is assigned
branch 0 and x branch 1. Returns
bubbles : bubble_id -> (lo, hi) coordinate span
node_bubble : node_id -> (bubble_id, branch_id)
"""
bubbles: Dict[int, Tuple[int, int]] = {}
node_bubble: Dict[int, Tuple[int, int]] = {}
bid = 0
for u in self.topological_order():
succ = sorted(set(self.adj[u]))
if len(succ) != 2:
continue
v, x = succ
joins = sorted(set(self.adj[v]) & set(self.adj[x]))
if not joins:
continue
w = joins[0]
members = {v: 0, x: 1}
for m, b in members.items():
node_bubble[m] = (bid, b)
lo = min(self.gidx[m] for m in members)
hi = max(self.gidx[m] + self.glen[m] for m in members)
bubbles[bid] = (lo, hi)
bid += 1
return bubbles, node_bubble
# ---------------------------------------------------------------------------
# Anchors
# ---------------------------------------------------------------------------
@dataclass(slots=True)
class Anchor:
read_id: int
ridx: int # read coordinate along the read's forward orientation
gidx: int # graph coordinate (topological prefix + in-node offset)
strand: int # +1 / -1 orientation of the read relative to the graph
weight: float # anchor weight
colors: int # sample provenance bitmask
node: int
pos: int # original read start position
bubble: int = -1
branch: int = -1
@property
def key(self):
return (self.read_id, self.strand, self.ridx, self.gidx)
def build_anchors(
read_id: int,
seq: str,
graph: ColoredVariationGraph,
index: Dict[str, List[Tuple[int, int, int]]],
k: int,
w: int,
node_bubble: Optional[Dict[int, Tuple[int, int]]] = None,
hash_fn=hash_kmer,
) -> List[Anchor]:
"""Project read minimizers onto the graph, producing anchors.
Reverse-complement handling
---------------------------
Let the read minimizer have strand ``s_r`` and the graph minimizer strand
``s_g`` relative to their common canonical k-mer. The read maps to the
graph forward direction iff ``s_r * s_g == +1``. For a reverse mapping we
mirror the read coordinate (``r = L - k - p``) so that a chain is still
strictly increasing in both r and g. The graph coordinate is always the
forward in-node coordinate ``gidx[node] + offset``.
"""
L = len(seq)
anchors: List[Anchor] = []
if node_bubble is None:
node_bubble = {}
for p, can, s_r, _h in minimizers(seq, k, w, hash_fn):
hits = index.get(can)
if not hits:
continue
for nid, npos, s_g in hits:
astrand = s_r * s_g
node = graph.nodes[nid]
gx = graph.gidx[nid] + npos
if astrand == 1:
rx = p
else:
rx = L - k - p
bb = node_bubble.get(nid, (-1, -1))
anchors.append(
Anchor(
read_id=read_id,
ridx=rx,
gidx=gx,
strand=astrand,
weight=1.0,
colors=node.colors,
node=nid,
pos=p,
bubble=bb[0],
branch=bb[1],
)
)
return anchors
# ---------------------------------------------------------------------------
# Segment tree: range maximum of V = f + lambda*(r + g)
# ---------------------------------------------------------------------------
class SegMax:
"""Segment tree over compressed graph coordinates.
Each node stores the pair (V, index) that is maximal under
``_better``: larger V wins; on an exact tie the smaller anchor index
wins. The deterministic tie-break makes the fast DP and the brute-force
reference select the identical optimal chain.
"""
__slots__ = ("n", "val", "idx")
def __init__(self, n: int):
self.n = 1
while self.n < max(1, n):
self.n <<= 1
self.val = [NEG] * (2 * self.n)
self.idx = [-1] * (2 * self.n)
@staticmethod
def _better(v, i, bv, bi) -> bool:
return v > bv or (v == bv and (bi < 0 or (i >= 0 and i < bi)))
def update(self, i: int, v: float, j: int) -> None:
p = i + self.n
if not self._better(v, j, self.val[p], self.idx[p]):
return
self.val[p] = v
self.idx[p] = j
p >>= 1
while p:
l, r = 2 * p, 2 * p + 1
if self._better(self.val[l], self.idx[l], self.val[r], self.idx[r]):
self.val[p], self.idx[p] = self.val[l], self.idx[l]
else:
self.val[p], self.idx[p] = self.val[r], self.idx[r]
p >>= 1
def query(self, lo: int, hi: int) -> Tuple[float, int]:
"""Max over inclusive compressed indices [lo, hi]."""
if lo > hi:
return (NEG, -1)
lo += self.n
hi += self.n
best_v, best_i = NEG, -1
while lo <= hi:
if lo & 1:
if self._better(self.val[lo], self.idx[lo], best_v, best_i):
best_v, best_i = self.val[lo], self.idx[lo]
lo += 1
if not (hi & 1):
if self._better(self.val[hi], self.idx[hi], best_v, best_i):
best_v, best_i = self.val[hi], self.idx[hi]
hi -= 1
lo >>= 1
hi >>= 1
return best_v, best_i
# ---------------------------------------------------------------------------
# The chaining DP
# ---------------------------------------------------------------------------
@dataclass
class ChainResult:
chain: List[Anchor]
score: float
colors: int
f: List[float]
def _coordinate_compress(anchors: Sequence[Anchor]) -> Tuple[List[int], Dict[int, int]]:
gs = sorted({a.gidx for a in anchors})
return gs, {g: i for i, g in enumerate(gs)}
def chain_one_orientation(
anchors: Sequence[Anchor],
k: int,
lam: float,
max_gap: int,
bubbles: Dict[int, Tuple[int, int]],
) -> ChainResult:
"""O(n log n) maximum-score chain for anchors of a single read+strand.
``max_gap`` bounds the allowed graph jump across an SV bubble. The DP is
exact for the Manhattan-gap objective; the only pruning performed is the
dominance pruning embedded in `SegMax.update` (a leaf keeps the single best
V seen at that graph coordinate) plus the delayed read-overlap insertion,
which is an exact feasibility constraint rather than a heuristic.
"""
n = len(anchors)
if n == 0:
return ChainResult([], 0.0, 0, [])
order = sorted(range(n), key=lambda i: (anchors[i].ridx, anchors[i].gidx))
a = [anchors[i] for i in order]
gs, gpos = _coordinate_compress(a)
tree = SegMax(len(gs))
f = [0.0] * n
prev = [-1] * n
col = [0] * n
# Per-(bubble,branch) inserted predecessors. Bubbles are bounded, so the
# number of anchors per branch is small and an exact range-filtered scan is
# cheap. (A per-branch Fenwick tree would give the same result with a
# larger constant; the scan is used for clarity.)
branch_buckets: Dict[Tuple[int, int], List[int]] = defaultdict(list)
Varr = [NEG] * n
inserted = 0
for i in range(n):
ai = a[i]
# --- delayed insertion: only non-overlapping read predecessors -----
while inserted < i and a[inserted].ridx + k <= ai.ridx:
aj = a[inserted]
V = f[inserted] + lam * (aj.ridx + aj.gidx)
Varr[inserted] = V
tree.update(gpos[aj.gidx], V, inserted)
if aj.bubble >= 0:
branch_buckets[(aj.bubble, aj.branch)].append(inserted)
inserted += 1
# --- predecessor query --------------------------------------------
g_hi = ai.gidx - 1
g_lo = ai.gidx - max_gap
best: Tuple[float, int] = (NEG, -1)
if ai.bubble >= 0:
blo, bhi = bubbles[ai.bubble]
# Outside the bubble on the left.
hi_idx = bisect_right(gs, min(g_hi, blo - 1)) - 1
lo_idx = bisect_left(gs, g_lo)
if lo_idx <= hi_idx:
cand = tree.query(lo_idx, hi_idx)
if SegMax._better(cand[0], cand[1], best[0], best[1]):
best = cand
# Outside the bubble on the right (kept for completeness).
lo2 = bisect_left(gs, max(g_lo, bhi + 1))
hi2 = bisect_right(gs, g_hi) - 1
if lo2 <= hi2:
cand = tree.query(lo2, hi2)
if SegMax._better(cand[0], cand[1], best[0], best[1]):
best = cand
# Same branch inside the bubble, restricted to gidx < gidx_i.
for j in branch_buckets.get((ai.bubble, ai.branch), ()):
if (
a[j].gidx < ai.gidx
and ai.gidx - a[j].gidx <= max_gap
and SegMax._better(Varr[j], j, best[0], best[1])
):
best = (Varr[j], j)
else:
lo_idx = bisect_left(gs, g_lo)
hi_idx = bisect_right(gs, g_hi) - 1
if lo_idx <= hi_idx:
best = tree.query(lo_idx, hi_idx)
base = best[0] - lam * (ai.ridx + ai.gidx)
if best[1] >= 0 and base > 0.0:
f[i] = ai.weight + base
prev[i] = best[1]
col[i] = ai.colors & col[best[1]]
else:
f[i] = ai.weight
prev[i] = -1
col[i] = ai.colors
end = max(range(n), key=lambda i: f[i])
chain: List[Anchor] = []
cur = end
while cur >= 0:
chain.append(a[cur])
cur = prev[cur]
chain.reverse()
return ChainResult(chain, f[end], col[end], f)
def chain_read(
read_id: int,
seq: str,
graph: ColoredVariationGraph,
index: Dict[str, List[Tuple[int, int, int]]],
k: int,
w: int,
lam: float = 0.01,
max_gap: int = 10000,
node_bubble: Optional[Dict[int, Tuple[int, int]]] = None,
bubbles: Optional[Dict[int, Tuple[int, int]]] = None,
hash_fn=hash_kmer,
) -> ChainResult:
"""Full pipeline for one read: build anchors, chain each strand, take best."""
if node_bubble is None:
node_bubble = {}
if bubbles is None:
bubbles = {}
anchors = build_anchors(
read_id, seq, graph, index, k, w, node_bubble=node_bubble, hash_fn=hash_fn
)
best = ChainResult([], 0.0, 0, [])
for strand in (+1, -1):
sub = [a for a in anchors if a.strand == strand]
if not sub:
continue
res = chain_one_orientation(sub, k, lam, max_gap, bubbles)
if res.score > best.score:
best = res
return best
# ---------------------------------------------------------------------------
# Brute force reference (O(n^2)) for verification
# ---------------------------------------------------------------------------
def chain_brute(
anchors: Sequence[Anchor],
k: int,
lam: float,
max_gap: int,
bubbles: Dict[int, Tuple[int, int]],
) -> ChainResult:
n = len(anchors)
if n == 0:
return ChainResult([], 0.0, 0, [])
order = sorted(range(n), key=lambda i: (anchors[i].ridx, anchors[i].gidx))
a = [anchors[i] for i in order]
f = [0.0] * n
prev = [-1] * n
col = [0] * n
for i in range(n):
ai = a[i]
f[i] = ai.weight
prev[i] = -1
col[i] = ai.colors
for j in range(i):
aj = a[j]
if aj.ridx + k > ai.ridx:
continue
if aj.gidx >= ai.gidx:
continue
if ai.gidx - aj.gidx > max_gap:
continue
if aj.bubble >= 0 and aj.bubble == ai.bubble and aj.branch != ai.branch:
continue
Vj = f[j] + lam * (aj.ridx + aj.gidx)
base = Vj - lam * (ai.ridx + ai.gidx)
if base <= 0.0:
continue
val = ai.weight + base
if val > f[i]:
f[i] = val
prev[i] = j
col[i] = ai.colors & col[j]
end = max(range(n), key=lambda i: f[i])
chain: List[Anchor] = []
cur = end
while cur >= 0:
chain.append(a[cur])
cur = prev[cur]
chain.reverse()
return ChainResult(chain, f[end], col[end], f)
# ---------------------------------------------------------------------------
# Sub-quadratic correctness argument
# ---------------------------------------------------------------------------
def proof_domination() -> str:
return (
"Dominance pruning is exact for the Manhattan-gap DP.\n"
"For a fixed future anchor i the contribution of a predecessor j is\n"
" f[j] - lam*((r_i-r_j)+(g_i-g_j))\n"
" = [f[j] + lam*(r_j+g_j)] - lam*(r_i+g_i)\n"
" = V_j - lam*(r_i+g_i).\n"
"The term -lam*(r_i+g_i) is identical for every predecessor, so a\n"
"predecessor is fully characterised by V_j. If two anchors j,k satisfy\n"
"r_j <= r_k and g_j <= g_k then whenever V_j <= V_k, j can never beat k\n"
"for any future i (future i has r_i >= r_k, g_i >= g_k and the gap term\n"
"is monotone). Hence only the per-graph-coordinate maximum V is needed.\n"
"The segment tree stores exactly those maxima, so it answers the true\n"
"maximum over all feasible predecessors and runs in O(n log n)."
)
# ---------------------------------------------------------------------------
# Memory-capped streaming driver
# ---------------------------------------------------------------------------
def stream_reads(
reads: Iterable[str],
graph: ColoredVariationGraph,
k: int = 15,
w: int = 10,
lam: float = 0.01,
max_gap: int = 10000,
max_occ: int = 1000,
hash_fn=hash_kmer,
):
"""Process reads one at a time under a memory cap.
Only one read's anchor list plus the graph minimizer index is resident.
Retaining the final chains for millions of reads would itself exceed the
cap, so callers should consume the generator and write results to disk.
"""
graph.linearize()
index, dropped = graph.index_minimizers(k, w, max_occ=max_occ, hash_fn=hash_fn)
bubbles, node_bubble = graph.detect_bubbles()
for rid, seq in enumerate(reads):
res = chain_read(
rid,
seq,
graph,
index,
k,
w,
lam=lam,
max_gap=max_gap,
node_bubble=node_bubble,
bubbles=bubbles,
hash_fn=hash_fn,
)
yield rid, res
__all__ = [
"revcomp",
"canonical_kmer",
"minimizers",
"minimizers_naive",
"Node",
"ColoredVariationGraph",
"Anchor",
"build_anchors",
"SegMax",
"ChainResult",
"chain_one_orientation",
"chain_read",
"chain_brute",
"proof_domination",
"stream_reads",
]
Claim. For the Manhattan-gap DP, keeping only the per-graph-coordinate maximum
V_j = f[j] + lambda*(r_j+g_j)is exact.
Proof. For a fixed future anchor i, predecessor j contributes V_j - lambda*(r_i+g_i). The subtracted term is independent of j, so the best predecessor maximises V_j. If r_j <= r_k and g_j <= g_k and V_j <= V_k, then every future i usable by j also has r_i >= r_k, g_i >= g_k; the gap is monotone, so j can never beat k and is dominated. Dropping dominated candidates preserves the true maximum. The segment tree stores exactly the per-coordinate maxima (and, in bubbles, one best per branch filtered by gidx and max_gap), so each query is exact. Time O(n log n), memory O(n) per read. ∎
Delayed insertion (r_j + k <= r_i) removes only read-overlap-infeasible transitions; branch exclusion removes only branch-switching transitions. Both are exact constraints.
test_pangenome_chain.py#!/usr/bin/env python3
"""Verification suite for pangenome_chain.py.
Run: python3 -m pytest -q test_pangenome_chain.py
or: python3 test_pangenome_chain.py
"""
from __future__ import annotations
import random
import sys
import time
from pangenome_chain import (
Anchor,
ColoredVariationGraph,
Node,
build_anchors,
chain_brute,
chain_one_orientation,
chain_read,
minimizers,
minimizers_naive,
proof_domination,
revcomp,
stream_reads,
)
COL = lambda *ids: sum(1 << i for i in ids)
# ---------------------------------------------------------------------------
# 1. minimizer equivalence
# ---------------------------------------------------------------------------
def test_minimizers_match_naive():
rng = random.Random(1)
for _ in range(200):
n = rng.randint(1, 60)
seq = "".join(rng.choice("ACGT") for _ in range(n))
k = rng.randint(1, 7)
w = rng.randint(1, 8)
got = [(p, c, s) for p, c, s, _ in minimizers(seq, k, w)]
exp = [(p, c, s) for p, c, s, _ in minimizers_naive(seq, k, w)]
assert got == exp, (seq, k, w, got, exp)
print("ok minimizers(deque) == minimizers(naive)")
# ---------------------------------------------------------------------------
# 2. fast DP == brute force on random anchors
# ---------------------------------------------------------------------------
def _random_anchors(rng, n, *, dup_ridx=True, bubbles=None):
bubbles = bubbles or {}
out = []
for i in range(n):
ridx = rng.randint(0, 50 if dup_ridx else n)
gidx = rng.randint(0, 200)
# Bubble tags are only valid for anchors whose graph coordinate is
# geometrically inside the bubble span (this is guaranteed for real
# graphs by detect_bubbles).
node_b, branch_b = -1, -1
for bid, (lo, hi) in bubbles.items():
if lo <= gidx <= hi:
node_b, branch_b = bid, rng.randint(0, 1)
break
out.append(
Anchor(
read_id=0,
ridx=ridx,
gidx=gidx,
strand=1,
weight=1.0,
colors=rng.randint(0, 15),
node=0,
pos=0,
bubble=node_b,
branch=branch_b,
)
)
return out
def test_fast_equals_brute_random():
rng = random.Random(7)
for trial in range(300):
n = rng.randint(1, 60)
k = rng.randint(1, 5)
lam = rng.choice([0.0, 0.01, 0.1, 1.0])
max_gap = rng.choice([10, 100, 100000])
bubbles = {0: (40, 120)} if rng.random() < 0.5 else {}
anchors = _random_anchors(rng, n, bubbles=bubbles)
fast = chain_one_orientation(anchors, k, lam, max_gap, bubbles)
slow = chain_brute(anchors, k, lam, max_gap, bubbles)
assert abs(fast.score - slow.score) < 1e-9, (
trial, fast.score, slow.score, k, lam, max_gap, bubbles
)
assert fast.colors == slow.colors or fast.score == 0.0
print("ok fast DP == O(n^2) brute force (300 random trials)")
# ---------------------------------------------------------------------------
# 3. read-overlap prevents double counting repeated sequence
# ---------------------------------------------------------------------------
def test_no_double_count_same_read_interval():
# 6 anchors all at the same read position, mapping to 6 distinct graph
# positions. A correct chain may use at most one of them.
anchors = [
Anchor(0, 10, 100 + 20 * i, 1, 1.0, COL(1), 0, 10, -1, -1)
for i in range(6)
]
res = chain_one_orientation(anchors, k=15, lam=0.01, max_gap=10**9, bubbles={})
assert len(res.chain) == 1, [a.gidx for a in res.chain]
print("ok repeated read interval counted at most once")
def test_tandem_repeat_not_double_counted():
# Read has one copy of a repeat; both graph copies are present. The chain
# must pick a single consistent copy, not mix them.
anchors = []
for i in range(8):
r = 5 * i
anchors.append(Anchor(0, r, 200 + r, 1, 1.0, COL(1), 0, r, -1, -1)) # copy 1
anchors.append(Anchor(0, r, 2000 + r, 1, 1.0, COL(1), 0, r, -1, -1)) # copy 2
res = chain_one_orientation(anchors, k=4, lam=0.01, max_gap=10**9, bubbles={})
# graph coordinates of chosen chain must all lie in one copy
gs = [a.gidx for a in res.chain]
assert all(g >= 2000 for g in gs) or all(g < 2000 for g in gs), gs
print("ok tandem-repeat chain stays on one graph copy")
# ---------------------------------------------------------------------------
# 4. reverse-complement read orientation
# ---------------------------------------------------------------------------
def test_reverse_complement_read():
rng = random.Random(3)
seq = "".join(rng.choice("ACGT") for _ in range(400))
g = ColoredVariationGraph({1: Node(seq, COL(0))}, [])
g.linearize()
idx, _ = g.index_minimizers(k=15, w=10)
fwd = chain_read(0, seq, g, idx, k=15, w=10, lam=0.01, max_gap=10**9)
rc = chain_read(1, revcomp(seq), g, idx, k=15, w=10, lam=0.01, max_gap=10**9)
assert fwd.chain, "forward chain must be non-empty"
assert rc.chain, "reverse chain must be non-empty"
assert abs(fwd.score - rc.score) < 1e-9, (fwd.score, rc.score)
assert all(a.strand == -1 for a in rc.chain), [a.strand for a in rc.chain]
assert all(a.strand == 1 for a in fwd.chain)
print(f"ok revcomp read chains symmetrically (score={fwd.score:.3f})")
# ---------------------------------------------------------------------------
# 5. bubble branch switching is forbidden
# ---------------------------------------------------------------------------
def test_bubble_no_branch_switch():
# bubble 0 spans gidx [100, 200); branch 0 at gidx 110..130, branch 1 at
# 160..180. A chain must not use both branches.
anchors = [
Anchor(0, 0, 50, 1, 1.0, COL(1), 0, 0, -1, -1), # before
Anchor(0, 5, 110, 1, 1.0, COL(1), 0, 5, 0, 0), # branch 0
Anchor(0, 10, 120, 1, 1.0, COL(1), 0, 10, 0, 0),
Anchor(0, 15, 160, 1, 1.0, COL(1), 0, 15, 0, 1), # branch 1
Anchor(0, 20, 170, 1, 1.0, COL(1), 0, 20, 0, 1),
Anchor(0, 30, 250, 1, 1.0, COL(1), 0, 30, -1, -1), # after
]
bubbles = {0: (100, 200)}
res = chain_one_orientation(anchors, k=1, lam=0.0, max_gap=10**9, bubbles=bubbles)
branches = {a.branch for a in res.chain if a.bubble == 0}
assert len(branches) <= 1, [a.gidx for a in res.chain]
print("ok bubble chain does not switch branches")
# ---------------------------------------------------------------------------
# 6. provenance is preserved (intersection of sample colours)
# ---------------------------------------------------------------------------
def test_provenance_intersection():
# sample masks: 0b0001, 0b0011, 0b0111
anchors = [
Anchor(0, 0, 0, 1, 1.0, COL(0, 1, 2), 0, 0, -1, -1),
Anchor(0, 10, 10, 1, 1.0, COL(0, 1), 0, 10, -1, -1),
Anchor(0, 20, 20, 1, 1.0, COL(0), 0, 20, -1, -1),
]
res = chain_one_orientation(anchors, k=1, lam=0.0, max_gap=10**9, bubbles={})
assert len(res.chain) == 3
assert res.colors == COL(0), bin(res.colors)
print("ok chain provenance = intersection of anchor colours")
# ---------------------------------------------------------------------------
# 7. end-to-end graph / sequence test
# ---------------------------------------------------------------------------
def test_end_to_end_colored_graph():
rng = random.Random(11)
prefix = "".join(rng.choice("ACGT") for _ in range(300))
allele_a = "".join(rng.choice("ACGT") for _ in range(300))
allele_b = "".join(rng.choice("ACGT") for _ in range(300))
suffix = "".join(rng.choice("ACGT") for _ in range(300))
# proper diamond bubble: 0 -> 3 -> 1 (sample 0) and 0 -> 2 -> 1 (sample 1)
g = ColoredVariationGraph(
{
0: Node(prefix, COL(0, 1)),
2: Node(allele_b, COL(1)),
3: Node(allele_a, COL(0)),
1: Node(suffix, COL(0, 1)),
},
[(0, 3), (0, 2), (3, 1), (2, 1)],
)
g.linearize()
index, dropped = g.index_minimizers(k=15, w=10)
bubbles, nb = g.detect_bubbles()
assert bubbles, "diamond bubble should be detected"
read = prefix + allele_a + suffix # follows sample 0 allele
res = chain_read(
0, read, g, index, k=15, w=10, lam=0.01, max_gap=10**9,
node_bubble=nb, bubbles=bubbles,
)
assert res.chain
nodes = {an.node for an in res.chain}
assert 2 not in nodes, "read must not use the other allele"
assert nodes <= {0, 1, 3}, nodes
print(f"ok end-to-end graph chain nodes={sorted(nodes)} score={res.score:.3f}")
# ---------------------------------------------------------------------------
# 8. sub-quadratic scaling
# ---------------------------------------------------------------------------
def test_subquadratic_scaling():
rng = random.Random(5)
times = {}
for n in (1000, 2000, 4000, 8000):
anchors = [
Anchor(0, rng.randint(0, 10**6), rng.randint(0, 10**6), 1, 1.0,
COL(0), 0, 0, -1, -1)
for _ in range(n)
]
t0 = time.perf_counter()
chain_one_orientation(anchors, k=15, lam=0.001, max_gap=10**9, bubbles={})
times[n] = time.perf_counter() - t0
ratio = times[8000] / max(times[1000], 1e-9)
# O(n^2) would give ratio 64; O(n log n) gives about 8*log ratio ~ 9-10.
assert ratio < 20, f"scaling ratio {ratio:.1f} looks quadratic: {times}"
print(f"ok scaling near n log n: ratio(8k/1k)={ratio:.1f} times={ {k: round(v,4) for k,v in times.items()} }")
# ---------------------------------------------------------------------------
# 9. memory-capped streaming driver
# ---------------------------------------------------------------------------
def test_stream_reads():
rng = random.Random(13)
seq = "".join(rng.choice("ACGT") for _ in range(200))
g = ColoredVariationGraph({1: Node(seq, COL(0))}, [])
reads = [seq, revcomp(seq), seq[:100]]
got = list(stream_reads(reads, g, k=15, w=10, lam=0.01, max_gap=10**9))
assert len(got) == 3
assert all(res.chain for _, res in got)
print("ok stream_reads yields one result per read under memory cap")
def main():
test_minimizers_match_naive()
test_fast_equals_brute_random()
test_no_double_count_same_read_interval()
test_tandem_repeat_not_double_counted()
test_reverse_complement_read()
test_bubble_no_branch_switch()
test_provenance_intersection()
test_end_to_end_colored_graph()
test_subquadratic_scaling()
test_stream_reads()
print()
print(proof_domination())
if __name__ == "__main__":
main()
cd ~
python3 -m pytest -q test_pangenome_chain.py # 10 passed
python3 test_pangenome_chain.py # verbose, prints evidence
ok minimizers(deque) == minimizers(naive)
ok fast DP == O(n^2) brute force (300 random trials)
ok repeated read interval counted at most once
ok tandem-repeat chain stays on one graph copy
ok revcomp read chains symmetrically (score=14.760)
ok bubble chain does not switch branches
ok chain provenance = intersection of anchor colours
ok end-to-end graph chain nodes=[0, 1, 3] score=25.760
ok scaling near n log n: ratio(8k/1k)=11.5 times={1000: 0.0034, 2000: 0.0077, 4000: 0.0172, 8000: 0.0392}
ok stream_reads yields one result per read under memory cap
10 passed in 0.15s
k, lambda, max_gap, checking both optimal score and provenance bitmask: mismatches: 0 / 5000.1000 -> 8000 gives ratio 11.5 (0.0034 -> 0.0392); quadratic would be 64. Matches n log n.strand == -1.import random
from pangenome_chain import Anchor, chain_one_orientation, chain_brute
rng = random.Random(12345)
bad = 0
for trial in range(5000):
nb = rng.randint(0, 3)
bubbles, x = {}, 0
for b in range(nb):
lo = x + rng.randint(0, 30); hi = lo + rng.randint(1, 60)
bubbles[b] = (lo, hi); x = hi + rng.randint(1, 20)
n = rng.randint(1, 80); k = rng.randint(1, 8)
lam = rng.choice([0.0, 0.001, 0.01, 0.1])
max_gap = rng.choice([1, 5, 30, 1000, 10**9])
anchors = []
for i in range(n):
ridx = rng.randint(0, 120); gidx = rng.randint(0, x + 40)
b = br = -1
for bid, (lo, hi) in bubbles.items():
if lo <= gidx <= hi: b, br = bid, rng.randint(0, 1); break
anchors.append(Anchor(0, ridx, gidx, 1, 1.0, rng.randint(0, 31), 0, 0, b, br))
fast = chain_one_orientation(anchors, k, lam, max_gap, bubbles)
slow = chain_brute(anchors, k, lam, max_gap, bubbles)
if abs(fast.score - slow.score) > 1e-9 or fast.colors != slow.colors:
bad += 1
print("mismatches:", bad, "/ 5000")
| Quantity | Cost |
|---|---|
| Minimizers per read | O(L) with a monotonic deque |
| Graph index | O(#graph minimizers); keys with occ > max_occ dropped |
| Chaining per read/strand | O(n log n) time, O(n) auxiliary |
| Bubble handling | O(n * m), m = anchors per branch of the bounded bubble |
| Provenance | one big-int AND per chosen anchor |
| Peak memory | graph index + one read's anchors (streaming) |
from pangenome_chain import Node, ColoredVariationGraph, stream_reads, revcomp
COL = lambda *ids: sum(1 << i for i in ids) # sample ids -> bitmask
g = ColoredVariationGraph(
{0: Node("ACGT...", COL(0, 1)), 1: Node("TTGA...", COL(1))},
[(0, 1)],
)
for read_id, result in stream_reads(open_reads(), g, k=15, w=10, lam=0.01, max_gap=10000):
# result.score, result.chain, result.colors
write_chain(read_id, result)
Because stream_reads yields lazily, the caller decides whether to write each chain to disk (true memory cap) or retain them.
gidx (otherwise it chains backwards inside a bubble), and (2) the same-branch scan must enforce max_gap. Both are now covered by the randomized brute-force comparison.|dr-dg|-style gap, the same dominance framework generalizes to a 2D (CDQ / Fenwick-of-Fenwick) structure, but the read-overlap and bubble invariants are unchanged.# Evidence - Problem class: python-20261012-pangenome-minimizer-chaining-with-structural-variants - Model: openrouter/deepseek/deepseek-v4.1-flash - Solved: 2026-10-12T16:14:48.568Z - Verification: solution produced by pi in sandbox; see signatures.json
{"description": "Given millions of noisy long reads, a colored variation graph, and a memory cap, compute a maximum-scoring chain of minimizer anchors that may cross bounded structural-variant bubbles without double-counting repeated sequence. The solution must preserve sample provenance, handle reverse-complement paths, and provide a subquadratic algorithm with a correctness argument for its pruning strategy.", "environment": "python3", "language": "python", "model": "openrouter/deepseek/deepseek-v4.1-flash", "problem_class": "python-20261012-pangenome-minimizer-chaining-with-structural-variants", "provider": "openrouter", "solved_at": "2026-10-12T16:14:48.568Z", "version": "3.11"}